Remove b0xcv (curvature term) from TokamakCoordinatesFactory

tomc271 · Nov 22, 2024 · f2a0ffd · f2a0ffd
1 parent f818a36
commit f2a0ffd
Show file tree

Hide file tree

Showing 19 changed files with 210 additions and 75 deletions.
diff --git a/examples/6field-simple/elm_6f.cxx b/examples/6field-simple/elm_6f.cxx
@@ -36,7 +36,8 @@ class Elm_6f : public PhysicsModel {
   // 2D inital profiles
   /// Current and pressure
   Field2D J0, P0;
-
+  /// Curvature term
+  Vector2D b0xcv;
   /// When diamagnetic terms used
   Field2D phi0;
 
@@ -388,6 +389,10 @@ class Elm_6f : public PhysicsModel {
     mesh->get(J0, "Jpar0");    // A / m^2
     mesh->get(P0, "pressure"); // Pascals
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     //////////////////////////////////////////////////////////////
     // Read parameters from the options file
     //
@@ -663,10 +668,36 @@ class Elm_6f : public PhysicsModel {
     phi_curv = options["phi_curv"].doc("Compressional ExB terms").withDefault(true);
     g = options["gamma"].doc("Ratio of specific heats").withDefault(5.0 / 3.0);
 
+    if (!include_curvature) {
+      b0xcv = 0.0;
+    }
+
     if (!include_jpar0) {
       J0 = 0.0;
     }
 
+    if (noshear) {
+      if (include_curvature) {
+        b0xcv.z += I * b0xcv.x;
+      }
+      I = 0.0;
+    }
+
+    //////////////////////////////////////////////////////////////
+    // SHIFTED RADIAL COORDINATES
+
+    if (mesh->IncIntShear) {
+      // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
+      coord->setIntShiftTorsion(I);
+
+    } else {
+      // Dimits style, using local coordinate system
+      if (include_curvature) {
+        b0xcv.z += I * b0xcv.x;
+      }
+      I = 0.0; // I disappears from metric
+    }
+
     //////////////////////////////////////////////////////////////
     // NORMALISE QUANTITIES
 
@@ -812,6 +843,10 @@ class Elm_6f : public PhysicsModel {
 
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
 
+    b0xcv.x /= Bbar;
+    b0xcv.y *= Lbar * Lbar;
+    b0xcv.z *= Lbar * Lbar;
+
     if ((!T0_fake_prof) && n0_fake_prof) {
       N0 = N0tanh(n0_height * Nbar, n0_ave * Nbar, n0_width, n0_center, n0_bottom_x);
 
@@ -999,7 +1034,7 @@ class Elm_6f : public PhysicsModel {
     }
 
     /**************** CALCULATE METRICS ******************/
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, include_curvature);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     //////////////////////////////////////////////////////////////
     // SHIFTED RADIAL COORDINATES
@@ -1358,7 +1393,7 @@ class Elm_6f : public PhysicsModel {
 
       ddt(U) = -SQ(tokamak_coordinates_factory.get_Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates_factory.get_Bxy(); // Grad j term
 
-      ddt(U) += 2.0 * Upara1 * tokamak_coordinates_factory.get_b0xcv() * Grad(P); // curvature term
+      ddt(U) += 2.0 * Upara1 * b0xcv * Grad(P); // curvature term
 
       ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar); // b dot grad j
 

diff --git a/examples/conducting-wall-mode/cwm.cxx b/examples/conducting-wall-mode/cwm.cxx
@@ -130,7 +130,7 @@ class CWM : public PhysicsModel {
 
     // Normalise geometry
     tokamak_coordinates_factory.normalise(rho_s, bmag / 1e4, ShearFactor);
-    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     // Set nu
     nu = nu_hat * Ni0 / pow(Te0, 1.5);

diff --git a/examples/constraints/alfven-wave/alfven.cxx b/examples/constraints/alfven-wave/alfven.cxx
@@ -172,7 +172,7 @@ class Alfven : public PhysicsModel {
 
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
     tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
   }
 };
 

diff --git a/examples/dalf3/dalf3.cxx b/examples/dalf3/dalf3.cxx
@@ -53,6 +53,7 @@ class DALF3 : public PhysicsModel {
   Field3D phi, apar, jpar;
 
   Field2D B0, Pe0, Jpar0;
+  Vector2D b0xcv;
 
   Field2D eta; // Collisional damping (resistivity)
   BoutReal beta_hat, mu_hat;
@@ -99,6 +100,10 @@ class DALF3 : public PhysicsModel {
     Pe0 = 2. * Charge * Ni0 * Te0; // Electron pressure in Pascals
     SAVE_ONCE(Pe0);
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     //////////////////////////////////////////////////////////////
     // Options
 
@@ -156,7 +161,9 @@ class DALF3 : public PhysicsModel {
         Options::root()["mesh"]["paralleltransform"]["type"].withDefault("identity");
 
     if (lowercase(ptstr) == "shifted") {
-      noshear = true;
+      // Dimits style, using local coordinate system
+      b0xcv.z += I * b0xcv.x;
+      I = 0.0; // I disappears from metric
     }
 
     ///////////////////////////////////////////////////
@@ -210,9 +217,13 @@ class DALF3 : public PhysicsModel {
     hyper_viscosity /= wci * SQ(SQ(rho_s));
     viscosity_par /= wci * SQ(rho_s);
 
+    b0xcv.x /= Bnorm;
+    b0xcv.y *= rho_s * rho_s;
+    b0xcv.z *= rho_s * rho_s;
+
     // Metrics
     tokamak_coordinates_factory.normalise(rho_s, Bnorm);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     SOLVE_FOR3(Vort, Pe, Vpar);
     comms.add(Vort, Pe, Vpar);

diff --git a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx
@@ -96,6 +96,7 @@ class ELMpb : public PhysicsModel {
 private:
   // 2D inital profiles
   Field2D J0, P0; // Current and pressure
+  Vector2D b0xcv; // Curvature term
   Field2D beta;   // Used for Vpar terms
   Coordinates::FieldMetric gradparB;
   Field2D phi0; // When diamagnetic terms used
@@ -368,6 +369,10 @@ class ELMpb : public PhysicsModel {
     mesh->get(J0, "Jpar0");    // A / m^2
     mesh->get(P0, "pressure"); // Pascals
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     mesh->get(Psixy, "psixy");        // get Psi
     mesh->get(Psiaxis, "psi_axis");   // axis flux
     mesh->get(Psibndry, "psi_bndry"); // edge flux
@@ -712,7 +717,7 @@ class ELMpb : public PhysicsModel {
       Lbar = 1.0;
     }
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
-    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, include_curvature);
+    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
     V0 = -tokamak_coordinates_factory.get_Rxy() * tokamak_coordinates_factory.get_Bpxy() * Dphi0 / tokamak_coordinates_factory.get_Bxy();
 
@@ -799,10 +804,34 @@ class ELMpb : public PhysicsModel {
       }
     }
 
+    if (!include_curvature) {
+      b0xcv = 0.0;
+    }
+
     if (!include_jpar0) {
       J0 = 0.0;
     }
 
+    if (noshear) {
+      if (include_curvature) {
+        b0xcv.z += I * b0xcv.x;
+      }
+      I = 0.0;
+    }
+
+    //////////////////////////////////////////////////////////////
+    // SHIFTED RADIAL COORDINATES
+
+    if (not mesh->IncIntShear) {
+      // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
+      metric->setIntShiftTorsion(I);
+      // Dimits style, using local coordinate system
+      if (include_curvature) {
+        b0xcv.z += I * b0xcv.x;
+      }
+      I = 0.0; // I disappears from metric
+    }
+
     //////////////////////////////////////////////////////////////
     // NORMALISE QUANTITIES
 
@@ -914,6 +943,11 @@ class ELMpb : public PhysicsModel {
     V0 = V0 / Va;
     Dphi0 *= Tbar;
 
+    b0xcv.x /= Bbar;
+    b0xcv.y *= Lbar * Lbar;
+    b0xcv.z *= Lbar * Lbar;
+
+
     if (constn0) {
       T0_fake_prof = false;
       n0_fake_prof = false;
@@ -1690,7 +1724,7 @@ class ELMpb : public PhysicsModel {
     //   ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
     // }
 
-    ddt(U) += tokamak_coordinates_factory.get_b0xcv() * Grad(P); // curvature term
+    ddt(U) += b0xcv * Grad(P); // curvature term
 
     // if (!nogradparj) { // Parallel current term
     //   ddt(U) -= SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
@@ -1866,7 +1900,7 @@ class ELMpb : public PhysicsModel {
       ddt(P) -= beta * Div_par(Vpar, CELL_CENTRE);
 
       if (phi_curv) {
-        ddt(P) -= 2. * beta * tokamak_coordinates_factory.get_b0xcv() * Grad(phi);
+        ddt(P) -= 2. * beta * b0xcv * Grad(phi);
       }
 
       // Vpar equation
@@ -1966,7 +2000,7 @@ class ELMpb : public PhysicsModel {
     Coordinates* metric = mesh->getCoordinates();
 
     Field3D U1 = ddt(U);
-    U1 += (gamma * metric->Bxy() * metric->Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * tokamak_coordinates_factory.get_b0xcv()) * Grad(P);
+    U1 += (gamma * metric->Bxy() * metric->Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P);
 
     // Second matrix, solving Alfven wave dynamics
     static std::unique_ptr<InvertPar> invU{nullptr};
@@ -2025,7 +2059,7 @@ class ELMpb : public PhysicsModel {
     JPsi.setBoundary("Psi");
     JPsi.applyBoundary();
 
-    Field3D JU = tokamak_coordinates_factory.get_b0xcv() * Grad(ddt(P)) - SQ(metric->Bxy()) * Grad_par(Jpar, CELL_CENTRE)
+    Field3D JU = b0xcv * Grad(ddt(P)) - SQ(metric->Bxy()) * Grad_par(Jpar, CELL_CENTRE)
                  + SQ(metric->Bxy()) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE);
     JU.setBoundary("U");
     JU.applyBoundary();

diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx
@@ -38,6 +38,7 @@ class ELMpb : public PhysicsModel {
 private:
   // 2D inital profiles
   Field2D J0, P0; // Current and pressure
+  Vector2D b0xcv; // Curvature term
   Field2D beta;   // Used for Vpar terms
   Coordinates::FieldMetric gradparB;
   Field2D phi0; // When diamagnetic terms used
@@ -310,6 +311,10 @@ class ELMpb : public PhysicsModel {
     mesh->get(J0, "Jpar0");    // A / m^2
     mesh->get(P0, "pressure"); // Pascals
 
+    // Load curvature term
+    b0xcv.covariant = false;  // Read contravariant components
+    mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2
+
     mesh->get(Psixy, "psixy");        // get Psi
     mesh->get(Psiaxis, "psi_axis");   // axis flux
     mesh->get(Psibndry, "psi_bndry"); // edge flux
@@ -725,13 +730,28 @@ class ELMpb : public PhysicsModel {
       }
     }
 
+    if (!include_curvature) {
+      b0xcv = 0.0;
+    }
+
     if (!include_jpar0) {
       J0 = 0.0;
     }
 
-    // Dimits style, using local coordinate system
-    if (include_curvature and !mesh->IncIntShear) {
-      noshear = true;
+    if (noshear) {
+      if (include_curvature) {
+        b0xcv.z += I * b0xcv.x;
+      }
+    }
+
+    //////////////////////////////////////////////////////////////
+    // SHIFTED RADIAL COORDINATES
+
+    if (not mesh->IncIntShear) {
+      // Dimits style, using local coordinate system
+      if (include_curvature) {
+        b0xcv.z += I * b0xcv.x;
+      }
     }
 
     //////////////////////////////////////////////////////////////
@@ -852,8 +872,12 @@ class ELMpb : public PhysicsModel {
     V0 = V0 / Va;
     Dphi0 *= Tbar;
 
+    b0xcv.x /= Bbar;
+    b0xcv.y *= Lbar * Lbar;
+    b0xcv.z *= Lbar * Lbar;
+
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
-    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, include_curvature);
+    const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
 
 
     //////////////////////////////////////////////////////////////
@@ -1499,7 +1523,7 @@ class ELMpb : public PhysicsModel {
       ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
     }
 
-    ddt(U) += tokamak_coordinates_factory.get_b0xcv() * Grad(P); // curvature term
+    ddt(U) += b0xcv * Grad(P); // curvature term
 
     if (!nogradparj) {                                 // Parallel current term
       ddt(U) -= SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
@@ -1673,7 +1697,7 @@ class ELMpb : public PhysicsModel {
       ddt(P) -= beta * Div_par(Vpar, CELL_CENTRE);
 
       if (phi_curv) {
-        ddt(P) -= 2. * beta * tokamak_coordinates_factory.get_b0xcv() * Grad(phi);
+        ddt(P) -= 2. * beta * b0xcv * Grad(phi);
       }
 
       // Vpar equation
@@ -1772,7 +1796,7 @@ class ELMpb : public PhysicsModel {
     mesh->communicate(Jrhs, ddt(P));
 
     Field3D U1 = ddt(U);
-    U1 += (gamma * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * tokamak_coordinates_factory.get_b0xcv()) * Grad(P);
+    U1 += (gamma * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P);
 
     // Second matrix, solving Alfven wave dynamics
     static std::unique_ptr<InvertPar> invU{nullptr};
@@ -1829,7 +1853,7 @@ class ELMpb : public PhysicsModel {
     JPsi.setBoundary("Psi");
     JPsi.applyBoundary();
 
-    Field3D JU = tokamak_coordinates_factory.get_b0xcv() * Grad(ddt(P)) - SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jpar, CELL_CENTRE)
+    Field3D JU = b0xcv * Grad(ddt(P)) - SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jpar, CELL_CENTRE)
                  + SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE);
     JU.setBoundary("U");
     JU.applyBoundary();

diff --git a/examples/gyro-gem/gem.cxx b/examples/gyro-gem/gem.cxx
@@ -253,7 +253,7 @@ class GEM : public PhysicsModel {
     SAVE_ONCE(Bbar);
 
     tokamak_coordinates_factory.normalise(Lbar, Bbar);
-    coord = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    coord = tokamak_coordinates_factory.make_tokamak_coordinates(true);
 
     beta_e = 4.e-7 * PI * max(p_e, true) / (Bbar * Bbar);
     SAVE_ONCE(beta_e);

diff --git a/examples/laplacexy/alfven-wave/alfven.cxx b/examples/laplacexy/alfven-wave/alfven.cxx
@@ -174,7 +174,7 @@ class Alfven : public PhysicsModel {
 
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
     tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true);
   }
 };
 

diff --git a/examples/laplacexy/laplace_perp/test.cxx b/examples/laplacexy/laplace_perp/test.cxx
@@ -15,7 +15,7 @@ int main(int argc, char** argv) {
   calc_metric = Options::root()["calc_metric"].withDefault(false);
   if (calc_metric) {
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
-    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(true);
   }
 
   ///////////////////////////////////////

diff --git a/examples/wave-slab/wave_slab.cxx b/examples/wave-slab/wave_slab.cxx
@@ -18,7 +18,7 @@ class WaveTest : public PhysicsModel {
   int init(bool UNUSED(restarting)) {
 
     auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
-    const auto& coords = tokamak_coordinates_factory.make_tokamak_coordinates(true, true);
+    const auto& coords = tokamak_coordinates_factory.make_tokamak_coordinates(true);
 
     solver->add(f, "f");
     solver->add(g, "g");