Refactor coordinates

CMakeLists.txt

@@ -352,9 +352,15 @@ set(BOUT_SOURCES
   ./src/sys/timer.cxx
   ./src/sys/type_name.cxx
   ./src/sys/utils.cxx
+  ./include/bout/metric_tensor.hxx
+  ./src/mesh/metric_tensor.cxx
+  ./include/bout/christoffel_symbols.hxx
+  ./src/mesh/christoffel_symbols.cxx
+  ./include/bout/g_values.hxx
+  ./src/mesh/g_values.cxx
   ${CMAKE_CURRENT_BINARY_DIR}/include/bout/revision.hxx
   ${CMAKE_CURRENT_BINARY_DIR}/include/bout/version.hxx
-  )
+)
 
 
 find_package(Python3)

examples/6field-simple/elm_6f.cxx

@@ -361,7 +361,7 @@ class Elm_6f : public PhysicsModel {
       result.allocate();
       for (auto i : result) {
         result[i] =
-            (fp[i.yp()] - fm[i.ym()]) / (2. * coord->dy[i] * sqrt(coord->g_22[i]));
+            (fp[i.yp()] - fm[i.ym()]) / (2. * coord->dy()[i] * sqrt(coord->g_22()[i]));
       }
     } else {
       result = Grad_par(f, loc);
@@ -705,7 +705,7 @@ class Elm_6f : public PhysicsModel {
 
     if (mesh->IncIntShear) {
       // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
-      coord->IntShiftTorsion = I;
+      coord->setIntShiftTorsion(I);
 
     } else {
       // Dimits style, using local coordinate system
@@ -867,7 +867,7 @@ class Elm_6f : public PhysicsModel {
     Btxy /= Bbar;
     B0 /= Bbar;
     hthe /= Lbar;
-    coord->dx /= Lbar * Lbar * Bbar;
+    coord->setDx(coord->dx() / (Lbar * Lbar * Bbar));
     I *= Lbar * Lbar * Bbar;
 
     if ((!T0_fake_prof) && n0_fake_prof) {
@@ -1051,24 +1051,25 @@ class Elm_6f : public PhysicsModel {
 
     /**************** CALCULATE METRICS ******************/
 
-    coord->g11 = SQ(Rxy * Bpxy);
-    coord->g22 = 1.0 / SQ(hthe);
-    coord->g33 = SQ(I) * coord->g11 + SQ(B0) / coord->g11;
-    coord->g12 = 0.0;
-    coord->g13 = -I * coord->g11;
-    coord->g23 = -Btxy / (hthe * Bpxy * Rxy);
+    const auto g11 = SQ(Rxy * Bpxy);
+    const auto g22 = 1.0 / SQ(hthe);
+    const auto g33 = SQ(I) * g11 + SQ(B0) / g11;
+    const auto g12 = 0.0;
+    const auto g13 = -I * g11;
+    const auto g23 = -Btxy / (hthe * Bpxy * Rxy);
 
-    coord->J = hthe / Bpxy;
-    coord->Bxy = B0;
+    const auto g_11 = 1.0 / g11 + SQ(I * Rxy);
+    const auto g_22 = SQ(B0 * hthe / Bpxy);
+    const auto g_33 = Rxy * Rxy;
+    const auto g_12 = Btxy * hthe * I * Rxy / Bpxy;
+    const auto g_13 = I * Rxy * Rxy;
+    const auto g_23 = Btxy * hthe * Rxy / Bpxy;
 
-    coord->g_11 = 1.0 / coord->g11 + SQ(I * Rxy);
-    coord->g_22 = SQ(B0 * hthe / Bpxy);
-    coord->g_33 = Rxy * Rxy;
-    coord->g_12 = Btxy * hthe * I * Rxy / Bpxy;
-    coord->g_13 = I * Rxy * Rxy;
-    coord->g_23 = Btxy * hthe * Rxy / Bpxy;
+    coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
+                           CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));
 
-    coord->geometry(); // Calculate quantities from metric tensor
+    coord->setJ(hthe / Bpxy);
+    coord->setBxy(B0);
 
     // Set B field vector
 
@@ -1446,11 +1447,11 @@ class Elm_6f : public PhysicsModel {
       if (hyperviscos > 0.0) {
         // Calculate coefficient.
 
-        hyper_mu_x = hyperviscos * coord->g_11 * SQ(coord->dx)
-                     * abs(coord->g11 * D2DX2(U)) / (abs(U) + 1e-3);
+        hyper_mu_x = hyperviscos * coord->g_11() * SQ(coord->dx())
+                     * abs(coord->g11() * D2DX2(U)) / (abs(U) + 1e-3);
         hyper_mu_x.applyBoundary("dirichlet"); // Set to zero on all boundaries
 
-        ddt(U) += hyper_mu_x * coord->g11 * D2DX2(U);
+        ddt(U) += hyper_mu_x * coord->g11() * D2DX2(U);
 
         if (first_run) {
           // Print out maximum values of viscosity used on this processor

examples/IMEX/drift-wave-constraint/test-drift.cxx

@@ -82,7 +82,7 @@ class DriftWave : public PhysicsModel {
     // ddt(phi) = Delp2(phi) - Vort;
 
     // This version uses central differencing for Delp2
-    ddt(phi) = (coord->g11 * D2DX2(phi) + coord->g33 * D2DZ2(phi)) - Vort;
+    ddt(phi) = (coord->g11() * D2DX2(phi) + coord->g33() * D2DZ2(phi)) - Vort;
 
     return 0;
   }

examples/conducting-wall-mode/cwm.cxx

@@ -65,7 +65,7 @@ class CWM : public PhysicsModel {
 
     // Load metrics
     GRID_LOAD(Rxy, Zxy, Bpxy, Btxy, hthe);
-    mesh->get(coord->dx, "dpsi");
+    coord->setDx(mesh->get("dpsi"));
     mesh->get(I, "sinty");
 
     // Load normalisation values
@@ -136,35 +136,36 @@ class CWM : public PhysicsModel {
     Rxy /= rho_s;
     hthe /= rho_s;
     I *= rho_s * rho_s * (bmag / 1e4) * ShearFactor;
-    coord->dx /= rho_s * rho_s * (bmag / 1e4);
+    coord->setDx(coord->dx() / (rho_s * rho_s * (bmag / 1e4)));
 
     // Normalise magnetic field
     Bpxy /= (bmag / 1.e4);
     Btxy /= (bmag / 1.e4);
-    coord->Bxy /= (bmag / 1.e4);
+    coord->setBxy(coord->Bxy() / (bmag / 1.e4));
 
     // Set nu
     nu = nu_hat * Ni0 / pow(Te0, 1.5);
 
     /**************** CALCULATE METRICS ******************/
 
-    coord->g11 = SQ(Rxy * Bpxy);
-    coord->g22 = 1.0 / SQ(hthe);
-    coord->g33 = SQ(I) * coord->g11 + SQ(coord->Bxy) / coord->g11;
-    coord->g12 = 0.0;
-    coord->g13 = -I * coord->g11;
-    coord->g23 = -Btxy / (hthe * Bpxy * Rxy);
+    const auto g11 = SQ(Rxy * Bpxy);
+    const auto g22 = 1.0 / SQ(hthe);
+    const auto g33 = SQ(I) * g11 + SQ(coord->Bxy()) / g11;
+    const auto g12 = 0.0;
+    const auto g13 = -I * g11;
+    const auto g23 = -Btxy / (hthe * Bpxy * Rxy);
 
-    coord->J = hthe / Bpxy;
+    const auto g_11 = 1.0 / g11 + SQ(I * Rxy);
+    const auto g_22 = SQ(coord->Bxy() * hthe / Bpxy);
+    const auto g_33 = Rxy * Rxy;
+    const auto g_12 = Btxy * hthe * I * Rxy / Bpxy;
+    const auto g_13 = I * Rxy * Rxy;
+    const auto g_23 = Btxy * hthe * Rxy / Bpxy;
 
-    coord->g_11 = 1.0 / coord->g11 + SQ(I * Rxy);
-    coord->g_22 = SQ(coord->Bxy * hthe / Bpxy);
-    coord->g_33 = Rxy * Rxy;
-    coord->g_12 = Btxy * hthe * I * Rxy / Bpxy;
-    coord->g_13 = I * Rxy * Rxy;
-    coord->g_23 = Btxy * hthe * Rxy / Bpxy;
+    coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
+                           CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));
 
-    coord->geometry();
+    coord->setJ(hthe / Bpxy);
 
     /**************** SET EVOLVING VARIABLES *************/
 
@@ -349,7 +350,7 @@ class CWM : public PhysicsModel {
       result = VDDX(DDZ(p), f);
     } else {
       // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
     }
     return result;
   }
@@ -361,7 +362,7 @@ class CWM : public PhysicsModel {
       result = VDDZ(-DDX(p), f);
     } else {
       // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
     }
     return result;
   }
@@ -373,7 +374,7 @@ class CWM : public PhysicsModel {
       result = VDDX(DDZ(p), f) + VDDZ(-DDX(p), f);
     } else {
       // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
     }
     return result;
   }

examples/constraints/alfven-wave/alfven.cxx

@@ -170,7 +170,7 @@ class Alfven : public PhysicsModel {
     Field2D dx;
     if (!mesh->get(dx, "dpsi")) {
       output << "\tUsing dpsi as the x grid spacing\n";
-      coord->dx = dx; // Only use dpsi if found
+      coord->setDx(dx); // Only use dpsi if found
     } else {
       // dx will have been read already from the grid
       output << "\tUsing dx as the x grid spacing\n";
@@ -179,11 +179,11 @@ class Alfven : public PhysicsModel {
     Rxy /= Lnorm;
     hthe /= Lnorm;
     sinty *= SQ(Lnorm) * Bnorm;
-    coord->dx /= SQ(Lnorm) * Bnorm;
+    coord->setDx(coord->dx() / (SQ(Lnorm) * Bnorm));
 
     Bpxy /= Bnorm;
     Btxy /= Bnorm;
-    coord->Bxy /= Bnorm;
+    coord->setBxy(coord->Bxy() / Bnorm);
 
     // Check type of parallel transform
     std::string ptstr =
@@ -201,23 +201,24 @@ class Alfven : public PhysicsModel {
 
     // Calculate metric components
 
-    coord->g11 = SQ(Rxy * Bpxy);
-    coord->g22 = 1.0 / SQ(hthe);
-    coord->g33 = SQ(sinty) * coord->g11 + SQ(coord->Bxy) / coord->g11;
-    coord->g12 = 0.0;
-    coord->g13 = -sinty * coord->g11;
-    coord->g23 = -sbp * Btxy / (hthe * Bpxy * Rxy);
+    const auto g11 = SQ(Rxy * Bpxy);
+    const auto g22 = 1.0 / SQ(hthe);
+    const auto g33 = SQ(sinty) * g11 + SQ(coord->Bxy()) / g11;
+    const auto g12 = 0.0;
+    const auto g13 = -sinty * g11;
+    const auto g23 = -sbp * Btxy / (hthe * Bpxy * Rxy);
 
-    coord->J = hthe / Bpxy;
+    const auto g_11 = 1.0 / g11 + SQ(sinty * Rxy);
+    const auto g_22 = SQ(coord->Bxy() * hthe / Bpxy);
+    const auto g_33 = Rxy * Rxy;
+    const auto g_12 = sbp * Btxy * hthe * sinty * Rxy / Bpxy;
+    const auto g_13 = sinty * Rxy * Rxy;
+    const auto g_23 = sbp * Btxy * hthe * Rxy / Bpxy;
 
-    coord->g_11 = 1.0 / coord->g11 + SQ(sinty * Rxy);
-    coord->g_22 = SQ(coord->Bxy * hthe / Bpxy);
-    coord->g_33 = Rxy * Rxy;
-    coord->g_12 = sbp * Btxy * hthe * sinty * Rxy / Bpxy;
-    coord->g_13 = sinty * Rxy * Rxy;
-    coord->g_23 = sbp * Btxy * hthe * Rxy / Bpxy;
+    coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
+                           CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));
 
-    coord->geometry();
+    coord->setJ(hthe / Bpxy);
   }
 };
 

examples/dalf3/dalf3.cxx

@@ -243,29 +243,30 @@ class DALF3 : public PhysicsModel {
     Btxy /= Bnorm;
     B0 /= Bnorm;
 
-    coord->dx /= rho_s * rho_s * Bnorm;
+    coord->setDx(coord->dx() / (rho_s * rho_s * Bnorm));
 
     ///////////////////////////////////////////////////
     // CALCULATE METRICS
 
-    coord->g11 = SQ(Rxy * Bpxy);
-    coord->g22 = 1.0 / SQ(hthe);
-    coord->g33 = SQ(I) * coord->g11 + SQ(B0) / coord->g11;
-    coord->g12 = 0.0;
-    coord->g13 = -I * coord->g11;
-    coord->g23 = -Btxy / (hthe * Bpxy * Rxy);
-
-    coord->J = hthe / Bpxy;
-    coord->Bxy = B0;
-
-    coord->g_11 = 1.0 / coord->g11 + SQ(I * Rxy);
-    coord->g_22 = SQ(B0 * hthe / Bpxy);
-    coord->g_33 = Rxy * Rxy;
-    coord->g_12 = Btxy * hthe * I * Rxy / Bpxy;
-    coord->g_13 = I * Rxy * Rxy;
-    coord->g_23 = Btxy * hthe * Rxy / Bpxy;
-
-    coord->geometry(); // Calculate quantities from metric tensor
+    const auto g11 = SQ(Rxy * Bpxy);
+    const auto g22 = 1.0 / SQ(hthe);
+    const auto g33 = SQ(I) * g11 + SQ(B0) / g11;
+    const auto g12 = 0.0;
+    const auto g13 = -I * g11;
+    const auto g23 = -Btxy / (hthe * Bpxy * Rxy);
+
+    const auto g_11 = 1.0 / g11 + SQ(I * Rxy);
+    const auto g_22 = SQ(B0 * hthe / Bpxy);
+    const auto g_33 = Rxy * Rxy;
+    const auto g_12 = Btxy * hthe * I * Rxy / Bpxy;
+    const auto g_13 = I * Rxy * Rxy;
+    const auto g_23 = Btxy * hthe * Rxy / Bpxy;
+
+    coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
+                           CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));
+
+    coord->setJ(hthe / Bpxy);
+    coord->setBxy(B0);
 
     SOLVE_FOR3(Vort, Pe, Vpar);
     comms.add(Vort, Pe, Vpar);

examples/elm-pb-outerloop/elm_pb_outerloop.cxx

@@ -825,7 +825,7 @@ class ELMpb : public PhysicsModel {
 
     if (mesh->IncIntShear) {
       // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
-      metric->IntShiftTorsion = I;
+      metric->setIntShiftTorsion(I);
 
     } else {
       // Dimits style, using local coordinate system
@@ -962,7 +962,7 @@ class ELMpb : public PhysicsModel {
     Btxy /= Bbar;
     B0 /= Bbar;
     hthe /= Lbar;
-    metric->dx /= Lbar * Lbar * Bbar;
+    metric->setDx(metric->dx() / (Lbar * Lbar * Bbar));
     I *= Lbar * Lbar * Bbar;
 
     if (constn0) {
@@ -1092,24 +1092,25 @@ class ELMpb : public PhysicsModel {
 
     /**************** CALCULATE METRICS ******************/
 
-    metric->g11 = SQ(Rxy * Bpxy);
-    metric->g22 = 1.0 / SQ(hthe);
-    metric->g33 = SQ(I) * metric->g11 + SQ(B0) / metric->g11;
-    metric->g12 = 0.0;
-    metric->g13 = -I * metric->g11;
-    metric->g23 = -Btxy / (hthe * Bpxy * Rxy);
+    const auto g11 = SQ(Rxy * Bpxy);
+    const auto g22 = 1.0 / SQ(hthe);
+    const auto g33 = SQ(I) * g11 + SQ(B0) / g11;
+    const auto g12 = 0.0;
+    const auto g13 = -I * g11;
+    const auto g23 = -Btxy / (hthe * Bpxy * Rxy);
 
-    metric->J = hthe / Bpxy;
-    metric->Bxy = B0;
+    const auto g_11 = 1.0 / g11 + SQ(I * Rxy);
+    const auto g_22 = SQ(B0 * hthe / Bpxy);
+    const auto g_33 = Rxy * Rxy;
+    const auto g_12 = Btxy * hthe * I * Rxy / Bpxy;
+    const auto g_13 = I * Rxy * Rxy;
+    const auto g_23 = Btxy * hthe * Rxy / Bpxy;
 
-    metric->g_11 = 1.0 / metric->g11 + SQ(I * Rxy);
-    metric->g_22 = SQ(B0 * hthe / Bpxy);
-    metric->g_33 = Rxy * Rxy;
-    metric->g_12 = Btxy * hthe * I * Rxy / Bpxy;
-    metric->g_13 = I * Rxy * Rxy;
-    metric->g_23 = Btxy * hthe * Rxy / Bpxy;
+    metric->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
+                            CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));
 
-    metric->geometry(); // Calculate quantities from metric tensor
+    metric->setJ(hthe / Bpxy);
+    metric->setBxy(B0);
 
     // Set B field vector
 
@@ -1792,11 +1793,11 @@ class ELMpb : public PhysicsModel {
     if (hyperviscos > 0.0) {
       // Calculate coefficient.
 
-      hyper_mu_x = hyperviscos * metric->g_11 * SQ(metric->dx)
-                   * abs(metric->g11 * D2DX2(U)) / (abs(U) + 1e-3);
+      hyper_mu_x = hyperviscos * metric->g_11() * SQ(metric->dx())
+                   * abs(metric->g11() * D2DX2(U)) / (abs(U) + 1e-3);
       hyper_mu_x.applyBoundary("dirichlet"); // Set to zero on all boundaries
 
-      ddt(U) += hyper_mu_x * metric->g11 * D2DX2(U);
+      ddt(U) += hyper_mu_x * metric->g11() * D2DX2(U);
 
       if (first_run) { // Print out maximum values of viscosity used on this processor
         output.write("   Hyper-viscosity values:\n");
@@ -1907,7 +1908,7 @@ class ELMpb : public PhysicsModel {
         BoutReal pnorm = P0(0, 0);
         ddt(P) += heating_P * source_expx2(P0, 2. * hp_width, 0.5 * hp_length)
                   * (Tbar / pnorm); // heat source
-        ddt(P) += (100. * source_tanhx(P0, hp_width, hp_length) + 0.01) * metric->g11
+        ddt(P) += (100. * source_tanhx(P0, hp_width, hp_length) + 0.01) * metric->g11()
                   * D2DX2(P) * (Tbar / Lbar / Lbar); // radial diffusion
       }