View Issue Details
|ID||Project||Category||View Status||Date Submitted||Last Update|
|0002980||OpenFOAM||[All Projects] Bug||public||2018-06-13 15:51||2018-06-14 11:32|
|Platform||Unix||OS||GNU/Linux Debian||OS Version||9|
|Fixed in Version||dev|
|Summary||0002980: Possible bug in src/finiteVolume/cfdTools/general/include/gh.H|
|Description||In 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 Reproduce||N/A|
|Additional Information||I would propose changing the code to:|
thus '-g' instead of 'cmptMag(g)'.
|Tags||No tags attached.|
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?
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.
||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.|
"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.
> 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.
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.
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.
|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|