From fa35e6a5a5f09032773db3b851a6d2b07e51cef4 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Mon, 16 Dec 2024 18:18:10 +0000 Subject: [PATCH] Rename class and method class `TokamakCoordinatesFactory` renamed to `TokamakCoordinates` method `tokamak_coordinates_factory` renamed to `tokamak_coordinates` --- CMakeLists.txt | 2 +- examples/6field-simple/elm_6f.cxx | 99 ++++++++------- examples/conducting-wall-mode/cwm.cxx | 16 +-- examples/constraints/alfven-wave/alfven.cxx | 6 +- examples/dalf3/dalf3.cxx | 8 +- .../elm-pb-outerloop/elm_pb_outerloop.cxx | 41 +++---- examples/elm-pb/elm_pb.cxx | 113 +++++++++--------- examples/gyro-gem/gem.cxx | 10 +- examples/laplacexy/alfven-wave/alfven.cxx | 7 +- examples/laplacexy/laplace_perp/test.cxx | 6 +- examples/wave-slab/wave_slab.cxx | 6 +- ...es_factory.hxx => tokamak_coordinates.hxx} | 12 +- tests/MMS/GBS/gbs.cxx | 4 +- tests/MMS/elm-pb/elm_pb.cxx | 8 +- tests/MMS/tokamak/tokamak.cxx | 4 +- .../test-drift-instability/2fluid.cxx | 9 +- .../test-interchange-instability/2fluid.cxx | 10 +- .../integrated/test-laplacexy/loadmetric.cxx | 4 +- 18 files changed, 182 insertions(+), 183 deletions(-) rename include/bout/{tokamak_coordinates_factory.hxx => tokamak_coordinates.hxx} (88%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6cb42c2f5a..c561e91fa5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -360,7 +360,7 @@ set(BOUT_SOURCES ./src/mesh/g_values.cxx ${CMAKE_CURRENT_BINARY_DIR}/include/bout/revision.hxx ${CMAKE_CURRENT_BINARY_DIR}/include/bout/version.hxx - ./include/bout/tokamak_coordinates_factory.hxx + ./include/bout/tokamak_coordinates.hxx ) diff --git a/examples/6field-simple/elm_6f.cxx b/examples/6field-simple/elm_6f.cxx index 6a51fc4221..88bf18b4b0 100644 --- a/examples/6field-simple/elm_6f.cxx +++ b/examples/6field-simple/elm_6f.cxx @@ -15,7 +15,7 @@ #include "bout/msg_stack.hxx" #include "bout/physicsmodel.hxx" #include "bout/sourcex.hxx" -#include "bout/tokamak_coordinates_factory.hxx" +#include "bout/tokamak_coordinates.hxx" #include @@ -235,7 +235,7 @@ class Elm_6f : public PhysicsModel { int damp_width; // Width of inner damped region BoutReal damp_t_const; // Timescale of damping - TokamakCoordinatesFactory tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + TokamakCoordinates tokamak_coordinates = TokamakCoordinates(*mesh); BoutReal LnLambda; // ln(Lambda) @@ -366,7 +366,7 @@ class Elm_6f : public PhysicsModel { result = Grad_par(f, loc); if (nonlinear) { - result -= bracket(Psi, f, bm_mag) * tokamak_coordinates_factory.Bxy(); + result -= bracket(Psi, f, bm_mag) * tokamak_coordinates.Bxy(); } } @@ -678,7 +678,7 @@ class Elm_6f : public PhysicsModel { if (noshear) { if (include_curvature) { - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; } } @@ -688,7 +688,7 @@ class Elm_6f : public PhysicsModel { if (not mesh->IncIntShear) { // Dimits style, using local coordinate system if (include_curvature) { - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; } } @@ -832,7 +832,7 @@ class Elm_6f : public PhysicsModel { dump.add(sp_length, "sp_length", 1); } - J0 = SI::mu0 * Lbar * J0 / tokamak_coordinates_factory.Bxy(); + J0 = SI::mu0 * Lbar * J0 / tokamak_coordinates.Bxy(); P0 = P0 / (SI::kb * (Tibar + Tebar) * eV_K / 2. * Nbar * density); b0xcv.x /= Bbar; @@ -904,9 +904,9 @@ class Elm_6f : public PhysicsModel { q95 = q95_input; // use a constant for test } else { if (local_q) { - q95 = abs(tokamak_coordinates_factory.hthe() - * tokamak_coordinates_factory.Btxy() - / tokamak_coordinates_factory.Bpxy()) + q95 = abs(tokamak_coordinates.hthe() + * tokamak_coordinates.Btxy() + / tokamak_coordinates.Bpxy()) * q_alpha; } else { output.write("\tUsing q profile from grid.\n"); @@ -1029,37 +1029,36 @@ class Elm_6f : public PhysicsModel { } /**************** CALCULATE METRICS ******************/ - const auto& coord = - tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lbar, Bbar); + const auto& coord = tokamak_coordinates.make_coordinates(noshear, Lbar, Bbar); ////////////////////////////////////////////////////////////// // SHIFTED RADIAL COORDINATES if (mesh->IncIntShear) { // BOUT-06 style, using d/dx = d/dpsi + I * d/dz - coord->setIntShiftTorsion(tokamak_coordinates_factory.ShearFactor()); + coord->setIntShiftTorsion(tokamak_coordinates.ShearFactor()); } // Set B field vector B0vec.covariant = false; B0vec.x = 0.; - B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe(); + B0vec.y = tokamak_coordinates.Bpxy() / tokamak_coordinates.hthe(); B0vec.z = 0.; // Set V0vec field vector V0vec.covariant = false; V0vec.x = 0.; - V0vec.y = Vp0 / tokamak_coordinates_factory.hthe(); - V0vec.z = Vt0 / tokamak_coordinates_factory.Rxy(); + V0vec.y = Vp0 / tokamak_coordinates.hthe(); + V0vec.z = Vt0 / tokamak_coordinates.Rxy(); // Set V0eff field vector V0eff.covariant = false; V0eff.x = 0.; - V0eff.y = -(tokamak_coordinates_factory.Btxy() / (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy())) * (Vp0 * tokamak_coordinates_factory.Btxy() - Vt0 * tokamak_coordinates_factory.Bpxy()) / tokamak_coordinates_factory.hthe(); - V0eff.z = (tokamak_coordinates_factory.Bpxy() / (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy())) * (Vp0 * tokamak_coordinates_factory.Btxy() - Vt0 * tokamak_coordinates_factory.Bpxy()) / tokamak_coordinates_factory.Rxy(); + V0eff.y = -(tokamak_coordinates.Btxy() / (tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy())) * (Vp0 * tokamak_coordinates.Btxy() - Vt0 * tokamak_coordinates.Bpxy()) / tokamak_coordinates.hthe(); + V0eff.z = (tokamak_coordinates.Bpxy() / (tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy())) * (Vp0 * tokamak_coordinates.Btxy() - Vt0 * tokamak_coordinates.Bpxy()) / tokamak_coordinates.Rxy(); Pe.setBoundary("P"); Pi.setBoundary("P"); @@ -1101,10 +1100,10 @@ class Elm_6f : public PhysicsModel { if (diamag && diamag_phi0) { if (experiment_Er) { // get phi0 from grid file mesh->get(phi0, "Phi_0"); - phi0 /= tokamak_coordinates_factory.Bxy() * Lbar * Va; + phi0 /= tokamak_coordinates.Bxy() * Lbar * Va; } else { // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift - phi0 = -Upara0 * Pi0 / tokamak_coordinates_factory.Bxy() / N0; + phi0 = -Upara0 * Pi0 / tokamak_coordinates.Bxy() / N0; } SAVE_ONCE(phi0); } @@ -1114,7 +1113,7 @@ class Elm_6f : public PhysicsModel { SAVE_ONCE(J0, P0); SAVE_ONCE(density, Lbar, Bbar, Tbar); SAVE_ONCE(Tibar, Tebar, Nbar); - Field2D tmp = tokamak_coordinates_factory.Bxy(); + Field2D tmp = tokamak_coordinates.Bxy(); SAVE_ONCE(Va, tmp); SAVE_ONCE(Ti0, Te0, N0); @@ -1138,13 +1137,13 @@ class Elm_6f : public PhysicsModel { Field2D logn0 = laplace_alpha * N0; Field3D Ntemp; Ntemp = N0; - ubyn = U * tokamak_coordinates_factory.Bxy() / Ntemp; + ubyn = U * tokamak_coordinates.Bxy() / Ntemp; // Phi should be consistent with U if (laplace_alpha <= 0.0) { - phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy(); + phi = phiSolver->solve(ubyn) / tokamak_coordinates.Bxy(); } else { phiSolver->setCoefC(logn0); - phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy(); + phi = phiSolver->solve(ubyn) / tokamak_coordinates.Bxy(); } } @@ -1219,9 +1218,9 @@ class Elm_6f : public PhysicsModel { // Field2D lap_temp=0.0; Field2D logn0 = laplace_alpha * N0; - ubyn = U * tokamak_coordinates_factory.Bxy() / N0; + ubyn = U * tokamak_coordinates.Bxy() / N0; if (diamag) { - ubyn -= Upara0 / N0 * Delp2(Pi) / tokamak_coordinates_factory.Bxy(); + ubyn -= Upara0 / N0 * Delp2(Pi) / tokamak_coordinates.Bxy(); mesh->communicate(ubyn); ubyn.applyBoundary(); } @@ -1229,7 +1228,7 @@ class Elm_6f : public PhysicsModel { if (laplace_alpha > 0.0) { phiSolver->setCoefC(logn0); } - phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy(); + phi = phiSolver->solve(ubyn) / tokamak_coordinates.Bxy(); mesh->communicate(phi); @@ -1325,9 +1324,9 @@ class Elm_6f : public PhysicsModel { if (compress0) { if (nonlinear) { - Vepar = Vipar - tokamak_coordinates_factory.Bxy() * (Jpar) / N_tmp * Vepara; + Vepar = Vipar - tokamak_coordinates.Bxy() * (Jpar) / N_tmp * Vepara; } else { - Vepar = Vipar - tokamak_coordinates_factory.Bxy() * (Jpar) / N0 * Vepara; + Vepar = Vipar - tokamak_coordinates.Bxy() * (Jpar) / N0 * Vepara; Vepar.applyBoundary(); mesh->communicate(Vepar); } @@ -1365,13 +1364,13 @@ class Elm_6f : public PhysicsModel { ddt(Psi) = 0.0; if (spitzer_resist) { - ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.Bxy() * phi) / tokamak_coordinates_factory.Bxy() - eta_spitzer * Jpar; + ddt(Psi) = -Grad_parP(tokamak_coordinates.Bxy() * phi) / tokamak_coordinates.Bxy() - eta_spitzer * Jpar; } else { - ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.Bxy() * phi) / tokamak_coordinates_factory.Bxy() - eta * Jpar; + ddt(Psi) = -Grad_parP(tokamak_coordinates.Bxy() * phi) / tokamak_coordinates.Bxy() - eta * Jpar; } if (diamag) { - ddt(Psi) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Psi, bm_exb); // Equilibrium flow + ddt(Psi) -= bracket(tokamak_coordinates.Bxy() * phi0, Psi, bm_exb); // Equilibrium flow } // Hyper-resistivity @@ -1388,18 +1387,18 @@ class Elm_6f : public PhysicsModel { ddt(U) = 0.0; - ddt(U) = -SQ(tokamak_coordinates_factory.Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates_factory.Bxy(); // Grad j term + ddt(U) = -SQ(tokamak_coordinates.Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates.Bxy(); // Grad j term ddt(U) += 2.0 * Upara1 * b0xcv * Grad(P); // curvature term - ddt(U) += SQ(tokamak_coordinates_factory.Bxy()) * Grad_parP(Jpar); // b dot grad j + ddt(U) += SQ(tokamak_coordinates.Bxy()) * Grad_parP(Jpar); // b dot grad j if (diamag) { - ddt(U) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, U, bm_exb); // Equilibrium flow + ddt(U) -= bracket(tokamak_coordinates.Bxy() * phi0, U, bm_exb); // Equilibrium flow } if (nonlinear) { - ddt(U) -= bracket(tokamak_coordinates_factory.Bxy() * phi, U, bm_exb); // Advection + ddt(U) -= bracket(tokamak_coordinates.Bxy() * phi, U, bm_exb); // Advection } // parallel hyper-viscous diffusion for vector potential @@ -1451,18 +1450,18 @@ class Elm_6f : public PhysicsModel { ddt(Ni) = 0.0; - ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi, N0, bm_exb); + ddt(Ni) -= bracket(tokamak_coordinates.Bxy() * phi, N0, bm_exb); if (diamag) { - ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Ni, bm_exb); // Equilibrium flow + ddt(Ni) -= bracket(tokamak_coordinates.Bxy() * phi0, Ni, bm_exb); // Equilibrium flow } if (nonlinear) { - ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Ni, bm_exb); // Advection + ddt(Ni) -= bracket(tokamak_coordinates.Bxy() * phi, Ni, bm_exb); // Advection } if (compress0) { - ddt(Ni) -= N0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.Bxy()); + ddt(Ni) -= N0 * tokamak_coordinates.Bxy() * Grad_parP(Vipar / tokamak_coordinates.Bxy()); } // 4th order Parallel diffusion terms @@ -1481,18 +1480,18 @@ class Elm_6f : public PhysicsModel { ddt(Ti) = 0.0; - ddt(Ti) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Ti0, bm_exb); + ddt(Ti) -= bracket(tokamak_coordinates.Bxy() * phi, Ti0, bm_exb); if (diamag) { - ddt(Ti) -= bracket(phi0 * tokamak_coordinates_factory.Bxy(), Ti, bm_exb); // Equilibrium flow + ddt(Ti) -= bracket(phi0 * tokamak_coordinates.Bxy(), Ti, bm_exb); // Equilibrium flow } if (nonlinear) { - ddt(Ti) -= bracket(phi * tokamak_coordinates_factory.Bxy(), Ti, bm_exb); // Advection + ddt(Ti) -= bracket(phi * tokamak_coordinates.Bxy(), Ti, bm_exb); // Advection } if (compress0) { - ddt(Ti) -= 2.0 / 3.0 * Ti0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.Bxy()); + ddt(Ti) -= 2.0 / 3.0 * Ti0 * tokamak_coordinates.Bxy() * Grad_parP(Vipar / tokamak_coordinates.Bxy()); } if (diffusion_par > 0.0) { @@ -1517,18 +1516,18 @@ class Elm_6f : public PhysicsModel { ddt(Te) = 0.0; - ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Te0, bm_exb); + ddt(Te) -= bracket(tokamak_coordinates.Bxy() * phi, Te0, bm_exb); if (diamag) { - ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Te, bm_exb); // Equilibrium flow + ddt(Te) -= bracket(tokamak_coordinates.Bxy() * phi0, Te, bm_exb); // Equilibrium flow } if (nonlinear) { - ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Te, bm_exb); // Advection + ddt(Te) -= bracket(tokamak_coordinates.Bxy() * phi, Te, bm_exb); // Advection } if (compress0) { - ddt(Te) -= 2.0 / 3.0 * Te0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vepar / tokamak_coordinates_factory.Bxy()); + ddt(Te) -= 2.0 / 3.0 * Te0 * tokamak_coordinates.Bxy() * Grad_parP(Vepar / tokamak_coordinates.Bxy()); } if (diffusion_par > 0.0) { @@ -1551,14 +1550,14 @@ class Elm_6f : public PhysicsModel { ddt(Vipar) = 0.0; ddt(Vipar) -= Vipara * Grad_parP(P) / N0; - ddt(Vipar) += Vipara * bracket(Psi, P0, bm_mag) * tokamak_coordinates_factory.Bxy() / N0; + ddt(Vipar) += Vipara * bracket(Psi, P0, bm_mag) * tokamak_coordinates.Bxy() / N0; if (diamag) { - ddt(Vipar) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Vipar, bm_exb); + ddt(Vipar) -= bracket(tokamak_coordinates.Bxy() * phi0, Vipar, bm_exb); } if (nonlinear) { - ddt(Vipar) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Vipar, bm_exb); + ddt(Vipar) -= bracket(tokamak_coordinates.Bxy() * phi, Vipar, bm_exb); } // parallel hyper-viscous diffusion for vector potential diff --git a/examples/conducting-wall-mode/cwm.cxx b/examples/conducting-wall-mode/cwm.cxx index c869ea69ec..6d9a3d5d21 100644 --- a/examples/conducting-wall-mode/cwm.cxx +++ b/examples/conducting-wall-mode/cwm.cxx @@ -5,7 +5,7 @@ * Model version in the code created by M. Umansky and J. Myra. *******************************************************************************/ #include -#include +#include #include #include @@ -118,7 +118,7 @@ class CWM : public PhysicsModel { hthe0 / rho_s); } - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + auto tokamak_coordinates = TokamakCoordinates(*mesh); /************** NORMALISE QUANTITIES *****************/ @@ -129,8 +129,8 @@ class CWM : public PhysicsModel { Te0 /= Te_x; // Normalise geometry - coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, rho_s, - bmag / 1e4, ShearFactor); + coord = tokamak_coordinates.make_coordinates(noshear, rho_s, bmag / 1e4, + ShearFactor); // Set nu nu = nu_hat * Ni0 / pow(Te0, 1.5); @@ -141,10 +141,10 @@ class CWM : public PhysicsModel { // add evolving variables to the communication object SOLVE_FOR(rho, te); - Field2D Rxy = tokamak_coordinates_factory.Rxy(); - Field2D Bpxy = tokamak_coordinates_factory.Bpxy(); - Field2D Btxy = tokamak_coordinates_factory.Btxy(); - Field2D hthe = tokamak_coordinates_factory.hthe(); + Field2D Rxy = tokamak_coordinates.Rxy(); + Field2D Bpxy = tokamak_coordinates.Bpxy(); + Field2D Btxy = tokamak_coordinates.Btxy(); + Field2D hthe = tokamak_coordinates.hthe(); SAVE_ONCE(Rxy, Bpxy, Btxy, Zxy, hthe); SAVE_ONCE(nu_hat, hthe0); diff --git a/examples/constraints/alfven-wave/alfven.cxx b/examples/constraints/alfven-wave/alfven.cxx index 93c5d1abc2..40fb0eecc9 100644 --- a/examples/constraints/alfven-wave/alfven.cxx +++ b/examples/constraints/alfven-wave/alfven.cxx @@ -2,7 +2,7 @@ #include #include #include -#include +#include /// Fundamental constants const BoutReal PI = 3.14159265; @@ -170,9 +170,9 @@ class Alfven : public PhysicsModel { noshear = true; } - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + auto tokamak_coordinates = TokamakCoordinates(*mesh); const auto& coord = - tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lnorm, Bnorm); + tokamak_coordinates.make_coordinates(noshear, Lnorm, Bnorm); } }; diff --git a/examples/dalf3/dalf3.cxx b/examples/dalf3/dalf3.cxx index fdfc48db2a..d2217926f5 100644 --- a/examples/dalf3/dalf3.cxx +++ b/examples/dalf3/dalf3.cxx @@ -19,7 +19,7 @@ ****************************************************************/ #include -#include +#include #include #include @@ -81,7 +81,7 @@ class DALF3 : public PhysicsModel { std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) Field2D phi2D; // Axisymmetric potential, used when split_n0=true - TokamakCoordinatesFactory tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + TokamakCoordinates tokamak_coordinates = TokamakCoordinates(*mesh); protected: int init(bool UNUSED(restarting)) override { @@ -162,7 +162,7 @@ class DALF3 : public PhysicsModel { if (lowercase(ptstr) == "shifted") { // Dimits style, using local coordinate system - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; noshear = true; } @@ -223,7 +223,7 @@ class DALF3 : public PhysicsModel { // Metrics const auto& coord = - tokamak_coordinates_factory.make_tokamak_coordinates(noshear, rho_s, Bnorm); + tokamak_coordinates.make_coordinates(noshear, rho_s, Bnorm); 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 index 4ac9774ac2..78b2c0d2fb 100644 --- a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx +++ b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx @@ -46,7 +46,7 @@ #include #include #include -#include +#include #include // Defines BOUT_FOR_RAJA @@ -288,7 +288,7 @@ class ELMpb : public PhysicsModel { int damp_width; // Width of inner damped region BoutReal damp_t_const; // Timescale of damping - TokamakCoordinatesFactory tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + TokamakCoordinates tokamak_coordinates = TokamakCoordinates(*mesh); const BoutReal MU0 = 4.0e-7 * PI; const BoutReal Mi = 2.0 * 1.6726e-27; // Ion mass @@ -716,9 +716,10 @@ class ELMpb : public PhysicsModel { if (mesh->get(Lbar, "rmag") != 0) { // Typical length scale Lbar = 1.0; } - const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lbar, Bbar); + const auto& metric = + tokamak_coordinates.make_coordinates(noshear, Lbar, Bbar); - V0 = -tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Bpxy() * Dphi0 / tokamak_coordinates_factory.Bxy(); + V0 = -tokamak_coordinates.Rxy() * tokamak_coordinates.Bpxy() * Dphi0 / tokamak_coordinates.Bxy(); if (simple_rmp) { include_rmp = true; @@ -813,7 +814,7 @@ class ELMpb : public PhysicsModel { if (noshear) { if (include_curvature) { - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; } } @@ -823,7 +824,7 @@ class ELMpb : public PhysicsModel { if (not mesh->IncIntShear) { // Dimits style, using local coordinate system if (include_curvature) { - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; } } @@ -933,7 +934,7 @@ class ELMpb : public PhysicsModel { Field2D Te; Te = P0 / (2.0 * density * 1.602e-19); // Temperature in eV - J0 = -MU0 * Lbar * J0 / tokamak_coordinates_factory.Bxy(); + J0 = -MU0 * Lbar * J0 / tokamak_coordinates.Bxy(); P0 = 2.0 * MU0 * P0 / (Bbar * Bbar); V0 = V0 / Va; Dphi0 *= Tbar; @@ -1072,16 +1073,16 @@ class ELMpb : public PhysicsModel { B0vec.covariant = false; B0vec.x = 0.; - B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe(); + B0vec.y = tokamak_coordinates.Bpxy() / tokamak_coordinates.hthe(); B0vec.z = 0.; V0net.covariant = false; // presentation for net flow V0net.x = 0.; V0net.y = - tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Btxy() * tokamak_coordinates_factory.Bpxy() / (tokamak_coordinates_factory.hthe() * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy()) * Dphi0; + tokamak_coordinates.Rxy() * tokamak_coordinates.Btxy() * tokamak_coordinates.Bpxy() / (tokamak_coordinates.hthe() * tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy()) * Dphi0; V0net.z = -Dphi0; - U0 = B0vec * Curl(V0net) / tokamak_coordinates_factory + U0 = B0vec * Curl(V0net) / tokamak_coordinates .Bxy(); // get 0th vorticity for Kelvin-Holmholtz term /**************** SET EVOLVING VARIABLES *************/ @@ -1105,8 +1106,8 @@ class ELMpb : public PhysicsModel { SOLVE_FOR(Vpar); comms.add(Vpar); - beta = tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (0.5 + (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (g * P0))); - gradparB = Grad_par(tokamak_coordinates_factory.Bxy()) / tokamak_coordinates_factory.Bxy(); + beta = tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy() / (0.5 + (tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy() / (g * P0))); + gradparB = Grad_par(tokamak_coordinates.Bxy()) / tokamak_coordinates.Bxy(); output.write("Beta in range {:e} -> {:e}\n", min(beta), max(beta)); } else { @@ -1139,10 +1140,10 @@ class ELMpb : public PhysicsModel { // Diamagnetic phi0 if (diamag_phi0) { if (constn0) { - phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy(); + phi0 = -0.5 * dnorm * P0 / tokamak_coordinates.Bxy(); } else { // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift - phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy() / N0; + phi0 = -0.5 * dnorm * P0 / tokamak_coordinates.Bxy() / N0; } SAVE_ONCE(phi0); } else { @@ -1153,7 +1154,7 @@ class ELMpb : public PhysicsModel { // everything needed to recover physical units SAVE_ONCE(J0, P0); SAVE_ONCE(density, Lbar, Bbar, Tbar); - Field2D Bxy = tokamak_coordinates_factory.Bxy(); + Field2D Bxy = tokamak_coordinates.Bxy(); SAVE_ONCE(Va, Bxy); SAVE_ONCE(Dphi0, U0); SAVE_ONCE(V0); @@ -1213,7 +1214,7 @@ class ELMpb : public PhysicsModel { if (mesh->IncIntShear) { // BOUT-06 style, using d/dx = d/dpsi + I * d/dz - metric->setIntShiftTorsion(tokamak_coordinates_factory.ShearFactor()); + metric->setIntShiftTorsion(tokamak_coordinates.ShearFactor()); } @@ -1527,7 +1528,7 @@ class ELMpb : public PhysicsModel { auto P0_acc = Field2DAccessor<>(P0); auto J0_acc = Field2DAccessor<>(J0); auto phi0_acc = Field2DAccessor<>(phi0); - Field2D Bxy = tokamak_coordinates_factory.Bxy(); + Field2D Bxy = tokamak_coordinates.Bxy(); auto B0_acc = Field2DAccessor<>(Bxy); // Evolving fields @@ -1715,16 +1716,16 @@ class ELMpb : public PhysicsModel { // Vorticity equation // Grad j term - // ddt(U) = SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(Psi, J0, CELL_CENTRE); + // ddt(U) = SQ(tokamak_coordinates.Bxy()) * b0xGrad_dot_Grad(Psi, J0, CELL_CENTRE); // if (include_rmp) { - // ddt(U) += SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE); + // ddt(U) += SQ(tokamak_coordinates.Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE); // } ddt(U) += b0xcv * Grad(P); // curvature term // if (!nogradparj) { // Parallel current term - // ddt(U) -= SQ(tokamak_coordinates_factory.Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j + // ddt(U) -= SQ(tokamak_coordinates.Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j // } if (withflow && K_H_term) { // K_H_term diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index 3a3ac0fe9c..e321d3672a 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -14,7 +14,7 @@ #include #include #include -#include +#include #include #include @@ -230,7 +230,7 @@ class ELMpb : public PhysicsModel { int damp_width; // Width of inner damped region BoutReal damp_t_const; // Timescale of damping - TokamakCoordinatesFactory tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + TokamakCoordinates tokamak_coordinates = TokamakCoordinates(*mesh); const BoutReal MU0 = 4.0e-7 * PI; const BoutReal Mi = 2.0 * 1.6726e-27; // Ion mass @@ -645,7 +645,7 @@ class ELMpb : public PhysicsModel { Dphi0 *= -1; } - V0 = -tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Bpxy() * Dphi0 / tokamak_coordinates_factory.Bxy(); + V0 = -tokamak_coordinates.Rxy() * tokamak_coordinates.Bpxy() * Dphi0 / tokamak_coordinates.Bxy(); if (simple_rmp) { include_rmp = true; @@ -740,7 +740,7 @@ class ELMpb : public PhysicsModel { if (noshear) { if (include_curvature) { - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; } } @@ -750,7 +750,7 @@ class ELMpb : public PhysicsModel { if (not mesh->IncIntShear) { // Dimits style, using local coordinate system if (include_curvature) { - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; } } @@ -867,7 +867,7 @@ class ELMpb : public PhysicsModel { Field2D Te; Te = P0 / (2.0 * density * 1.602e-19); // Temperature in eV - J0 = -MU0 * Lbar * J0 / tokamak_coordinates_factory.Bxy(); + J0 = -MU0 * Lbar * J0 / tokamak_coordinates.Bxy(); P0 = 2.0 * MU0 * P0 / (Bbar * Bbar); V0 = V0 / Va; Dphi0 *= Tbar; @@ -876,7 +876,8 @@ class ELMpb : public PhysicsModel { b0xcv.y *= Lbar * Lbar; b0xcv.z *= Lbar * Lbar; - const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lbar, Bbar); + const auto& metric = + tokamak_coordinates.make_coordinates(noshear, Lbar, Bbar); ////////////////////////////////////////////////////////////// @@ -884,7 +885,7 @@ class ELMpb : public PhysicsModel { if (mesh->IncIntShear) { // BOUT-06 style, using d/dx = d/dpsi + I * d/dz - metric->setIntShiftTorsion(tokamak_coordinates_factory.ShearFactor()); + metric->setIntShiftTorsion(tokamak_coordinates.ShearFactor()); } if (constn0) { @@ -1016,16 +1017,16 @@ class ELMpb : public PhysicsModel { B0vec.covariant = false; B0vec.x = 0.; - B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe(); + B0vec.y = tokamak_coordinates.Bpxy() / tokamak_coordinates.hthe(); B0vec.z = 0.; V0net.covariant = false; // presentation for net flow V0net.x = 0.; V0net.y = - tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Btxy() * tokamak_coordinates_factory.Bpxy() / (tokamak_coordinates_factory.hthe() * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy()) * Dphi0; + tokamak_coordinates.Rxy() * tokamak_coordinates.Btxy() * tokamak_coordinates.Bpxy() / (tokamak_coordinates.hthe() * tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy()) * Dphi0; V0net.z = -Dphi0; - U0 = B0vec * Curl(V0net) / tokamak_coordinates_factory + U0 = B0vec * Curl(V0net) / tokamak_coordinates .Bxy(); // get 0th vorticity for Kelvin-Holmholtz term /**************** SET EVOLVING VARIABLES *************/ @@ -1049,8 +1050,8 @@ class ELMpb : public PhysicsModel { SOLVE_FOR(Vpar); comms.add(Vpar); - beta = tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (0.5 + (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (g * P0))); - gradparB = Grad_par(tokamak_coordinates_factory.Bxy()) / tokamak_coordinates_factory.Bxy(); + beta = tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy() / (0.5 + (tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy() / (g * P0))); + gradparB = Grad_par(tokamak_coordinates.Bxy()) / tokamak_coordinates.Bxy(); output.write("Beta in range {:e} -> {:e}\n", min(beta), max(beta)); } else { @@ -1083,10 +1084,10 @@ class ELMpb : public PhysicsModel { // Diamagnetic phi0 if (diamag_phi0) { if (constn0) { - phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy(); + phi0 = -0.5 * dnorm * P0 / tokamak_coordinates.Bxy(); } else { // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift - phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy() / N0; + phi0 = -0.5 * dnorm * P0 / tokamak_coordinates.Bxy() / N0; } SAVE_ONCE(phi0); } @@ -1095,7 +1096,7 @@ class ELMpb : public PhysicsModel { // everything needed to recover physical units SAVE_ONCE(J0, P0); SAVE_ONCE(density, Lbar, Bbar, Tbar); - Field2D tmp = tokamak_coordinates_factory.Bxy(); + Field2D tmp = tokamak_coordinates.Bxy(); SAVE_ONCE(Va, tmp); SAVE_ONCE(Dphi0, U0); SAVE_ONCE(V0); @@ -1163,7 +1164,7 @@ class ELMpb : public PhysicsModel { Field3D result = Grad_par(f, loc); if (nonlinear) { - const auto Bxy = tokamak_coordinates_factory.Bxy(); + const auto Bxy = tokamak_coordinates.Bxy(); result -= bracket(interp_to(Psi, loc), f, bm_mag) * Bxy; if (include_rmp) { @@ -1303,12 +1304,12 @@ class ELMpb : public PhysicsModel { } if (diamag) { - phi -= 0.5 * dnorm * P / tokamak_coordinates_factory.Bxy(); + phi -= 0.5 * dnorm * P / tokamak_coordinates.Bxy(); } } else { ubyn = U / N0; if (diamag) { - ubyn -= 0.5 * dnorm / (N0 * tokamak_coordinates_factory.Bxy()) * Delp2(P); + ubyn -= 0.5 * dnorm / (N0 * tokamak_coordinates.Bxy()) * Delp2(P); mesh->communicate(ubyn); } // Invert laplacian for phi @@ -1453,9 +1454,9 @@ class ELMpb : public PhysicsModel { if (evolve_jpar) { // Jpar - Field3D B0U = tokamak_coordinates_factory.Bxy() * U; + Field3D B0U = tokamak_coordinates.Bxy() * U; mesh->communicate(B0U); - ddt(Jpar) = -Grad_parP(B0U, loc) / tokamak_coordinates_factory.Bxy() + eta * Delp2(Jpar); + ddt(Jpar) = -Grad_parP(B0U, loc) / tokamak_coordinates.Bxy() + eta * Delp2(Jpar); if (relax_j_vac) { // Make ddt(Jpar) relax to zero. @@ -1464,16 +1465,16 @@ class ELMpb : public PhysicsModel { } } else { // Vector potential - ddt(Psi) = -Grad_parP(phi * tokamak_coordinates_factory.Bxy(), loc) / tokamak_coordinates_factory.Bxy() + eta * Jpar; + ddt(Psi) = -Grad_parP(phi * tokamak_coordinates.Bxy(), loc) / tokamak_coordinates.Bxy() + eta * Jpar; if (eHall) { // electron parallel pressure ddt(Psi) += 0.25 * delta_i - * (Grad_parP(tokamak_coordinates_factory.Bxy() * P, loc) / tokamak_coordinates_factory.Bxy() - + bracket(interp_to(P0, loc), Psi, bm_mag) * tokamak_coordinates_factory.Bxy()); + * (Grad_parP(tokamak_coordinates.Bxy() * P, loc) / tokamak_coordinates.Bxy() + + bracket(interp_to(P0, loc), Psi, bm_mag) * tokamak_coordinates.Bxy()); } if (diamag_phi0) { // Equilibrium flow - ddt(Psi) -= bracket(interp_to(phi0, loc), Psi, bm_exb) * tokamak_coordinates_factory.Bxy(); + ddt(Psi) -= bracket(interp_to(phi0, loc), Psi, bm_exb) * tokamak_coordinates.Bxy(); } if (withflow) { // net flow @@ -1481,7 +1482,7 @@ class ELMpb : public PhysicsModel { } if (diamag_grad_t) { // grad_par(T_e) correction - ddt(Psi) += 1.71 * dnorm * 0.5 * Grad_parP(P, loc) / tokamak_coordinates_factory.Bxy(); + ddt(Psi) += 1.71 * dnorm * 0.5 * Grad_parP(P, loc) / tokamak_coordinates.Bxy(); } if (hyperresist > 0.0) { // Hyper-resistivity @@ -1519,15 +1520,15 @@ class ELMpb : public PhysicsModel { Psi_loc = interp_to(Psi, CELL_CENTRE, "RGN_ALL"); Psi_loc.applyBoundary(); // Grad j term - ddt(U) = SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(Psi_loc, J0, CELL_CENTRE); + ddt(U) = SQ(tokamak_coordinates.Bxy()) * b0xGrad_dot_Grad(Psi_loc, J0, CELL_CENTRE); if (include_rmp) { - ddt(U) += SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE); + ddt(U) += SQ(tokamak_coordinates.Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE); } ddt(U) += b0xcv * Grad(P); // curvature term if (!nogradparj) { // Parallel current term - ddt(U) -= SQ(tokamak_coordinates_factory.Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j + ddt(U) -= SQ(tokamak_coordinates.Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j } if (withflow && K_H_term) { // K_H_term @@ -1543,7 +1544,7 @@ class ELMpb : public PhysicsModel { } if (nonlinear) { // Advection - ddt(U) -= bracket(phi, U, bm_exb) * tokamak_coordinates_factory.Bxy(); + ddt(U) -= bracket(phi, U, bm_exb) * tokamak_coordinates.Bxy(); } // Viscosity terms @@ -1587,11 +1588,11 @@ class ELMpb : public PhysicsModel { Pi = 0.5 * P; Pi0 = 0.5 * P0; - Dperp2Phi0 = Field3D(Delp2(tokamak_coordinates_factory.Bxy() * phi0)); + Dperp2Phi0 = Field3D(Delp2(tokamak_coordinates.Bxy() * phi0)); Dperp2Phi0.applyBoundary(); mesh->communicate(Dperp2Phi0); - Dperp2Phi = Delp2(tokamak_coordinates_factory.Bxy() * phi); + Dperp2Phi = Delp2(tokamak_coordinates.Bxy() * phi); Dperp2Phi.applyBoundary(); mesh->communicate(Dperp2Phi); @@ -1603,35 +1604,35 @@ class ELMpb : public PhysicsModel { Dperp2Pi.applyBoundary(); mesh->communicate(Dperp2Pi); - bracketPhi0P = bracket(tokamak_coordinates_factory.Bxy() * phi0, Pi, bm_exb); + bracketPhi0P = bracket(tokamak_coordinates.Bxy() * phi0, Pi, bm_exb); bracketPhi0P.applyBoundary(); mesh->communicate(bracketPhi0P); - bracketPhiP0 = bracket(tokamak_coordinates_factory.Bxy() * phi, Pi0, bm_exb); + bracketPhiP0 = bracket(tokamak_coordinates.Bxy() * phi, Pi0, bm_exb); bracketPhiP0.applyBoundary(); mesh->communicate(bracketPhiP0); - ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi0, bm_exb) / tokamak_coordinates_factory.Bxy(); - ddt(U) -= 0.5 * Upara2 * bracket(Pi0, Dperp2Phi, bm_exb) / tokamak_coordinates_factory.Bxy(); - Field3D B0phi = tokamak_coordinates_factory.Bxy() * phi; + ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi0, bm_exb) / tokamak_coordinates.Bxy(); + ddt(U) -= 0.5 * Upara2 * bracket(Pi0, Dperp2Phi, bm_exb) / tokamak_coordinates.Bxy(); + Field3D B0phi = tokamak_coordinates.Bxy() * phi; mesh->communicate(B0phi); - Field3D B0phi0 = tokamak_coordinates_factory.Bxy() * phi0; + Field3D B0phi0 = tokamak_coordinates.Bxy() * phi0; mesh->communicate(B0phi0); - ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / tokamak_coordinates_factory.Bxy(); - ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / tokamak_coordinates_factory.Bxy(); - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / tokamak_coordinates_factory.Bxy(); - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / tokamak_coordinates_factory.Bxy(); + ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / tokamak_coordinates.Bxy(); + ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / tokamak_coordinates.Bxy(); + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / tokamak_coordinates.Bxy(); + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / tokamak_coordinates.Bxy(); if (nonlinear) { - Field3D B0phi = tokamak_coordinates_factory.Bxy() * phi; + Field3D B0phi = tokamak_coordinates.Bxy() * phi; mesh->communicate(B0phi); bracketPhiP = bracket(B0phi, Pi, bm_exb); bracketPhiP.applyBoundary(); mesh->communicate(bracketPhiP); - ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / tokamak_coordinates_factory.Bxy(); - ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / tokamak_coordinates_factory.Bxy(); - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / tokamak_coordinates_factory.Bxy(); + ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / tokamak_coordinates.Bxy(); + ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / tokamak_coordinates.Bxy(); + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / tokamak_coordinates.Bxy(); } } @@ -1661,7 +1662,7 @@ class ELMpb : public PhysicsModel { } if (nonlinear) { // Advection - ddt(P) -= bracket(phi, P, bm_exb) * tokamak_coordinates_factory.Bxy(); + ddt(P) -= bracket(phi, P, bm_exb) * tokamak_coordinates.Bxy(); } } @@ -1706,7 +1707,7 @@ class ELMpb : public PhysicsModel { ddt(Vpar) = -0.5 * (Grad_par(P, loc) + Grad_par(P0, loc)); if (nonlinear) { - ddt(Vpar) -= bracket(interp_to(phi, loc), Vpar, bm_exb) * tokamak_coordinates_factory.Bxy(); // Advection + ddt(Vpar) -= bracket(interp_to(phi, loc), Vpar, bm_exb) * tokamak_coordinates.Bxy(); // Advection } } @@ -1797,7 +1798,7 @@ class ELMpb : public PhysicsModel { mesh->communicate(Jrhs, ddt(P)); Field3D U1 = ddt(U); - U1 += (gamma * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P); + U1 += (gamma * tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P); // Second matrix, solving Alfven wave dynamics static std::unique_ptr invU{nullptr}; @@ -1806,7 +1807,7 @@ class ELMpb : public PhysicsModel { } invU->setCoefA(1.); - invU->setCoefB(-SQ(gamma) * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy()); + invU->setCoefB(-SQ(gamma) * tokamak_coordinates.Bxy() * tokamak_coordinates.Bxy()); ddt(U) = invU->solve(U1); ddt(U).applyBoundary(); @@ -1814,9 +1815,9 @@ class ELMpb : public PhysicsModel { Field3D phi3 = phiSolver->solve(ddt(U)); mesh->communicate(phi3); phi3.applyBoundary("neumann"); - Field3D B0phi3 = tokamak_coordinates_factory.Bxy() * phi3; + Field3D B0phi3 = tokamak_coordinates.Bxy() * phi3; mesh->communicate(B0phi3); - ddt(Psi) = ddt(Psi) - gamma * Grad_par(B0phi3, loc) / tokamak_coordinates_factory.Bxy(); + ddt(Psi) = ddt(Psi) - gamma * Grad_par(B0phi3, loc) / tokamak_coordinates.Bxy(); ddt(Psi).applyBoundary(); return 0; @@ -1848,14 +1849,14 @@ class ELMpb : public PhysicsModel { Field3D JP = -b0xGrad_dot_Grad(phi, P0); JP.setBoundary("P"); JP.applyBoundary(); - Field3D B0phi = tokamak_coordinates_factory.Bxy() * phi; + Field3D B0phi = tokamak_coordinates.Bxy() * phi; mesh->communicate(B0phi); - Field3D JPsi = -Grad_par(B0phi, loc) / tokamak_coordinates_factory.Bxy(); + Field3D JPsi = -Grad_par(B0phi, loc) / tokamak_coordinates.Bxy(); JPsi.setBoundary("Psi"); JPsi.applyBoundary(); - Field3D JU = b0xcv * Grad(ddt(P)) - SQ(tokamak_coordinates_factory.Bxy()) * Grad_par(Jpar, CELL_CENTRE) - + SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE); + Field3D JU = b0xcv * Grad(ddt(P)) - SQ(tokamak_coordinates.Bxy()) * Grad_par(Jpar, CELL_CENTRE) + + SQ(tokamak_coordinates.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 index be92ddea0a..5f401528f9 100644 --- a/examples/gyro-gem/gem.cxx +++ b/examples/gyro-gem/gem.cxx @@ -10,7 +10,7 @@ ****************************************************************/ #include -#include +#include #include #include @@ -242,17 +242,17 @@ class GEM : public PhysicsModel { Tbar = options["Tbar"].withDefault(Tbar); // Override in options file SAVE_ONCE(Tbar); // Timescale in seconds - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + auto tokamak_coordinates = TokamakCoordinates(*mesh); if (mesh->get(Bbar, "Bbar")) { if (mesh->get(Bbar, "bmag")) { - Bbar = max(tokamak_coordinates_factory.Bxy(), true); + Bbar = max(tokamak_coordinates.Bxy(), true); } } Bbar = options["Bbar"].withDefault(Bbar); // Override in options file SAVE_ONCE(Bbar); - coord = tokamak_coordinates_factory.make_tokamak_coordinates(true, Lbar, Bbar); + coord = tokamak_coordinates.make_coordinates(true, Lbar, Bbar); beta_e = 4.e-7 * PI * max(p_e, true) / (Bbar * Bbar); SAVE_ONCE(beta_e); @@ -349,7 +349,7 @@ class GEM : public PhysicsModel { B0vec.covariant = false; B0vec.x = 0.; - B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe(); + B0vec.y = tokamak_coordinates.Bpxy() / tokamak_coordinates.hthe(); B0vec.z = 0.; // Precompute this for use in RHS diff --git a/examples/laplacexy/alfven-wave/alfven.cxx b/examples/laplacexy/alfven-wave/alfven.cxx index fbe1de3e9d..375c213132 100644 --- a/examples/laplacexy/alfven-wave/alfven.cxx +++ b/examples/laplacexy/alfven-wave/alfven.cxx @@ -4,7 +4,7 @@ #include #include #include -#include +#include /// Fundamental constants const BoutReal PI = 3.14159265; @@ -172,9 +172,8 @@ class Alfven : public PhysicsModel { void LoadMetric(BoutReal Lnorm, BoutReal Bnorm) { - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); - const auto& coord = - tokamak_coordinates_factory.make_tokamak_coordinates(true, Lnorm, Bnorm); + auto tokamak_coordinates = TokamakCoordinates(*mesh); + const auto& coord = tokamak_coordinates.make_coordinates(true, Lnorm, Bnorm); } }; diff --git a/examples/laplacexy/laplace_perp/test.cxx b/examples/laplacexy/laplace_perp/test.cxx index 7d3c8972c4..485e321eb6 100644 --- a/examples/laplacexy/laplace_perp/test.cxx +++ b/examples/laplacexy/laplace_perp/test.cxx @@ -1,4 +1,4 @@ -#include +#include #include #include @@ -14,8 +14,8 @@ int main(int argc, char** argv) { bool calc_metric; 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); + auto tokamak_coordinates = TokamakCoordinates(*mesh); + const auto& coord = tokamak_coordinates.make_coordinates(true); } /////////////////////////////////////// diff --git a/examples/wave-slab/wave_slab.cxx b/examples/wave-slab/wave_slab.cxx index 6ac11568e3..efc3977084 100644 --- a/examples/wave-slab/wave_slab.cxx +++ b/examples/wave-slab/wave_slab.cxx @@ -11,14 +11,14 @@ */ #include -#include +#include class WaveTest : public PhysicsModel { public: int init(bool UNUSED(restarting)) { - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); - const auto& coords = tokamak_coordinates_factory.make_tokamak_coordinates(true); + auto tokamak_coordinates = TokamakCoordinates(*mesh); + const auto& coords = tokamak_coordinates.make_coordinates(true); solver->add(f, "f"); solver->add(g, "g"); diff --git a/include/bout/tokamak_coordinates_factory.hxx b/include/bout/tokamak_coordinates.hxx similarity index 88% rename from include/bout/tokamak_coordinates_factory.hxx rename to include/bout/tokamak_coordinates.hxx index 5690d4fdab..7a9833c025 100644 --- a/include/bout/tokamak_coordinates_factory.hxx +++ b/include/bout/tokamak_coordinates.hxx @@ -1,6 +1,6 @@ -#ifndef BOUT_TOKAMAK_COORDINATES_FACTORY_HXX -#define BOUT_TOKAMAK_COORDINATES_FACTORY_HXX +#ifndef BOUT_TOKAMAK_COORDINATES_HXX +#define BOUT_TOKAMAK_COORDINATES_HXX #include "bout.hxx" #include "bout/bout_types.hxx" @@ -10,7 +10,7 @@ #include "bout/utils.hxx" #include "bout/vector2d.hxx" -class TokamakCoordinatesFactory { +class TokamakCoordinates { private: Mesh& mesh_m; @@ -34,7 +34,7 @@ private: } public: - TokamakCoordinatesFactory(Mesh& mesh) : mesh_m(mesh) { + TokamakCoordinates(Mesh& mesh) : mesh_m(mesh) { mesh.get(Rxy_m, "Rxy"); // mesh->get(Zxy, "Zxy"); @@ -53,7 +53,7 @@ public: return 1.0; } - Coordinates* make_tokamak_coordinates(const bool noshear, BoutReal Lbar, BoutReal Bbar, + Coordinates* make_coordinates(const bool noshear, BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor = 1.0) { normalise(Lbar, Bbar, ShearFactor); @@ -98,4 +98,4 @@ public: const FieldMetric& ShearFactor() const { return ShearFactor_m; } }; -#endif //BOUT_TOKAMAK_COORDINATES_FACTORY_HXX +#endif //BOUT_TOKAMAK_COORDINATES_HXX diff --git a/tests/MMS/GBS/gbs.cxx b/tests/MMS/GBS/gbs.cxx index 34fa681945..eb14696cc4 100644 --- a/tests/MMS/GBS/gbs.cxx +++ b/tests/MMS/GBS/gbs.cxx @@ -315,8 +315,8 @@ int GBS::init(bool restarting) { void GBS::LoadMetric(BoutReal Lnorm, BoutReal Bnorm) { - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); - coords = tokamak_coordinates_factory.make_tokamak_coordinates(true, Lnorm, Bnorm); + auto tokamak_coordinates = TokamakCoordinates(*mesh); + coords = tokamak_coordinates.make_coordinates(true, Lnorm, Bnorm); // just define a macro for V_E dot Grad #define vE_Grad(f, p) (bracket(p, f, bm_exb)) diff --git a/tests/MMS/elm-pb/elm_pb.cxx b/tests/MMS/elm-pb/elm_pb.cxx index 2342b6c69e..ae97d7d347 100644 --- a/tests/MMS/elm-pb/elm_pb.cxx +++ b/tests/MMS/elm-pb/elm_pb.cxx @@ -122,7 +122,7 @@ class ELMpb : public PhysicsModel { // Communication objects FieldGroup comms; - TokamakCoordinatesFactory tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + TokamakCoordinates tokamak_coordinates = TokamakCoordinates(*mesh); // Parallel gradient along perturbed field-line @@ -297,7 +297,7 @@ class ELMpb : public PhysicsModel { if (not mesh->IncIntShear) { // Dimits style, using local coordinate system if (include_curvature) { - b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.get_ShearFactor() * b0xcv.x; } } } @@ -390,7 +390,7 @@ class ELMpb : public PhysicsModel { noshear = true; } - coords = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lbar, Bbar); + coords = tokamak_coordinates.make_coordinates(noshear, Lbar, Bbar); ////////////////////////////////////////////////////////////// // SHIFTED RADIAL COORDINATES @@ -400,7 +400,7 @@ class ELMpb : public PhysicsModel { if (ShiftXderivs) { if (mesh->IncIntShear) { // BOUT-06 style, using d/dx = d/dpsi + I * d/dz - coord->setIntShiftTorsion(tokamak_coordinates_factory.get_ShearFactor()); + coord->setIntShiftTorsion(tokamak_coordinates.get_ShearFactor()); } } diff --git a/tests/MMS/tokamak/tokamak.cxx b/tests/MMS/tokamak/tokamak.cxx index 2559fbcbd1..acffdcf882 100644 --- a/tests/MMS/tokamak/tokamak.cxx +++ b/tests/MMS/tokamak/tokamak.cxx @@ -45,9 +45,9 @@ class TokamakMMS : public PhysicsModel { } void LoadMetric(BoutReal Lnorm, BoutReal Bnorm) { - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + auto tokamak_coordinates = TokamakCoordinates(*mesh); auto coords = - tokamak_coordinates_factory.make_tokamak_coordinates(true, Lnorm, Bnorm); + tokamak_coordinates.make_coordinates(true, Lnorm, Bnorm); } private: diff --git a/tests/integrated/test-drift-instability/2fluid.cxx b/tests/integrated/test-drift-instability/2fluid.cxx index 23dbb7e837..f8c38cebe0 100644 --- a/tests/integrated/test-drift-instability/2fluid.cxx +++ b/tests/integrated/test-drift-instability/2fluid.cxx @@ -4,7 +4,7 @@ *******************************************************************************/ #include -#include +#include #include #include @@ -59,7 +59,7 @@ class TwoFluid : public PhysicsModel { FieldGroup comms; // Group of variables for communications - TokamakCoordinatesFactory tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + TokamakCoordinates tokamak_coordinates = TokamakCoordinates(*mesh); Coordinates* coord; // Coordinate system @@ -134,7 +134,7 @@ class TwoFluid : public PhysicsModel { const bool ShiftXderivs = (*globalOptions)["ShiftXderivs"].withDefault(false); if (ShiftXderivs) { ShearFactor = 0.0; // I disappears from metric - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; noshear = true; } @@ -193,8 +193,7 @@ class TwoFluid : public PhysicsModel { pei0 = (Ti0 + Te0) * Ni0; pe0 = Te0 * Ni0; - coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, rho_s, - bmag / 1e4, ShearFactor); + coord = tokamak_coordinates.make_coordinates(noshear, rho_s, bmag / 1e4, ShearFactor); /**************** SET EVOLVING VARIABLES *************/ diff --git a/tests/integrated/test-interchange-instability/2fluid.cxx b/tests/integrated/test-interchange-instability/2fluid.cxx index 843797e365..3845e7d050 100644 --- a/tests/integrated/test-interchange-instability/2fluid.cxx +++ b/tests/integrated/test-interchange-instability/2fluid.cxx @@ -3,7 +3,7 @@ * Same as Maxim's version of BOUT - simplified 2-fluid for benchmarking *******************************************************************************/ -#include +#include #include #include @@ -34,7 +34,7 @@ class Interchange : public PhysicsModel { Coordinates* coord; - TokamakCoordinatesFactory tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); + TokamakCoordinates tokamak_coordinates = TokamakCoordinates(*mesh); protected: int init(bool UNUSED(restarting)) override { @@ -100,11 +100,11 @@ class Interchange : public PhysicsModel { hthe0 / rho_s); } - coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, rho_s, - bmag / 1e4, ShearFactor); + coord = tokamak_coordinates.make_coordinates(noshear, rho_s, bmag / 1e4, + ShearFactor); if (ShiftXderivs) { - b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x; + b0xcv.z += tokamak_coordinates.ShearFactor() * b0xcv.x; } /************** NORMALISE QUANTITIES *****************/ diff --git a/tests/integrated/test-laplacexy/loadmetric.cxx b/tests/integrated/test-laplacexy/loadmetric.cxx index 437c42b3b8..a25c92baae 100644 --- a/tests/integrated/test-laplacexy/loadmetric.cxx +++ b/tests/integrated/test-laplacexy/loadmetric.cxx @@ -17,6 +17,6 @@ void LoadMetric(BoutReal Lnorm, BoutReal Bnorm) { noshear = true; } - auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh); - auto coords = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lnorm, Bnorm); + auto tokamak_coordinates = TokamakCoordinates(*mesh); + auto coords = tokamak_coordinates.make_coordinates(noshear, Lnorm, Bnorm); }