diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index a090dcd959..cfcd9368d4 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -564,17 +564,26 @@ std::string parallelSliceFieldName(std::string_view field, int offset); /// tensor, which is multiplied by the normalisation factor. struct MetricNormaliser { std::optional g = std::nullopt; + bool g_mul = false; std::optional g11 = std::nullopt; + bool g11_mul = false; std::optional g22 = std::nullopt; + bool g22_mul = false; std::optional g33 = std::nullopt; + bool g33_mul = false; std::optional g12 = std::nullopt; + bool g12_mul = false; std::optional g13 = std::nullopt; + bool g13_mul = false; std::optional g23 = std::nullopt; + bool g23_mul = false; std::optional dx = std::nullopt; std::optional dy = std::nullopt; std::optional dz = std::nullopt; std::optional J = std::nullopt; + bool J_mul = false; std::optional Bxy = std::nullopt; + bool Bxy_mul = false; }; #endif // BOUT_COORDINATES_H diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index bad8ebbd84..66f4195ca1 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1313,16 +1313,32 @@ void Coordinates::normaliseMetric(const MetricNormaliser& norm) { if (norm.J.has_value()) { if (J().hasParallelSlices()) { - setJ(FieldMetricParallel{J() / *norm.J}); + if (norm.J_mul) { + setJ(FieldMetricParallel{J() * *norm.J}); + } else { + setJ(FieldMetricParallel{J() / *norm.J}); + } } else { - setJ(J() / *norm.J); + if (norm.J_mul) { + setJ(J() * *norm.J); + } else { + setJ(J() / *norm.J); + } } } if (norm.Bxy.has_value()) { if (Bxy().hasParallelSlices()) { - setBxy(FieldMetricParallel{Bxy() / *norm.Bxy}); + if (norm.Bxy_mul) { + setBxy(FieldMetricParallel{Bxy() * *norm.Bxy}); + } else { + setBxy(FieldMetricParallel{Bxy() / *norm.Bxy}); + } } else { - setBxy(Bxy() / *norm.Bxy); + if (norm.Bxy_mul) { + setBxy(Bxy() * *norm.Bxy); + } else { + setBxy(Bxy() / *norm.Bxy); + } } } if (norm.dx.has_value()) { @@ -1334,7 +1350,7 @@ void Coordinates::normaliseMetric(const MetricNormaliser& norm) { if (norm.dz.has_value()) { setDz(dz() / *norm.dz); } - invalidateMetricCaches(); + recalculateAndReset(false, false); if (norm.g.has_value() or norm.g22.has_value()) { if (Bxy().isFci()) { // No we compute g_22_* - they must not be cleared. If they get @@ -1343,9 +1359,14 @@ void Coordinates::normaliseMetric(const MetricNormaliser& norm) { g_22_ylow(); g_22_yhigh(); ASSERT2(_g_22_ylow.has_value()); - (*_g_22_ylow) /= g22; ASSERT2(_g_22_yhigh.has_value()); - (*_g_22_yhigh) /= g22; + if (norm.g.has_value() ? norm.g_mul : norm.g22_mul) { + (*_g_22_ylow) *= g22; + (*_g_22_yhigh) *= g22; + } else { + (*_g_22_ylow) /= g22; + (*_g_22_yhigh) /= g22; + } } } } diff --git a/src/mesh/metric_tensor.cxx b/src/mesh/metric_tensor.cxx index 1347f88a7e..29b511e2cf 100644 --- a/src/mesh/metric_tensor.cxx +++ b/src/mesh/metric_tensor.cxx @@ -152,40 +152,56 @@ auto ContravariantMetricTensor::inverse(const std::string& region, bool communic template void MetricTensor::normaliseMetric(const MetricNormaliser& norm, const F& op) { if (norm.g.has_value()) { - op(g11_m, norm.g); - op(g22_m, norm.g); - op(g33_m, norm.g); - op(g12_m, norm.g); - op(g13_m, norm.g); - op(g23_m, norm.g); + op(g11_m, norm.g, norm.g_mul); + op(g22_m, norm.g, norm.g_mul); + op(g33_m, norm.g, norm.g_mul); + op(g12_m, norm.g, norm.g_mul); + op(g13_m, norm.g, norm.g_mul); + op(g23_m, norm.g, norm.g_mul); } else { - op(g11_m, norm.g11); - op(g22_m, norm.g22); - op(g33_m, norm.g33); - op(g12_m, norm.g12); - op(g13_m, norm.g13); - op(g23_m, norm.g23); + op(g11_m, norm.g11, norm.g11_mul); + op(g22_m, norm.g22, norm.g22_mul); + op(g33_m, norm.g33, norm.g33_mul); + op(g12_m, norm.g12, norm.g12_mul); + op(g13_m, norm.g13, norm.g13_mul); + op(g23_m, norm.g23, norm.g23_mul); } } void ContravariantMetricTensor::normaliseMetric(const MetricNormaliser& norm) { - MetricTensor::normaliseMetric(norm, [](FieldMetric& f, auto fac) { + MetricTensor::normaliseMetric(norm, [](FieldMetric& f, auto fac, bool fac_mul) { if (fac.has_value()) { if (f.hasParallelSlices()) { - f.asField3DParallel() *= fac.value(); + if (fac_mul) { + f.asField3DParallel() /= fac.value(); + } else { + f.asField3DParallel() *= fac.value(); + } } else { - f *= fac.value(); + if (fac_mul) { + f /= fac.value(); + } else { + f *= fac.value(); + } } } }); } void CovariantMetricTensor::normaliseMetric(const MetricNormaliser& norm) { - MetricTensor::normaliseMetric(norm, [](FieldMetric& f, auto fac) { + MetricTensor::normaliseMetric(norm, [](FieldMetric& f, auto fac, bool fac_mul) { if (fac.has_value()) { if (f.hasParallelSlices()) { - f.asField3DParallel() /= fac.value(); + if (fac_mul) { + f.asField3DParallel() *= fac.value(); + } else { + f.asField3DParallel() /= fac.value(); + } } else { - f /= fac.value(); + if (fac_mul) { + f *= fac.value(); + } else { + f /= fac.value(); + } } } }); diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index 8f651eec23..dd5e2f0bf2 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -81,14 +81,18 @@ MetricNormaliser TokamakOrFCIMetricNormaliser(const Mesh* mesh, BoutReal Bnorm, if (mesh->isFci()) { return {.g{SQ(rho_s0)}, .J{rho_s0 * rho_s0 * rho_s0}, .Bxy{Bnorm}}; } - return {.g11{1 / SQ(Bnorm * rho_s0)}, + return {.g11{SQ(Bnorm * rho_s0)}, + .g11_mul{true}, .g22{SQ(rho_s0)}, .g33{SQ(rho_s0)}, - .g12{1 / Bnorm}, - .g13{1 / Bnorm}, + .g12{Bnorm}, + .g12_mul{true}, + .g13{Bnorm}, + .g13_mul{true}, .g23{SQ(rho_s0)}, .dx{rho_s0 * rho_s0 * Bnorm}, - .J{rho_s0 / Bnorm}, + .J{Bnorm / rho_s0}, + .J_mul{true}, .Bxy{Bnorm}}; }