2018-06-23 13:04 BST

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0002980OpenFOAM[All Projects] Bugpublic2018-06-14 11:32
Reportercvr 
Assigned Towill 
PrioritynormalSeveritymajorReproducibilityN/A
StatusresolvedResolutionfixed 
PlatformUnixOSGNU/Linux DebianOS Version9
Product Versiondev 
Target VersionFixed 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:

(cmptMag(g.value())/mag(g.value()))*hRe

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:

(-g.value()/mag(g.value()))*hRe

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

-Relationships
+Relationships

-Notes

~0009754

will (manager)

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?

~0009755

cvr (reporter)

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.

~0009756

henry (manager)

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.

~0009760

cvr (reporter)

Hello,

"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.

~0009761

henry (manager)

> 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.

~0009762

cvr (reporter)

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

~0009763

will (manager)

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.
+Notes

-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
+Issue History