View Issue Details
| ID | Project | Category | View Status | Date Submitted | Last Update |
|---|---|---|---|---|---|
| 0002663 | OpenFOAM | Bug | public | 2017-08-11 11:10 | 2017-08-16 11:18 |
| Reporter | volker1 | Assigned To | henry | ||
| Priority | normal | Severity | crash | Reproducibility | always |
| Status | resolved | Resolution | fixed | ||
| OS | Ubuntu 16.04.3 LTS | ||||
| Summary | 0002663: Possible bug in chtMultiRegionFoam (chtMultiRegionFoam/fluid/pEqn.H) of V5.0 | ||||
| Description | I 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. Thanks! | ||||
| Tags | No tags attached. | ||||
|
|
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? |
|
|
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. |
|
|
Hy Henry, thanks for taking care. I confirm that your commit proposed at https://github.com/OpenFOAM/OpenFOAM-dev/commit/f2746b55362d7876c8ba50d3dd9a73241ba3a832 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. |
|
|
- 2/3 part of the zipfile |
|
|
- 3/3 part of the zipfile |
|
|
|
|
|
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. |
|
|
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. |
|
|
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! |
|
|
> 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 https://bugs.openfoam.org/view.php?id=2512 |
| 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: example_chTMRF_1.zip | |
| 2017-08-15 09:45 | volker1 | Note Added: 0008585 | |
| 2017-08-15 09:46 | volker1 | File Added: example_chTMRF_2.zip | |
| 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: example_chTMRF_3.zip | |
| 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 |
