From 0e0781bf004e572390d96693d65a5cc906ec1e91 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 4 Jun 2025 15:48:55 +0200 Subject: [PATCH 1/2] fix the energetic of ShearingBox+Fargo The work done by the gravity force was not properly removed when using fargo --- src/fluid/addSourceTerms.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/fluid/addSourceTerms.hpp b/src/fluid/addSourceTerms.hpp index a0cb5a71..6451c99a 100644 --- a/src/fluid/addSourceTerms.hpp +++ b/src/fluid/addSourceTerms.hpp @@ -85,6 +85,9 @@ struct Fluid_AddSourceTermsFunctor { } if(haveFargo) { Uc(MX1,k,j,i) += TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * sbS * x1(i); + if constexpr(Phys::pressure) { + Uc(ENG,k,j,i) += TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * sbS * x1(i) * Vc(VX1,k,j,i); + } } #endif // fetch fargo velocity when required From 688241a57f0db12e672b971f3dd4b8b5a42257e3 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 11 Jun 2025 17:06:53 +0200 Subject: [PATCH 2/2] - compensate for body force straight into RHS instead of removing what was added in addSourceTerms - fix reference file in Shearing box tests following this bug fix. --- reference | 2 +- src/fluid/addSourceTerms.hpp | 6 ------ src/fluid/calcRightHandSide.hpp | 12 +++++++++--- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/reference b/reference index c3289420..c4082b99 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit c328942083b618a85b84eb16ce9fd35e67c00597 +Subproject commit c4082b99a4c542def3177c96cb35b1c9d9002f18 diff --git a/src/fluid/addSourceTerms.hpp b/src/fluid/addSourceTerms.hpp index 6451c99a..e3f74586 100644 --- a/src/fluid/addSourceTerms.hpp +++ b/src/fluid/addSourceTerms.hpp @@ -83,12 +83,6 @@ struct Fluid_AddSourceTermsFunctor { Uc(MX1,k,j,i) += TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * Vc(VX2,k,j,i); Uc(MX2,k,j,i) += - TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * Vc(VX1,k,j,i); } - if(haveFargo) { - Uc(MX1,k,j,i) += TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * sbS * x1(i); - if constexpr(Phys::pressure) { - Uc(ENG,k,j,i) += TWO_F * dt * Vc(RHO,k,j,i) * OmegaZ * sbS * x1(i) * Vc(VX1,k,j,i); - } - } #endif // fetch fargo velocity when required [[maybe_unused]] real fargoV = ZERO_F; diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp index 38ef76d5..7777d1c2 100644 --- a/src/fluid/calcRightHandSide.hpp +++ b/src/fluid/calcRightHandSide.hpp @@ -472,12 +472,18 @@ struct Fluid_CalcRHSFunctor { // Body force if(needBodyForce) { - rhs[MX1+dir] += dt * Vc(RHO,k,j,i) * bodyForce(dir,k,j,i); + real bf = bodyForce(dir,k,j,i); + #if GEOMETRY == CARTESIAN + if(haveFargo && dir==IDIR) { + // Remove Body force that is already compensated by Fargo shear*Rotation + bf -= -2*Omega*sbS * x1(i); + } + #endif + rhs[MX1+dir] += dt * Vc(RHO,k,j,i) * bf; if constexpr(Phys::pressure) { // rho * v . f, where rhov is taken as a volume average of Flux(RHO) rhs[ENG] += HALF_F * dtdV * dl * - (Flux(RHO,k,j,i) + Flux(RHO, k+koffset, j+joffset, i+ioffset)) * - bodyForce(dir,k,j,i); + (Flux(RHO,k,j,i) + Flux(RHO, k+koffset, j+joffset, i+ioffset)) * bf; } // Pressure // Particular cases if we do not sweep all of the components