diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index f01e3883e0..dc88d20b3c 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -372,6 +372,13 @@ public: void communicateMetricTensor(); +private: + void normaliseMetricTokamak(BoutReal rho_s0, BoutReal Bnorm); + void normaliseMetricFCI(BoutReal rho_s0, BoutReal Bnorm); + +public: + void normaliseMetric(BoutReal rho_s0, BoutReal Bnorm); + void communicateDz(); ///< Coordinate system Jacobian, so volume of cell is J*dx*dy*dz diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 5e461602d0..bacbbe81c4 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -823,6 +823,7 @@ Coordinates::FieldMetric Coordinates::recalculateJacobian() const { } Coordinates::FieldMetric Coordinates::recalculateBxy() const { + ASSERT2(not J().isFci()); return sqrt(g_22()) / J(); } @@ -1277,6 +1278,71 @@ void Coordinates::communicateMetricTensor() { void Coordinates::communicateDz() { localmesh->communicate(dz_); } void Coordinates::splitBxyParallelSlices() { - Bxy_.splitParallelSlices(); - Bxy_.yup() = Bxy_.ydown() = Bxy_; + if (not Bxy_.hasParallelSlices()) { + auto copy = Bxy_; + Bxy_.splitParallelSlices(); + Bxy_.yup() = Bxy_.ydown() = copy; + } +} + +void Coordinates::normaliseMetricTokamak(const BoutReal rho_s0, const BoutReal Bnorm) { + contravariantMetricTensor = ContravariantMetricTensor( + contravariantMetricTensor.g11() / SQ(Bnorm * rho_s0), + contravariantMetricTensor.g22() * SQ(rho_s0), + contravariantMetricTensor.g33() * SQ(rho_s0), + contravariantMetricTensor.g12() / Bnorm, contravariantMetricTensor.g13() / Bnorm, + contravariantMetricTensor.g23() * SQ(rho_s0)); + + covariantMetricTensor = CovariantMetricTensor( + covariantMetricTensor.g11() * SQ(Bnorm * rho_s0), + covariantMetricTensor.g22() / SQ(rho_s0), covariantMetricTensor.g33() / SQ(rho_s0), + covariantMetricTensor.g12() * Bnorm, covariantMetricTensor.g13() * Bnorm, + covariantMetricTensor.g23() / SQ(rho_s0)); + + setDx(dx() / (rho_s0 * rho_s0 * Bnorm)); + setBxy(Bxy() / Bnorm); + setJ(J() * Bnorm / rho_s0); + invalidateMetricCaches(); +} + +void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm) { + BoutReal rhoSQ = SQ(rho_s0); + + // g??.asField3DParallel() *= rhoSQ; + contravariantMetricTensor.map([&](auto f) -> FieldMetric { + if (f.hasParallelSlices()) { + return f.asField3DParallel() * rhoSQ; + }; + return f * rhoSQ; + }); + + // g_??.asField3DParallel() /= rhoSQ; + covariantMetricTensor.map([&](auto f) -> FieldMetric { + if (f.hasParallelSlices()) { + return f.asField3DParallel() / rhoSQ; + }; + return f / rhoSQ; + }); + // J.asField3DParallel() /= rho_s0 * rho_s0 * rho_s0; +#if BOUT_USE_METRIC_3D + if (J().hasParallelSlices()) { + setJ(Field3DParallel{J() / (rho_s0 * rho_s0 * rho_s0)}); + } else { + setJ(J() / (rho_s0 * rho_s0 * rho_s0)); + } + setBxy(Field3DParallel(Bxy() / Bnorm)); +#else + setJ(J() / (rho_s0 * rho_s0 * rho_s0)); + setBxy(Bxy() / Bnorm); +#endif + invalidateMetricCaches(); +} + +void Coordinates::normaliseMetric(const BoutReal rho_s0, const BoutReal Bnorm) { + ASSERT2(hasParallelTransform()); + if (Bxy().isFci()) { + normaliseMetricFCI(rho_s0, Bnorm); + return; + } + normaliseMetricTokamak(rho_s0, Bnorm); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index e6913a67eb..00e6aecf68 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -762,10 +762,13 @@ Field3D Laplace(const Field3D& f, CELL_LOC outloc, * Inverse of Laplacian operator in LaplaceXY solver *******************************************************************************/ -Field2D Laplace_perpXY(const Field2D& A, const Field2D& f) { #if BOUT_USE_METRIC_3D +Field2D Laplace_perpXY([[maybe_unused]] const Field2D& A, + [[maybe_unused]] const Field2D& f) { throw BoutException("Coordinates::Laplace_perpXY for 3D metric not implemented"); +} #else +Field2D Laplace_perpXY(const Field2D& A, const Field2D& f) { const auto& coords = *f.getCoordinates(); Field2D result; @@ -821,8 +824,8 @@ Field2D Laplace_perpXY(const Field2D& A, const Field2D& f) { } return result; -#endif } +#endif /******************************************************************************* * b0xGrad_dot_Grad