From 34d9366527d8870651f268735a2874e9bc38c875 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 19 Feb 2025 21:37:54 +0100 Subject: [PATCH 1/6] fix a bug that could lead to Nans when using the RKL scheme on grids with small grid spacing (<1e-5) --- src/rkl/rkl.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp index a2f6734f..6f679620 100644 --- a/src/rkl/rkl.hpp +++ b/src/rkl/rkl.hpp @@ -763,11 +763,6 @@ void RKLegendre::CalcParabolicRHS(real t) { KOKKOS_LAMBDA (int n, int k, int j, int i) { real Ax = A(k,j,i); -#if GEOMETRY != CARTESIAN - if(Ax Date: Wed, 19 Feb 2025 22:22:19 +0100 Subject: [PATCH 2/6] change the way polar singularities are handled --- src/fluid/calcRightHandSide.hpp | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp index c8738b4b..b2e964ee 100644 --- a/src/fluid/calcRightHandSide.hpp +++ b/src/fluid/calcRightHandSide.hpp @@ -162,9 +162,7 @@ struct Fluid_CorrectFluxFunctor { // Conserve angular momentum, hence flux is R*Bphi Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); if constexpr(Phys::mhd) { - if(Ax0) Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax; //avoid singularity around poles } } #endif // GEOMETRY==POLAR OR CYLINDRICAL @@ -175,17 +173,17 @@ struct Fluid_CorrectFluxFunctor { Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); #endif // COMPONENTS == 3 if constexpr(Phys::mhd) { - if(Ax0) { // avoid singularity around poles + EXPAND( , + Flux(iBTH,k,j,i) = Flux(iBTH,k,j,i) * x1m(i) / Ax; , + Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) * x1m(i) / Ax; ) + } } } else if constexpr (dir==JDIR) { #if COMPONENTS == 3 Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(sinx2m(j)); if constexpr(Phys::mhd) { - if(Ax0) Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax; // avoid singularity around poles } #endif // COMPONENTS = 3 } From 2874fda4bd24157089d1b6940562a8c4a96d35ff Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 19 Feb 2025 23:18:46 +0100 Subject: [PATCH 3/6] -refactor curvature terms without using a small number on all Areas -properly handles curvature on Bphi in 2.5D spherical with RKL --- src/rkl/rkl.hpp | 64 ++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 50 insertions(+), 14 deletions(-) diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp index 6f679620..0bc220c9 100644 --- a/src/rkl/rkl.hpp +++ b/src/rkl/rkl.hpp @@ -770,19 +770,35 @@ void RKLegendre::CalcParabolicRHS(real t) { // Curvature terms #if (GEOMETRY == POLAR && COMPONENTS >= 2) \ || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) - if(dir==IDIR && nv==iMPHI) { - // Conserve angular momentum, hence flux is R*Vphi - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); + if(dir==IDIR) { + if(nv==iMPHI) { + Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1(i)); + } + if constexpr(Phys::mhd) { + if(nv==iBPHI) { + // No area for this one + if(Ax > 0) Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax;//avoid singularity around poles + } + } } #endif // GEOMETRY==POLAR OR CYLINDRICAL -#if GEOMETRY == SPHERICAL && COMPONENTS == 3 - if(dir==IDIR && nv==iMPHI) { - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - } else if(dir==JDIR && nv==iMPHI) { - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(sm(j)); +#if GEOMETRY == SPHERICAL + if(dir==IDIR && nv==VX3) { + Flux(nv,k,j,i) = Flux(nv,k,j,i) * FABS(x1m(i)); + } else if(dir==JDIR && nv==VX3) { + Flux(nv,k,j,i) = Flux(nv,k,j,i) * FABS(sm(j)); + } + + if constexpr(Phys::mhd) { + if(dir == IDIR && Ax > 0 && (nv==BX3 || nv == BX2)) { + Flux(nv,k,j,i) = Flux(nv,k,j,i) * x1m(i) / Ax; + } + if(dir==JDIR && Ax > 0 && nv==BX3) { + Flux(nv,k,j,i) = Flux(nv,k,j,i) / Ax; + } } -#endif // GEOMETRY == SPHERICAL && COMPONENTS == 3 +#endif // GEOMETRY == SPHERICAL$ } ); @@ -813,17 +829,37 @@ void RKLegendre::CalcParabolicRHS(real t) { } #if GEOMETRY != CARTESIAN - #ifdef iMPHI - if((dir==IDIR) && (nv == iMPHI)) { + if(dir==IDIR) { + #ifdef iMPHI + if((nv == iMPHI)) { rhs /= x1(i); } + #endif // iMPHI + if constexpr(Phys::mhd) { + #if (GEOMETRY == POLAR || GEOMETRY == CYLINDRICAL) && (defined iBPHI) + if(nv==iBPHI) rhs = - 1 / dx(i) * (Flux(iBPHI, k, j, i+1) - Flux(iBPHI, k, j, i) ); + + #elif (GEOMETRY == SPHERICAL) + real q = 1 / (x1(i)*dx(i)); + if(nv == BX2 || nv == BX3) { + rhs = -q * ((Flux(nv, k, j, i+1) - Flux(nv, k, j, i) )); + } + #endif // GEOMETRY + } // MHD + } // dir==IDIR + if(dir==JDIR) { #if (GEOMETRY == SPHERICAL) && (COMPONENTS == 3) - if((dir==JDIR) && (nv == iMPHI)) { + if(nv == iMPHI) { rhs /= FABS(s(j)); } + if constexpr(Phys::mhd) { + if(nv == iBPHI) { + rhs = - 1 / (rt(i)*dx(j)) * (Flux(nv, k, j+1, i) - Flux(nv, k, j, i)); + } + } #endif // GEOMETRY - // Nothing for KDIR - #endif // iMPHI + } // dir==JDIR + // Nothing for KDIR #endif // GEOMETRY != CARTESIAN // store the field components From 461ae58a13c2a6734ae0e0faece47e828cc2d32b Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 20 Feb 2025 09:14:14 +0100 Subject: [PATCH 4/6] fix typo --- src/rkl/rkl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp index 0bc220c9..1ebebe5a 100644 --- a/src/rkl/rkl.hpp +++ b/src/rkl/rkl.hpp @@ -772,7 +772,7 @@ void RKLegendre::CalcParabolicRHS(real t) { || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) if(dir==IDIR) { if(nv==iMPHI) { - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1(i)); + Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); } if constexpr(Phys::mhd) { if(nv==iBPHI) { From a4ed690b6acf4f4e23351a0cda7dcf2270291b15 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 20 Feb 2025 09:45:34 +0100 Subject: [PATCH 5/6] fix capture on Cuda --- src/rkl/rkl.hpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp index 1ebebe5a..3a7d3767 100644 --- a/src/rkl/rkl.hpp +++ b/src/rkl/rkl.hpp @@ -835,12 +835,15 @@ void RKLegendre::CalcParabolicRHS(real t) { rhs /= x1(i); } #endif // iMPHI + real dx_ = dx(i); + real x1_ = x1(i); + if constexpr(Phys::mhd) { #if (GEOMETRY == POLAR || GEOMETRY == CYLINDRICAL) && (defined iBPHI) - if(nv==iBPHI) rhs = - 1 / dx(i) * (Flux(iBPHI, k, j, i+1) - Flux(iBPHI, k, j, i) ); + if(nv==iBPHI) rhs = - 1 / dx_ * (Flux(iBPHI, k, j, i+1) - Flux(iBPHI, k, j, i) ); #elif (GEOMETRY == SPHERICAL) - real q = 1 / (x1(i)*dx(i)); + real q = 1 / (x1_*dx_); if(nv == BX2 || nv == BX3) { rhs = -q * ((Flux(nv, k, j, i+1) - Flux(nv, k, j, i) )); } @@ -852,9 +855,11 @@ void RKLegendre::CalcParabolicRHS(real t) { if(nv == iMPHI) { rhs /= FABS(s(j)); } + real dx_ = dx(j); + real rt_ = rt(i); if constexpr(Phys::mhd) { if(nv == iBPHI) { - rhs = - 1 / (rt(i)*dx(j)) * (Flux(nv, k, j+1, i) - Flux(nv, k, j, i)); + rhs = - 1 / (rt_*dx_) * (Flux(nv, k, j+1, i) - Flux(nv, k, j, i)); } } #endif // GEOMETRY From b3481f52365ce515c733f7b6e66a36c1214f451f Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 20 Feb 2025 13:02:33 +0100 Subject: [PATCH 6/6] do not make an exception when Area=0, as this artificiallykills off one of the EMF component when axisymmetric --- src/fluid/calcRightHandSide.hpp | 66 +++++++++++++++++---------------- src/rkl/rkl.hpp | 57 ++++++++++++++-------------- 2 files changed, 63 insertions(+), 60 deletions(-) diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp index b2e964ee..38ef76d5 100644 --- a/src/fluid/calcRightHandSide.hpp +++ b/src/fluid/calcRightHandSide.hpp @@ -149,45 +149,49 @@ struct Fluid_CorrectFluxFunctor { Flux(MX1+meanDir,k,j,i) += meanV * Flux(RHO,k,j,i); } // Fargo & Rotation corrections - real Ax = A(k,j,i); + ////////////////////////////////////////////// + // Define correction factor for the fluxes + ////////////////////////////////////////////// - for(int nv = 0 ; nv < Phys::nvar ; nv++) { - Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax; - } + real Ax[Phys::nvar]; // corrected Area + for(int nv = 0 ; nv < Phys::nvar ; nv++) Ax[nv] = A(k,j,i); // Curvature terms -#if (GEOMETRY == POLAR && COMPONENTS >= 2) \ - || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) - if constexpr (dir==IDIR) { - // Conserve angular momentum, hence flux is R*Bphi - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - if constexpr(Phys::mhd) { - if(Ax>0) Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax; //avoid singularity around poles + #if (GEOMETRY == POLAR && COMPONENTS >= 2) \ + || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) + if constexpr (dir==IDIR) { + // Conserve angular momentum, hence flux is R*Bphi + Ax[iMPHI] *= FABS(x1m(i)); + if constexpr(Phys::mhd) { + Ax[iBPHI] = 1.0; // corrected Flux is simply Bphi + } } - } -#endif // GEOMETRY==POLAR OR CYLINDRICAL + #endif // GEOMETRY==POLAR OR CYLINDRICAL -#if GEOMETRY == SPHERICAL - if constexpr(dir==IDIR) { - #if COMPONENTS == 3 - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - #endif // COMPONENTS == 3 - if constexpr(Phys::mhd) { - if(Ax>0) { // avoid singularity around poles + #if GEOMETRY == SPHERICAL + if constexpr(dir==IDIR) { + #if COMPONENTS == 3 + Ax[iMPHI] *= FABS(x1m(i)); + #endif // COMPONENTS == 3 + if constexpr(Phys::mhd) { EXPAND( , - Flux(iBTH,k,j,i) = Flux(iBTH,k,j,i) * x1m(i) / Ax; , - Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) * x1m(i) / Ax; ) - } - } - } else if constexpr (dir==JDIR) { - #if COMPONENTS == 3 - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(sinx2m(j)); - if constexpr(Phys::mhd) { - if(Ax>0) Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax; // avoid singularity around poles + Ax[iBTH] = x1m(i); , + Ax[iBPHI] = x1m(i); ) + }// Phys::mhd + } else if constexpr (dir==JDIR) { + #if COMPONENTS == 3 + Ax[iMPHI] *= FABS(sinx2m(j)); + if constexpr(Phys::mhd) { + Ax[iBPHI] = 1.0; // corrected Flux is simply Bphi + } + #endif // COMPONENTS = 3 } - #endif // COMPONENTS = 3 + #endif // GEOMETRY == SPHERICAL + + // Finally correct the flux + for(int nv = 0 ; nv < Phys::nvar ; nv++) { + Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax[nv]; } -#endif // GEOMETRY == SPHERICAL } }; diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp index 3a7d3767..c8dd8cb9 100644 --- a/src/rkl/rkl.hpp +++ b/src/rkl/rkl.hpp @@ -762,43 +762,42 @@ void RKLegendre::CalcParabolicRHS(real t) { data->beg[IDIR],data->end[IDIR]+ioffset, KOKKOS_LAMBDA (int n, int k, int j, int i) { real Ax = A(k,j,i); - const int nv = varList(n); - Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax; - // Curvature terms -#if (GEOMETRY == POLAR && COMPONENTS >= 2) \ - || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) - if(dir==IDIR) { - if(nv==iMPHI) { - Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i)); - } - if constexpr(Phys::mhd) { - if(nv==iBPHI) { - // No area for this one - if(Ax > 0) Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax;//avoid singularity around poles + #if (GEOMETRY == POLAR && COMPONENTS >= 2) \ + || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3) + if(dir==IDIR) { + if(nv==iMPHI) { + Ax *= FABS(x1m(i)); + } + if constexpr(Phys::mhd) { + if(nv==iBPHI) { + // No area for this one + Ax = 1; + } } } - } -#endif // GEOMETRY==POLAR OR CYLINDRICAL - -#if GEOMETRY == SPHERICAL - if(dir==IDIR && nv==VX3) { - Flux(nv,k,j,i) = Flux(nv,k,j,i) * FABS(x1m(i)); - } else if(dir==JDIR && nv==VX3) { - Flux(nv,k,j,i) = Flux(nv,k,j,i) * FABS(sm(j)); - } + #endif // GEOMETRY==POLAR OR CYLINDRICAL - if constexpr(Phys::mhd) { - if(dir == IDIR && Ax > 0 && (nv==BX3 || nv == BX2)) { - Flux(nv,k,j,i) = Flux(nv,k,j,i) * x1m(i) / Ax; + #if GEOMETRY == SPHERICAL + if(dir==IDIR && nv==VX3) { + Ax *= FABS(x1m(i)); + } else if(dir==JDIR && nv==VX3) { + Ax *= FABS(sm(j)); } - if(dir==JDIR && Ax > 0 && nv==BX3) { - Flux(nv,k,j,i) = Flux(nv,k,j,i) / Ax; + + if constexpr(Phys::mhd) { + if(dir == IDIR && (nv==BX3 || nv == BX2)) { + Ax = x1m(i); + } + if(dir==JDIR && nv==BX3) { + Ax = 1.0; + } } - } -#endif // GEOMETRY == SPHERICAL$ + #endif // GEOMETRY == SPHERICAL + + Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax; } );