From 2d857c4e9b5ac70ef71c054bf76aadb86ff4799d Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 22 Apr 2025 14:43:32 +0200 Subject: [PATCH 1/4] keep Flux array in each direction (increase memory footprint, but easier post-treatment) --- src/fluid/RiemannSolver/riemannSolver.hpp | 4 ---- src/fluid/addNonIdealMHDFlux.hpp | 3 +-- src/fluid/calcParabolicFlux.hpp | 8 ++++---- src/fluid/calcRightHandSide.hpp | 4 ++-- src/fluid/evolveStage.hpp | 6 +++--- src/fluid/fluid.hpp | 12 ++++++++---- src/rkl/rkl.hpp | 12 +++++------- 7 files changed, 23 insertions(+), 26 deletions(-) diff --git a/src/fluid/RiemannSolver/riemannSolver.hpp b/src/fluid/RiemannSolver/riemannSolver.hpp index abb70eed5..b6c1e6f11 100644 --- a/src/fluid/RiemannSolver/riemannSolver.hpp +++ b/src/fluid/RiemannSolver/riemannSolver.hpp @@ -70,7 +70,6 @@ class RiemannSolver { IdefixArray4D Vc; IdefixArray4D Vs; - IdefixArray4D Flux; IdefixArray3D cMax; Fluid* hydro; DataBlock *data; @@ -90,7 +89,6 @@ class RiemannSolver { template RiemannSolver::RiemannSolver(Input &input, Fluid* hydro) : Vc{hydro->Vc}, Vs{hydro->Vs}, - Flux{hydro->FluxRiemann}, cMax{hydro->cMax}, hydro{hydro}, data{hydro->data} @@ -148,8 +146,6 @@ RiemannSolver::RiemannSolver(Input &input, Fluid* hydro) : Vc{hydro- mySolver = HLL_DUST; } - - // Shock flattening this->haveShockFlattening = input.CheckEntry(std::string(Phys::prefix),"shockFlattening")>=0; // Init shock flattening diff --git a/src/fluid/addNonIdealMHDFlux.hpp b/src/fluid/addNonIdealMHDFlux.hpp index 52c281d2d..9ae6af304 100644 --- a/src/fluid/addNonIdealMHDFlux.hpp +++ b/src/fluid/addNonIdealMHDFlux.hpp @@ -19,8 +19,7 @@ void Fluid::AddNonIdealMHDFlux(const real t) { if constexpr(Phys::mhd) { int ioffset,joffset,koffset; - - IdefixArray4D Flux = this->FluxRiemann; + IdefixArray4D Flux = this->FluxRiemann[dir]; IdefixArray4D Vc = this->Vc; IdefixArray4D Vs = this->Vs; IdefixArray3D dMax = this->dMax; diff --git a/src/fluid/calcParabolicFlux.hpp b/src/fluid/calcParabolicFlux.hpp index a8c6761c0..8713f9107 100644 --- a/src/fluid/calcParabolicFlux.hpp +++ b/src/fluid/calcParabolicFlux.hpp @@ -42,7 +42,7 @@ void Fluid::CalcParabolicFlux(const real t) { if(data->haveFargo && viscosityStatus.isExplicit) { data->fargo->AddVelocityFluid(t,this); } - this->viscosity->AddViscousFlux(dir,t, this->FluxRiemann); + this->viscosity->AddViscousFlux(dir,t, this->FluxRiemann[dir]); // Remove back Fargo velocity if(data->haveFargo && viscosityStatus.isExplicit) { @@ -53,18 +53,18 @@ void Fluid::CalcParabolicFlux(const real t) { // Add thermal diffusion if( (thermalDiffusionStatus.isExplicit && (!data->rklCycle)) || (thermalDiffusionStatus.isRKL && data->rklCycle)) { - this->thermalDiffusion->AddDiffusiveFlux(dir,t, this->FluxRiemann); + this->thermalDiffusion->AddDiffusiveFlux(dir,t, this->FluxRiemann[dir]); } if( (bragViscosityStatus.isExplicit && (!data->rklCycle)) || (bragViscosityStatus.isRKL && data->rklCycle)) { - this->bragViscosity->AddBragViscousFlux(dir,t, this->FluxRiemann); + this->bragViscosity->AddBragViscousFlux(dir,t, this->FluxRiemann[dir]); } // Add braginskii thermal diffusion if( (bragThermalDiffusionStatus.isExplicit && (!data->rklCycle)) || (bragThermalDiffusionStatus.isRKL && data->rklCycle)) { - this->bragThermalDiffusion->AddBragDiffusiveFlux(dir,t, this->FluxRiemann); + this->bragThermalDiffusion->AddBragDiffusiveFlux(dir,t, this->FluxRiemann[dir]); } idfx::popRegion(); diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp index 38ef76d5d..357fcaaca 100644 --- a/src/fluid/calcRightHandSide.hpp +++ b/src/fluid/calcRightHandSide.hpp @@ -21,7 +21,7 @@ struct Fluid_CorrectFluxFunctor { explicit Fluid_CorrectFluxFunctor (Fluid *hydro, real dt) { Uc = hydro->Uc; Vc = hydro->Vc; - Flux = hydro->FluxRiemann; + Flux = hydro->FluxRiemann[dir]; A = hydro->data->A[dir]; dV = hydro->data->dV; x1m = hydro->data->xl[IDIR]; @@ -205,7 +205,7 @@ struct Fluid_CalcRHSFunctor { explicit Fluid_CalcRHSFunctor (Fluid *hydro, real dt) { Uc = hydro->Uc; Vc = hydro->Vc; - Flux = hydro->FluxRiemann; + Flux = hydro->FluxRiemann[dir]; A = hydro->data->A[dir]; dV = hydro->data->dV; x1m = hydro->data->xl[IDIR]; diff --git a/src/fluid/evolveStage.hpp b/src/fluid/evolveStage.hpp index 31dae5b7d..b0a467ba4 100644 --- a/src/fluid/evolveStage.hpp +++ b/src/fluid/evolveStage.hpp @@ -14,20 +14,20 @@ template template void Fluid::LoopDir(const real t, const real dt) { // Step 2: compute the intercell flux with our Riemann solver, store the resulting InvDt - this->rSolver->template CalcFlux(this->FluxRiemann); + this->rSolver->template CalcFlux(this->FluxRiemann[dir]); // Step 2.5: compute intercell parabolic flux when needed if(haveExplicitParabolicTerms) CalcParabolicFlux(t); // If we have tracers, compute the tracer intercell flux if(haveTracer) { - this->tracer->template CalcFlux(this->FluxRiemann); + this->tracer->template CalcFlux(this->FluxRiemann[dir]); } // Step 3: compute the resulting evolution of the conserved variables, stored in Uc CalcRightHandSide(t,dt); if(haveTracer) { - this->tracer->template CalcRightHandSide(this->FluxRiemann,t ,dt); + this->tracer->template CalcRightHandSide(this->FluxRiemann[dir],t ,dt); } // Recursive: do next dimension diff --git a/src/fluid/fluid.hpp b/src/fluid/fluid.hpp index 44cace6bf..e9c9f7d59 100644 --- a/src/fluid/fluid.hpp +++ b/src/fluid/fluid.hpp @@ -67,7 +67,7 @@ class Fluid { void EvolveStage(const real, const real); void ResetStage(); void ShowConfig(); - IdefixArray4D GetFlux() {return this->FluxRiemann;} + IdefixArray4D GetFlux(int dir) {return this->FluxRiemann.at(dir);} int CheckNan(); // Our boundary conditions @@ -182,7 +182,7 @@ class Fluid { // Required by time integrator IdefixArray3D InvDt; - IdefixArray4D FluxRiemann; + std::array,DIMENSIONS> FluxRiemann; IdefixArray3D dMax; // Maximum diffusion speed std::unique_ptr> rSolver; @@ -539,8 +539,12 @@ Fluid::Fluid(Grid &grid, Input &input, DataBlock *datain, int n) { data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); dMax = IdefixArray3D(prefix+"_dMax", data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); - FluxRiemann = IdefixArray4D(prefix+"_FluxRiemann", Phys::nvar+nTracer, - data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); + + for(int i = 0 ; i < DIMENSIONS ; i++) { + FluxRiemann[i] = IdefixArray4D(prefix+"_FluxRiemann_X"+std::to_string(i), + Phys::nvar+nTracer, + data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); + } if constexpr(Phys::mhd) { Vs = IdefixArray4D(prefix+"_Vs", DIMENSIONS, diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp index c8dd8cb96..a2636edaa 100644 --- a/src/rkl/rkl.hpp +++ b/src/rkl/rkl.hpp @@ -30,7 +30,7 @@ class RKLegendre { RKLegendre(Input &, Fluid*); void Cycle(); void ResetStage(); - void ResetFlux(); + void ResetFlux(int dir); void EvolveStage(real); template void CalcParabolicRHS(real); void ComputeDt(); @@ -555,9 +555,9 @@ void RKLegendre::Cycle() { template -void RKLegendre::ResetFlux() { +void RKLegendre::ResetFlux(int dir) { idfx::pushRegion("RKLegendre::ResetFlux"); - IdefixArray4D Flux = hydro->FluxRiemann; + IdefixArray4D Flux = hydro->FluxRiemann[dir]; IdefixArray1D vars = this->varList; idefix_for("RKL_ResetFlux", 0,nvarRKL, @@ -576,7 +576,6 @@ template struct RKLegendre_ResetStageFunctor { explicit RKLegendre_ResetStageFunctor(RKLegendre *rkl) { dU = rkl->dU; - Flux = rkl->hydro->FluxRiemann; vars = rkl->varList; stage = rkl->stage; nvar = rkl->nvarRKL; @@ -596,7 +595,6 @@ struct RKLegendre_ResetStageFunctor { } IdefixArray4D dU; - IdefixArray4D Flux; IdefixArray1D vars; IdefixArray4D dA, dB; IdefixArray3D ex,ey,ez; @@ -680,7 +678,7 @@ void RKLegendre::ComputeDt() { template template void RKLegendre::LoopDir(real t) { - ResetFlux(); + ResetFlux(dir); // CalcParabolicFlux hydro->template CalcParabolicFlux(t); @@ -724,7 +722,7 @@ template void RKLegendre::CalcParabolicRHS(real t) { idfx::pushRegion("RKLegendre::CalcParabolicRHS"); - IdefixArray4D Flux = hydro->FluxRiemann; + IdefixArray4D Flux = hydro->FluxRiemann[dir]; IdefixArray3D A = data->A[dir]; IdefixArray3D dV = data->dV; IdefixArray1D x1m = data->xl[IDIR]; From 04932c13ca290f8386d6356489231273174578fb Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 22 Apr 2025 15:31:08 +0200 Subject: [PATCH 2/4] fix linting issues --- src/fluid/fluid.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fluid/fluid.hpp b/src/fluid/fluid.hpp index e9c9f7d59..9e0615dee 100644 --- a/src/fluid/fluid.hpp +++ b/src/fluid/fluid.hpp @@ -539,7 +539,7 @@ Fluid::Fluid(Grid &grid, Input &input, DataBlock *datain, int n) { data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); dMax = IdefixArray3D(prefix+"_dMax", data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); - + for(int i = 0 ; i < DIMENSIONS ; i++) { FluxRiemann[i] = IdefixArray4D(prefix+"_FluxRiemann_X"+std::to_string(i), Phys::nvar+nTracer, From b32367c1f64857fc1ac07134428b453b2fc740eb Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 22 Apr 2025 16:51:56 +0200 Subject: [PATCH 3/4] fix Uniform collapse test --- test/MHD/AmbipolarWind/setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/MHD/AmbipolarWind/setup.cpp b/test/MHD/AmbipolarWind/setup.cpp index 78620f6e2..bde4ada97 100644 --- a/test/MHD/AmbipolarWind/setup.cpp +++ b/test/MHD/AmbipolarWind/setup.cpp @@ -335,7 +335,7 @@ void EmfBoundary(DataBlock& data, const real t) { } void FluxBoundary(DataBlock & data, int dir, BoundarySide side, const real t) { - IdefixArray4D Flux = data.hydro->FluxRiemann; + IdefixArray4D Flux = data.hydro->FluxRiemann[dir]; if( dir==IDIR && side == left) { int iref = data.beg[IDIR]; From 813108bd7c334e96b8bac676bbe1e22bd4cbf08a Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 23 Apr 2025 09:25:42 +0200 Subject: [PATCH 4/4] fix uniform collapse test (forgot to save...) --- test/SelfGravity/UniformCollapse/setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/SelfGravity/UniformCollapse/setup.cpp b/test/SelfGravity/UniformCollapse/setup.cpp index 6f773a153..22f4a34d6 100644 --- a/test/SelfGravity/UniformCollapse/setup.cpp +++ b/test/SelfGravity/UniformCollapse/setup.cpp @@ -65,7 +65,7 @@ void FluxBoundary(Fluid *hydro, int dir, BoundarySide side, cons if((dir==IDIR) && (side == left)) { // Loading needed data DataBlock &data = *hydro->data; - IdefixArray4D Flux = hydro->FluxRiemann; + IdefixArray4D Flux = hydro->FluxRiemann[dir]; real halfDt = data.dt/2.; // RK2, dt is actually half at each flux calculation int iref = data.nghost[IDIR]; real rin = data.xbeg[IDIR];