Skip to content

fix RKL in grids with small spacing #323

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
glesur merged 6 commits into from
Feb 21, 2025
Merged

fix RKL in grids with small spacing #323

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
34d9366
fix a bug that could lead to Nans when using the RKL scheme on grids
glesur Feb 19, 2025
cb7b999
change the way polar singularities are handled
glesur Feb 19, 2025
2874fda
-refactor curvature terms without using a small number on all Areas
glesur Feb 19, 2025
461ae58
fix typo
glesur Feb 20, 2025
a4ed690
fix capture on Cuda
glesur Feb 20, 2025
b3481f5
do not make an exception when Area=0, as this artificiallykills off one
glesur Feb 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 36 additions & 34 deletions src/fluid/calcRightHandSide.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,47 +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<SMALL_NUMBER) Ax=SMALL_NUMBER; //avoid singularity around poles
// No area for this one
Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax;
#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<SMALL_NUMBER) Ax=SMALL_NUMBER; // 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(Ax<SMALL_NUMBER) Ax=SMALL_NUMBER; // avoid singularity around poles
Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax;
#if GEOMETRY == SPHERICAL
if constexpr(dir==IDIR) {
#if COMPONENTS == 3
Ax[iMPHI] *= FABS(x1m(i));
#endif // COMPONENTS == 3
if constexpr(Phys::mhd) {
EXPAND( ,
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
}
};

Expand Down
89 changes: 62 additions & 27 deletions src/rkl/rkl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -762,32 +762,42 @@ void RKLegendre<Phys>::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);

#if GEOMETRY != CARTESIAN
if(Ax<SMALL_NUMBER)
Ax=SMALL_NUMBER; // Essentially to avoid singularity around poles
#endif

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 && nv==iMPHI) {
// Conserve angular momentum, hence flux is R*Vphi
Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i));
}
#endif // GEOMETRY==POLAR OR CYLINDRICAL
#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) {
Ax *= FABS(x1m(i));
} else if(dir==JDIR && nv==VX3) {
Ax *= FABS(sm(j));
}

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

#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));
}
#endif // GEOMETRY == SPHERICAL && COMPONENTS == 3
Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax;
}
);

Expand Down Expand Up @@ -818,17 +828,42 @@ void RKLegendre<Phys>::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
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_ * (Flux(iBPHI, k, j, i+1) - Flux(iBPHI, k, j, i) );

#elif (GEOMETRY == SPHERICAL)
real q = 1 / (x1_*dx_);
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));
}
real dx_ = dx(j);
real rt_ = rt(i);
if constexpr(Phys::mhd) {
if(nv == iBPHI) {
rhs = - 1 / (rt_*dx_) * (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
Expand Down