Fargo#127
Conversation
… remaps in the X2 and X3 directions rather than try to reuse one routine but pass k-slices or j-slices.
|
@pdmullen tests now pass (I fixed one issue), so I think this is ready for real review. |
| using ArtemisUtils::PLM; | ||
| using ArtemisUtils::VI; | ||
|
|
||
| #define CUB(x) ((x) * (x) * (x)) |
|
|
||
| //---------------------------------------------------------------------------------------- | ||
| //! template instantiations | ||
| template TaskListStatus Advect<Coordinates::cartesian>(Mesh *, const SimTime &); |
There was a problem hiding this comment.
Totally unrelated... but in the most recent parthenon swarm pack PR, I enrolled macros for template instantiations. We have these in almost every cpp file for different geometries... maybe we could remove a lot of lines of code by similarly adopting such a macro.
| //! \fn Real RotatingFrame::OmegaKep | ||
| //! \brief Returns Keplerian angular velocity at spherical radius R | ||
| KOKKOS_FORCEINLINE_FUNCTION Real OmegaKep(const Real gm, const Real R) { | ||
| return std::sqrt(gm / (R * R * R)); |
| ((GEOM == Coordinates::cartesian) || (GEOM == Coordinates::cylindrical)) ? 1 : 2; | ||
|
|
||
| KOKKOS_FORCEINLINE_FUNCTION Real OmegaKep(const Real R) const { | ||
| return std::sqrt(gm / (R * R * R)); |
There was a problem hiding this comment.
CUBE? Or just get get rid of CUBE?
| if (j >= js && j <= je) v0(b, n, k, j, i) += dq / r.vol; | ||
| if (jp >= js && jp <= je) v0(b, n, k, jp, i) -= dq / rp.vol; |
There was a problem hiding this comment.
Can you help me understand the conditionals here?
| if (k >= ks && k <= ke) v0(b, n, k, j, i) += dq / r.vol; | ||
| if (kp >= ks && kp <= ke) v0(b, n, kp, j, i) -= dq / rp.vol; |
| oa.ApplyGradSkew(rd); | ||
| } | ||
| RemapUpdate(v0, ru, rc, vb, dwdt, three_d, b, n, k, j, j + joff, i); | ||
| RemapUpdateX2<GEOM>(oa, v0, ru, rc, vb, three_d, b, n, k, j, j + joff, i, jb.s, |
There was a problem hiding this comment.
Also need some help here. Are we still doing a sweep as we had done it before. It looks like it, but more indexing enters the RemapUpdate function call.
| // v_stored = v_full - v_bg(R) | ||
| // where v_bg = R*(OmegaKep(R) - omega_f) in the phi direction. The Riemann solver | ||
| // transports the stored momentum, so we must correct for the missing background | ||
| // angular momentum flux: d(rho*v_stored)/dt -= (1/V) * div(F_rho * v_bg). |
There was a problem hiding this comment.
I worry about consistency. Is the Riemann solve mass flux the right quantity to be used for the correction? We are inside an operator split update, it feels a bit odd to be using the mass flux computed from the final stage of e.g., an RK2 integration.
2 participants
Background
This extends our linear advection machinery to work in spherical and cylindrical coordinates where the background velocity is assumed to be keplerian.
To keep the overlap integrals analytic, we make the approximation that in cylindrical/axisymmetric coordinates the radius in the keplerian Omega is the cylindrical R but in spherical is the spherical r.
This was done with the help of GPT (I still need to put the AI statement in the files).
This more or less works, but it's not quite ready yet.
Description of Changes
Checklist
// This file was created in part or in whole by generative AI