View Issue Details

IDProjectCategoryView StatusLast Update
0002663OpenFOAMBugpublic2017-08-16 11:18
Reportervolker1 Assigned Tohenry  
Status resolvedResolutionfixed 
OSUbuntu 16.04.3 LTS 
Summary0002663: Possible bug in chtMultiRegionFoam (chtMultiRegionFoam/fluid/pEqn.H) of V5.0
DescriptionI recently upgraded from 4.1 to 5.0 and stumbled into an issue related to chtMultiRegionFoam.
As usual when it comes to an upgrade, I verified that I get the same results from 4.1 and 5.0 for the "multiRegionHeater" tutorial and a few simple own test cases.
In a more complex flow simulation I noticed that the flow field behaved unstable and the calculation crashed quite early with chtMultiRegionFoam of 5.0. However the same simulation was no problem with 4.1 before.
One cannot fix the problem by using the PBiCGStab solver instead of PBiCG (which had been suggested by an error message in one case).

Making hybrid solvers with files from both versions I was able to track my problem down to a modification in $solverDir/chtMultiRegionFoam/fluid/pEqn.H.
There the thermodynamic density update calculation now uses "p" (5.0) instead of "p_rgh" (4.1 and before).
Using "p_rgh" again in 5.0 (in lines 37 and 74 of pEqn.H) stabilized the flow field and chtMultiRegionFoam did not crash any more. Anyway the flow field still looked more unstable for 5.0 than for 4.1 at the early time steps but both U fields seem to become close to equal with time as it should be (it should become a steady state flow field in this case and gravity is zero).
As also the formalism of the density calculation changed a bit, I tend to suppose that the variable change p_rgh -> p might have been an unintended by-product of the intended formal change.
Nevertheless I did not go into the equations that deep to decide which one really is the proper variable ( I would of course prefer "p_rgh" :-) ).

If this should not just be a bug, but the change of the pressure variable at this place in fact should have had a physical reason, one should maybe try to find a more stable implementation.

TagsNo tags attached.



2017-08-11 12:15

manager   ~0008560

Using either p or p_rgh to update density are both valid but in principle p is better as it includes any limiters or other post-solution manipulations of p. However, there is an ordering issue with the current implementation.

Could you provide a test-case which demonstrates the problem with the current implementation? Can it be reproduced with any of the tutorial cases?


2017-08-11 14:29

manager   ~0008561

I have just pushed an update to OpenFOAM-dev:

commit f2746b55362d7876c8ba50d3dd9a73241ba3a832

which should resolve the issue. Could you test it on your case and if it runs OK I will make the same change to 5.x.


2017-08-15 09:45

reporter   ~0008585

Hy Henry,

thanks for taking care.

I confirm that your commit proposed at
fixes my crash problem in my test case.

The bug seemed not to have any severe effect the "multiRegionHeater" tutorial case or any other of my more complex test cases.
"U" differs a bit in low decimal places, thats it.

I am sorry that the zipfile of my example calculation case to demonstrate the bug is about 12MB.
(The mesh to this case had been generated by snappyHexMesh from a CAD file via *.stl files.)
However file upload here seems limited to 5MB, so I used zipsplit to make three files. I hope this is all right.
Even with the crash now avoided with V5.0, the following observations remained for my example:

1. Having a look at the velocity field "U" in the x-y plane view in paraview on the region "fluid_parts", "U" looks and develops much more stable for the first time steps if one uses chtMultiRegionFoam of V4.1.
In the same view with (fixed) version 5.0 intermediate "eddy states" occur.
For version 2.2.2 the same x-y plane view shows a similar development of U like for 4.1 without "eddy states".

2. For the V4.1 version also the time step becomes larger (1.4E-5 to 1.E-5) and the courant number reaches the prescribed value of 0.45 while it remains about 0.23 for (the already fixed) V5.0.
I did not check that out yet, it might be related to the new residual control feature.

In summary the fixed V5 seems to work, too, but uses about 50% more computation time for this case.

For these reasons the V5.0 version still seems suspicious to me beyond this fixed bug. (4,718,395 bytes)


2017-08-15 09:46

reporter   ~0008586

- 2/3 part of the zipfile (4,781,313 bytes)


2017-08-15 09:47

reporter   ~0008587

- 3/3 part of the zipfile


2017-08-15 09:51

reporter (3,390,164 bytes)


2017-08-15 10:03

manager   ~0008588

I have applied the change to pEqn.H to OpenFOAM-5.x.

Please let me know what changes between OpenFOAM-4.1 and 5.x/dev are causing the differences in the velocity field.


2017-08-15 11:16

reporter   ~0008589

A test with a fixed time step (1.E-5s) for both versions indicated that the physics equations of 4.1 and 5.x/dev do almost digitally the same for "U". A few low decimal places in "U" differ but field transient development is basically the same including that I see temporary "eddy states" now for "U" with V4.1, too.
Just for completeness: the original 5.0 version crashes early with this fixed time step, too.
This means basically the time step seems responsible for the remaining differences in "U", not the physics equations.

As already mentioned, adaptive time step control by Courant number (e.g. maxCo = 0.45;) results in different time steps for both versions in my case.

So I would think one has to look into the time step control.


2017-08-16 11:05

reporter   ~0008594

I had a look at the time step control.

In line 39 of chtMultiregionFoam/solid/solidRegionDiffNo.C the field "kapparhoCpbyDelta" is used to evaluate the Diffusion number.

In V4.1 the (inverse) mesh dimensions enter this value (as "(mesh.surfaceInterpolation::deltaCoeffs())" ) while in V5.0
the square of the (inverse) mesh dimensions ( sqr(mesh.surfaceInterpolation::deltaCoeffs()) is used.

The dimensions of kappa/Cp/rho are "m²/s". The time step multiplies "s" so we still need "1/m²" to get a dimensionless number.
So I think it is V5.0 that calculates the Diffusion number properly.

In my isothermal case I did not feel a need to prescribe the maximum diffusion number "maxDi" explicitely. Consequently the default value "maxDi=10;" as prescribed in line 32 of chtMultiregionFoam/solid/readSolidTimeControls.H was used in all my calculations.
As the diffusion number calculation itself however differed by the above "sqr()" evaluation I caught the observed time step reduction by the diffusion number.

One might consider if "maxDi=10" is still a proper default value for the thermal diffusion number but for my example the V41-V50 differences seem clarified now.

Thanks again!


2017-08-16 11:17

manager   ~0008595

> So I think it is V5.0 that calculates the Diffusion number properly.

Yes, see

commit bda41c441c60f27698608130c2b190881e7e866e

    Diffusion number: Corrected in chtMultiRegionFoam and pyrolysisModels::reactingOneDim
    Resolves bug-report

Issue History

Date Modified Username Field Change
2017-08-11 11:10 volker1 New Issue
2017-08-11 12:04 henry Category Contribution => Bug
2017-08-11 12:15 henry Note Added: 0008560
2017-08-11 14:29 henry Note Added: 0008561
2017-08-15 09:45 volker1 File Added:
2017-08-15 09:45 volker1 Note Added: 0008585
2017-08-15 09:46 volker1 File Added:
2017-08-15 09:46 volker1 Note Added: 0008586
2017-08-15 09:47 volker1 Note Added: 0008587
2017-08-15 09:51 volker1 File Added:
2017-08-15 10:03 henry Note Added: 0008588
2017-08-15 11:16 volker1 Note Added: 0008589
2017-08-16 11:05 volker1 Note Added: 0008594
2017-08-16 11:17 henry Note Added: 0008595
2017-08-16 11:18 henry Assigned To => henry
2017-08-16 11:18 henry Status new => resolved
2017-08-16 11:18 henry Resolution open => fixed
2017-08-16 11:18 henry Fixed in Version => 5.x