Add API to normalise metrics

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

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);
 }

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