View Issue Details

IDProjectCategoryView StatusLast Update
0002980OpenFOAMBugpublic2018-06-14 11:32
Reportercvr Assigned Towill  
Status resolvedResolutionfixed 
PlatformUnixOSGNU/Linux DebianOS Version9
Product Versiondev 
Fixed in Versiondev 
Summary0002980: Possible bug in src/finiteVolume/cfdTools/general/include/gh.H
DescriptionIn file src/finiteVolume/cfdTools/general/include/gh.H, a geopotential 'gh' is computed together with its reference 'ghRef'.
If the magnitude of 'g' is > 0 (or small), the code for 'ghRef' becomes:


meaning 'G ยท A href' (where upper case are vectors) where versor 'A' is given by 'cmptMag(G)/mag(G)'. If 'G' is pointing downwards (in the -y or -z direction), 'A' will point upwards in the opposed direction.

For the sake of simplicity I will refer only to an x-y plane. If 'G=(-5,-10)', i.e. a vector in the 3rd quadrant, then 'A' will again be in its opposite direction in the 1st quadrant. However if 'G=(-5,+10)' or 'G=(+5,-10)', meaning 2nd and 4th quadrants, A will not be on the opposite direction, as it will remain in the 1st quadrant. If 'G' is in the 1st quadrant like 'G=(+5,+10)', 'G' and 'A' will be colinear and have equal directions.

The conclusions are the following:
1. If I perform a simulation with 'hRef = 30 m' (some value different than 0) and 'g = (0,-9.8,0)', then 'ghRef = -294 m^2/s^2'. The geopotential at an height of 'y' of 0 m is 'gh = 0 - gHref = 294 m^2/s^2'.

2. If I flip the simulation around 0, having 'hRef = -30 m' and 'g = (0,+9.8,0)', I obtain 'ghRef = -294 m^2/s^2' again. For the same y=0 surface I obtain a geopotential 'gh = 294 m^2/s^2', instead of 'gh = -294 m^2/s^2' which I was expecting.

3. If I rotate my referential by 30 degrees keeping 'hRef = 30 m', I obtain 'ghRef = -294 m^2/s^2'.

4. If I rotate my referential by -30 degrees keeping 'hRef = 30 m', I obtain 'ghRef = -147 m^2/s^2'.
Steps To ReproduceN/A
Additional InformationI would propose changing the code to:


thus '-g' instead of 'cmptMag(g)'.
TagsNo tags attached.



2018-06-13 16:56

manager   ~0009754

g & g reduces to mag(g)*mag(g), so presumably we could simplify the entire gh.H file to the following:

    Info<< "Calculating field g.h\n" << endl;
    dimensionedScalar ghRef(- mag(g)*hRef);
    volScalarField gh("gh", (g & mesh.C()) - ghRef);
    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);

Could you confirm if this works correctly with your test?


2018-06-13 17:24

reporter   ~0009755

Yes, your formulation is better as is even simpler. The value of 'g & -g/mag(g)' will always be negative because the vectors are colinear but with opposite directions.

Now I believe both of these formulations are what was intended in the 'ghRef' calculation. I would like to see some discussion on this to understand what formulation is indeed correct (the old or the proposed one) and, if it is the old one, why is that so.

I also presume this has not produced that many problems as the hRef is usually undefined (e.g. on most tutorials) and in those cases its value becomes 0.


2018-06-13 17:41

manager   ~0009756

The issue is not what is correct or incorrect but what is a convenient way to define hRef; any of the forms would work given an appropriate value for hRef. I chose to define hRef in the positive quadrant which was fine for the cases we have used it for but I see it is not so convenient if you want to transform the case without changing hRef and get the same results.


2018-06-14 10:28

reporter   ~0009760


"The issue is not what is correct or incorrect but what is a convenient way to define hRef"

Sorry but I disagree. I believe the intent in computing ghRef and gh is:
1. To compute a geopotential.
2. Have simulation results which are invariant of the referential choosen.

To me, the present formulation in gh.H appears to fail in point 2. One may argue that its innocuous in terms of the final result, because ghRef is merely a reference value and gh is still correct (i.e. in 'g & mesh.C() - ghRef' the 'g & mesh.C()' term is correct, only ghRef is incorrectly shifted). However conceptually if you want a gh field whose 0 value respects to a particular iso-surface, you will have to correct that after the simulation.

This will not affect most people because I presume one does not even bother in setting constant/hRef, so that ghRef will become 0. But for people developing a solver that needs a specific ghRef pegged to a particular surface (e.g. oceanography or meteorology), having to account for an arbitary shift in gh is undesirable.


2018-06-14 10:41

manager   ~0009761

> Sorry but I disagree. I believe the intent in computing ghRef and gh is:
> 1. To compute a geopotential.
> 2. Have simulation results which are invariant of the referential choosen.

That was not the intent, it was purely numerical to avoid the large pressure offset caused by a domain being far from the origin.


2018-06-14 11:23

reporter   ~0009762

Thank you for clarifying this. Putting it in such way its purpose is merely numerical.

Either way, if you want to give it another physical meaning (again, perhaps only useful for meteorology and oceanography people) the formulation I proposed doesn't seem unreasonable.

Best regards


2018-06-14 11:32

manager   ~0009763

I think we are in agreement that the new form is a more natural definition, even if it is fundamentally just a matter of convention. The change can also be made without too much disruption. Cases without hRef (the majority) are unaffected, and cases in which g is aligned with a negative axis (again, the majority) are also unaffected, as -g == cmptMag(g).

The change has been made in dev as commit 268f1f61.

Issue History

Date Modified Username Field Change
2018-06-13 15:51 cvr New Issue
2018-06-13 16:56 will Note Added: 0009754
2018-06-13 17:24 cvr Note Added: 0009755
2018-06-13 17:41 henry Note Added: 0009756
2018-06-14 10:28 cvr Note Added: 0009760
2018-06-14 10:41 henry Note Added: 0009761
2018-06-14 11:23 cvr Note Added: 0009762
2018-06-14 11:32 will Assigned To => will
2018-06-14 11:32 will Status new => resolved
2018-06-14 11:32 will Resolution open => fixed
2018-06-14 11:32 will Fixed in Version => dev
2018-06-14 11:32 will Note Added: 0009763