diff --git a/external/BOUT-dev b/external/BOUT-dev index a511a6632..ae5e990c2 160000 --- a/external/BOUT-dev +++ b/external/BOUT-dev @@ -1 +1 @@ -Subproject commit a511a6632640ebfe2df2f1fd483dcd72c6702319 +Subproject commit ae5e990c2fc3fef4bad202d45e88cdd436acd245 diff --git a/hermes-3.cxx b/hermes-3.cxx index eb0dd72ff..d3b71e2c2 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -121,8 +121,8 @@ class DecayLengthBoundary : public BoundaryOp { // Get cell radial length Coordinates* coord = mesh->getCoordinates(); - auto dx = coord->dx; - auto g11 = coord->g11; + auto dx = coord->dx(); + auto g11 = coord->g11(); Coordinates::FieldMetric dr = dx / sqrt(g11); // cell radial length. dr = dx/(Bpol * R) and g11 = (Bpol*R)**2 @@ -254,28 +254,7 @@ int Hermes::init(bool restarting) { .doc("Normalise input metric tensor? (assumes input is in SI units)") .withDefault(true)) { Coordinates* coord = mesh->getCoordinates(); - // To use non-orthogonal metric - // Normalise - coord->dx /= rho_s0 * rho_s0 * Bnorm; - coord->Bxy /= Bnorm; - // Metric is in grid file - just need to normalise - coord->g11 /= SQ(Bnorm * rho_s0); - coord->g22 *= SQ(rho_s0); - coord->g33 *= SQ(rho_s0); - coord->g12 /= Bnorm; - coord->g13 /= Bnorm; - coord->g23 *= SQ(rho_s0); - - coord->J *= Bnorm / rho_s0; - - coord->g_11 *= SQ(Bnorm * rho_s0); - coord->g_22 /= SQ(rho_s0); - coord->g_33 /= SQ(rho_s0); - coord->g_12 *= Bnorm; - coord->g_13 *= Bnorm; - coord->g_23 /= SQ(rho_s0); - - coord->geometry(); // Calculate other metrics + coord->normaliseMetric(rho_s0, Bnorm); } } diff --git a/include/braginskii_conduction.hxx b/include/braginskii_conduction.hxx index 6356a9c8b..f40f96096 100644 --- a/include/braginskii_conduction.hxx +++ b/include/braginskii_conduction.hxx @@ -79,7 +79,7 @@ private: /// void transform_impl(GuardedOptions& state) override; - std::string conduction_method{"orginal"}; + bout::ConductionMethod conduction_method{bout::ConductionMethod::Original}; }; namespace { diff --git a/include/div_ops.hxx b/include/div_ops.hxx index d02f054c3..31fcb462f 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -43,7 +43,7 @@ /*! * Diffusion in index space * - * Similar to using Div_par_diffusion(SQ(mesh->dy)*mesh->g_22, f) + * Similar to using Div_par_diffusion(SQ(mesh->dy())*mesh->g_22(), f) * * @param[in] The field to be differentiated * @param[in] bndry_flux Are fluxes through the boundary calculated? diff --git a/include/integrate.hxx b/include/integrate.hxx index 09a10114d..af9a88c70 100644 --- a/include/integrate.hxx +++ b/include/integrate.hxx @@ -75,7 +75,7 @@ auto cellAverage(Function func, const RegionType& region) { result.allocate(); // Get the coordinate Jacobian - auto J = result.getCoordinates()->J; + auto J = result.getCoordinates()->J(); BOUT_FOR(i, region) { // Offset indices auto yp = i.yp(); diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 6b164bc0a..feb89425b 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -82,7 +82,7 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? - std::string conduction_method{"orginal"}; + bout::ConductionMethod conduction_method{bout::ConductionMethod::Original}; // Flow diagnostics Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; diff --git a/include/polarisation_drift.hxx b/include/polarisation_drift.hxx index a302ab22a..c59ec7d5b 100644 --- a/include/polarisation_drift.hxx +++ b/include/polarisation_drift.hxx @@ -69,7 +69,7 @@ struct PolarisationDrift : public NamedComponent { private: std::unique_ptr phiSolver; // Laplacian solver in X-Z - Coordinates::FieldMetric Bsq; // Cached SQ(coord->Bxy) + Coordinates::FieldMetric Bsq; // Cached SQ(coord->Bxy()) // Diagnostic outputs bool diagnose; ///< Save diagnostic outputs? diff --git a/include/rate_helper.hxx b/include/rate_helper.hxx index ffc985eb6..88b1054ea 100644 --- a/include/rate_helper.hxx +++ b/include/rate_helper.hxx @@ -156,7 +156,7 @@ struct RateHelper { // Populate cell_props differently according to the rate function type std::visit( [this, &cell_props, &do_averaging, &result](auto&& rate_calc_func) { - auto J = result.rate.getCoordinates()->J; + auto J = result.rate.getCoordinates()->J(); BOUT_FOR(i, region) { // Get densities for this cell collect_densities(i, cell_props, do_averaging); diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index ed7717497..ca9877af4 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -116,7 +116,7 @@ private: BoutReal phi_boundary_last_update; ///< Time when last updated bool phi_core_averagey; ///< Average phi core boundary in Y? - Coordinates::FieldMetric Bsq; ///< SQ(coord->Bxy) + Coordinates::FieldMetric Bsq; ///< SQ(coord->Bxy()) VectorMetric Curlb_B; ///< Curvature vector Curl(b/B) BoutReal hyper_z; ///< Hyper-viscosity in Z Coordinates::FieldMetric viscosity; ///< Perpendicular Kinematic viscosity diff --git a/include/scale_timederivs.hxx b/include/scale_timederivs.hxx index 5889d4ace..a2a1352f5 100644 --- a/include/scale_timederivs.hxx +++ b/include/scale_timederivs.hxx @@ -40,7 +40,7 @@ private: void transform_impl(GuardedOptions& state) override { auto* coord = bout::globals::mesh->getCoordinates(); - auto dl2 = coord->g_22 * SQ(coord->dy); + auto dl2 = coord->g_22() * SQ(coord->dy()); // Scale by parallel heat conduction CFL timescale auto Te = get(state["species"]["e"]["temperature"]); diff --git a/include/vorticity.hxx b/include/vorticity.hxx index d461bae58..5ef018b3d 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -127,7 +127,7 @@ private: bool split_n0; // Split phi into n=0 and n!=0 components std::unique_ptr laplacexy; // Laplacian solver in X-Y (n=0) - Coordinates::FieldMetric Bsq; // SQ(coord->Bxy) + Coordinates::FieldMetric Bsq; // SQ(coord->Bxy()) VectorMetric Curlb_B; // Curvature vector Curl(b/B) BoutReal hyper_z; ///< Hyper-viscosity in Z Coordinates::FieldMetric viscosity; ///< Kinematic viscosity diff --git a/src/braginskii_conduction.cxx b/src/braginskii_conduction.cxx index c29bd45b6..58d3a4e66 100644 --- a/src/braginskii_conduction.cxx +++ b/src/braginskii_conduction.cxx @@ -114,8 +114,7 @@ BraginskiiConduction::BraginskiiConduction(const std::string& name, Options& all } substitutePermissions("sp", species); - conduction_method = - alloptions["conduction_method"].withDefault(conduction_method); + conduction_method = alloptions["conduction_method"].withDefault(conduction_method); } void BraginskiiConduction::transform_impl(GuardedOptions& state) { diff --git a/src/braginskii_electron_viscosity.cxx b/src/braginskii_electron_viscosity.cxx index 36cc64c14..362700697 100644 --- a/src/braginskii_electron_viscosity.cxx +++ b/src/braginskii_electron_viscosity.cxx @@ -50,7 +50,7 @@ void BraginskiiElectronViscosity::transform_impl(GuardedOptions& state) { const Field3D V = get(species["velocity"]); Coordinates* coord = P.getCoordinates(); - const Field3D Bxy = coord->Bxy; + const Field3D Bxy = coord->Bxy(); const Field3D sqrtB = sqrt(Bxy); // Parallel electron viscosity diff --git a/src/braginskii_ion_viscosity.cxx b/src/braginskii_ion_viscosity.cxx index b104a3184..6884f72c7 100644 --- a/src/braginskii_ion_viscosity.cxx +++ b/src/braginskii_ion_viscosity.cxx @@ -120,7 +120,7 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name, Curlb_B.y *= SQ(Lnorm); Curlb_B.z *= SQ(Lnorm); - Curlb_B *= 2. / coord->Bxy; + Curlb_B *= 2. / coord->Bxy(); } if (bounce_frequency) { const Options& units = alloptions["units"]; @@ -146,7 +146,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { GuardedOptions allspecies = state["species"]; auto coord = mesh->getCoordinates(); - const auto Bxy = coord->Bxy; + const auto Bxy = coord->Bxy(); const auto sqrtB = sqrt(Bxy); const auto Grad_par_logB = Grad_par(log(Bxy)); diff --git a/src/classical_diffusion.cxx b/src/classical_diffusion.cxx index d23cf72bb..ea2713ff6 100644 --- a/src/classical_diffusion.cxx +++ b/src/classical_diffusion.cxx @@ -12,7 +12,7 @@ ClassicalDiffusion::ClassicalDiffusion(std::string name, Options& alloptions, So readWrite("species:{all_species}:{output}")}) { Options& options = alloptions[name]; - Bsq = SQ(bout::globals::mesh->getCoordinates()->Bxy); + Bsq = SQ(bout::globals::mesh->getCoordinates()->Bxy()); diagnose = options["diagnose"].doc("Output additional diagnostics?").withDefault(false); diff --git a/src/diamagnetic_drift.cxx b/src/diamagnetic_drift.cxx index 041e4a6f1..e8f242f29 100644 --- a/src/diamagnetic_drift.cxx +++ b/src/diamagnetic_drift.cxx @@ -47,7 +47,7 @@ DiamagneticDrift::DiamagneticDrift(std::string name, Options& alloptions, Curlb_B.y *= SQ(Lnorm); Curlb_B.z *= SQ(Lnorm); - Curlb_B *= 2. / mesh->getCoordinates()->Bxy; + Curlb_B *= 2. / mesh->getCoordinates()->Bxy(); // Set drift to zero through sheath boundaries. // Flux through those cell faces should be set by sheath. diff --git a/src/div_ops.cxx b/src/div_ops.cxx index 2160e6606..c7809f5ba 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -310,16 +310,16 @@ Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flu for (int k = 0; k < mesh->LocalNz; k++) { // Average dfdy to right X boundary - BoutReal f_R = - 0.5 - * ((coord->g11(i + 1, j, k) * coord->g23(i + 1, j, k) - / SQ(coord->Bxy(i + 1, j, k))) - * dfdy(i + 1, j, k) - + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) - * dfdy(i, j, k)); + BoutReal f_R = 0.5 + * ((coord->g11(i + 1, j, k) * coord->g23(i + 1, j, k) + / SQ(coord->Bxy()(i + 1, j, k))) + * dfdy(i + 1, j, k) + + (coord->g11()(i, j, k) * coord->g23()(i, j, k) + / SQ(coord->Bxy()(i, j, k))) + * dfdy(i, j, k)); // Advection velocity across cell face - BoutReal Vx = 0.5 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) * f_R; + BoutReal Vx = 0.5 * (coord->J()(i + 1, j, k) + coord->J()(i, j, k)) * f_R; // Fromm method BoutReal flux = Vx; @@ -381,15 +381,15 @@ Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flu // Y flow // Average dfdx to upper Y boundary - BoutReal f_U = - 0.5 - * ((coord->g11(i, j + 1, k) * coord->g23(i, j + 1, k) - / SQ(coord->Bxy(i, j + 1, k))) - * dfdx(i, j + 1, k) - + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) - * dfdx(i, j, k)); + BoutReal f_U = 0.5 + * ((coord->g11(i, j + 1, k) * coord->g23(i, j + 1, k) + / SQ(coord->Bxy()(i, j + 1, k))) + * dfdx(i, j + 1, k) + + (coord->g11()(i, j, k) * coord->g23()(i, j, k) + / SQ(coord->Bxy()(i, j, k))) + * dfdx(i, j, k)); - BoutReal Vy = -0.5 * (coord->J(i, j + 1, k) + coord->J(i, j, k)) * f_U; + BoutReal Vy = -0.5 * (coord->J()(i, j + 1, k) + coord->J()(i, j, k)) * f_U; if (mesh->firstY(i) && !mesh->periodicY(i) && (j == mesh->ystart - 1)) { // Lower y boundary. Allow flows out of the domain only @@ -745,12 +745,12 @@ Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, Field3D& flow_ Field3D fc = f; Field3D ac = a; - Field3D g23c = coord->g23; - Field3D g_23c = coord->g_23; - Field3D Jc = coord->J; - Field3D dyc = coord->dy; - Field3D dzc = coord->dz; - Field3D Bxyc = coord->Bxy; + Field3D g23c = coord->g23(); + Field3D g_23c = coord->g_23(); + Field3D Jc = coord->J(); + Field3D dyc = coord->dy(); + Field3D dzc = coord->dz(); + Field3D Bxyc = coord->Bxy(); // Result of the Y and Z fluxes Field3D yzresult(mesh); @@ -779,39 +779,39 @@ Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, Field3D& flow_ // 3D Metric, need yup/ydown fields. // Requires previous communication of metrics // -- should insert communication here? - if (!coord->g23.hasParallelSlices() || !coord->g_23.hasParallelSlices() - || !coord->dy.hasParallelSlices() || !coord->dz.hasParallelSlices() - || !coord->Bxy.hasParallelSlices() || !coord->J.hasParallelSlices()) { + if (!coord->g23().hasParallelSlices() || !coord->g_23().hasParallelSlices() + || !coord->dy().hasParallelSlices() || !coord->dz().hasParallelSlices() + || !coord->Bxy().hasParallelSlices() || !coord->J().hasParallelSlices()) { throw BoutException("metrics have no yup/down: Maybe communicate in init?"); } - g23up = coord->g23.yup(); - g23down = coord->g23.ydown(); + g23up = coord->g23().yup(); + g23down = coord->g23().ydown(); - g_23up = coord->g_23.yup(); - g_23down = coord->g_23.ydown(); + g_23up = coord->g_23().yup(); + g_23down = coord->g_23().ydown(); - Jup = coord->J.yup(); - Jdown = coord->J.ydown(); + Jup = coord->J().yup(); + Jdown = coord->J().ydown(); - dyup = coord->dy.yup(); - dydown = coord->dy.ydown(); + dyup = coord->dy().yup(); + dydown = coord->dy().ydown(); - dzup = coord->dz.yup(); - dzdown = coord->dz.ydown(); + dzup = coord->dz().yup(); + dzdown = coord->dz().ydown(); - Bxyup = coord->Bxy.yup(); - Bxydown = coord->Bxy.ydown(); + Bxyup = coord->Bxy().yup(); + Bxydown = coord->Bxy().ydown(); } else { // No 3D metrics // Need to shift to/from field aligned coordinates - g23up = g23down = g23c = toFieldAligned(coord->g23); - g_23up = g_23down = g_23c = toFieldAligned(coord->g_23); - Jup = Jdown = Jc = toFieldAligned(coord->J); - dyup = dydown = dyc = toFieldAligned(coord->dy); - dzup = dzdown = dzc = toFieldAligned(coord->dz); - Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy); + g23up = g23down = g23c = toFieldAligned(coord->g23()); + g_23up = g_23down = g_23c = toFieldAligned(coord->g_23()); + Jup = Jdown = Jc = toFieldAligned(coord->J()); + dyup = dydown = dyc = toFieldAligned(coord->dy()); + dzup = dzdown = dzc = toFieldAligned(coord->dz()); + Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy()); } // Y flux @@ -984,15 +984,15 @@ Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { for (int k = zstart; k <= zend; k++) { BoutReal coef_u = 0.5 - * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) - + coord->g_23(i, j + 1, k) - / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k))); + * (coord->g_23()(i, j, k) / SQ(coord->J()(i, j, k) * coord->Bxy()(i, j, k)) + + coord->g_23()(i, j + 1, k) + / SQ(coord->J()(i, j + 1, k) * coord->Bxy()(i, j + 1, k))); BoutReal coef_d = 0.5 - * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) - + coord->g_23(i, j - 1, k) - / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); + * (coord->g_23()(i, j, k) / SQ(coord->J()(i, j, k) * coord->Bxy()(i, j, k)) + + coord->g_23()(i, j - 1, k) + / SQ(coord->J()(i, j - 1, k) * coord->Bxy()(i, j - 1, k))); const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; const auto zend = bout::build::use_metric_3d ? k : mesh->zend; @@ -1044,10 +1044,10 @@ Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = zstart; k <= zend; k++) { // Coefficient in front of df/dy term - BoutReal coef = - coord->g_23(i, j, k) - / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k)) - / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); + BoutReal coef = coord->g_23(i, j, k) + / (coord->dy()(i, j + 1, k) + 2. * coord->dy()(i, j, k) + + coord->dy()(i, j - 1, k)) + / SQ(coord->J()(i, j, k) * coord->Bxy()(i, j, k)); const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; const auto zend = bout::build::use_metric_3d ? k : mesh->zend; @@ -1117,15 +1117,15 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f Field3D fc = f; Field3D ac = a; - Field3D g23c = coord->g23; - Field3D g_23c = coord->g_23; - Field3D g12c = coord->g12; - Field3D g_12c = coord->g_12; - Field3D Jc = coord->J; - Field3D dxc = coord->dx; - Field3D dyc = coord->dy; - Field3D dzc = coord->dz; - Field3D Bxyc = coord->Bxy; + Field3D g23c = coord->g23(); + Field3D g_23c = coord->g_23(); + Field3D g12c = coord->g12(); + Field3D g_12c = coord->g_12(); + Field3D Jc = coord->J(); + Field3D dxc = coord->dx(); + Field3D dyc = coord->dy(); + Field3D dzc = coord->dz(); + Field3D Bxyc = coord->Bxy(); // Calculate the X derivative at cell edge (X + 1/2), including in Y guard cells // This is used to calculate Y flux contribution from g21 * d/dx @@ -1176,48 +1176,48 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f // 3D Metric, need yup/ydown fields. // Requires previous communication of metrics // -- should insert communication here? - if (!coord->g23.hasParallelSlices() || !coord->g_23.hasParallelSlices() - || !coord->dy.hasParallelSlices() || !coord->dz.hasParallelSlices() - || !coord->Bxy.hasParallelSlices() || !coord->J.hasParallelSlices()) { + if (!coord->g23().hasParallelSlices() || !coord->g_23().hasParallelSlices() + || !coord->dy().hasParallelSlices() || !coord->dz().hasParallelSlices() + || !coord->Bxy().hasParallelSlices() || !coord->J().hasParallelSlices()) { throw BoutException("metrics have no yup/down"); } - g23up = coord->g23.yup(); - g23down = coord->g23.ydown(); + g23up = coord->g23().yup(); + g23down = coord->g23().ydown(); - g_23up = coord->g_23.yup(); - g_23down = coord->g_23.ydown(); + g_23up = coord->g_23().yup(); + g_23down = coord->g_23().ydown(); - g12up = coord->g12.yup(); - g12down = coord->g12.ydown(); + g12up = coord->g12().yup(); + g12down = coord->g12().ydown(); - g_12up = coord->g_12.yup(); - g_12down = coord->g_12.ydown(); + g_12up = coord->g_12().yup(); + g_12down = coord->g_12().ydown(); - Jup = coord->J.yup(); - Jdown = coord->J.ydown(); + Jup = coord->J().yup(); + Jdown = coord->J().ydown(); - dyup = coord->dy.yup(); - dydown = coord->dy.ydown(); + dyup = coord->dy().yup(); + dydown = coord->dy().ydown(); - dzup = coord->dz.yup(); - dzdown = coord->dz.ydown(); + dzup = coord->dz().yup(); + dzdown = coord->dz().ydown(); - Bxyup = coord->Bxy.yup(); - Bxydown = coord->Bxy.ydown(); + Bxyup = coord->Bxy().yup(); + Bxydown = coord->Bxy().ydown(); } else { // No 3D metrics // Need to shift to/from field aligned coordinates - g23up = g23down = g23c = toFieldAligned(coord->g23); - g_23up = g_23down = g_23c = toFieldAligned(coord->g_23); - g12up = g12down = g12c = toFieldAligned(coord->g12); - g_12up = g_12down = g_12c = toFieldAligned(coord->g_12); - Jup = Jdown = Jc = toFieldAligned(coord->J); - dxc = toFieldAligned(coord->dx); - dyup = dydown = dyc = toFieldAligned(coord->dy); - dzup = dzdown = dzc = toFieldAligned(coord->dz); - Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy); + g23up = g23down = g23c = toFieldAligned(coord->g23()); + g_23up = g_23down = g_23c = toFieldAligned(coord->g_23()); + g12up = g12down = g12c = toFieldAligned(coord->g12()); + g_12up = g_12down = g_12c = toFieldAligned(coord->g_12()); + Jup = Jdown = Jc = toFieldAligned(coord->J()); + dxc = toFieldAligned(coord->dx()); + dyup = dydown = dyc = toFieldAligned(coord->dy()); + dzup = dzdown = dzc = toFieldAligned(coord->dz()); + Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy()); } // Y flux @@ -1542,15 +1542,15 @@ Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, for (int k = zstart; k <= zend; k++) { BoutReal coef_u = 0.5 - * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) - + coord->g_23(i, j + 1, k) - / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k))); + * (coord->g_23()(i, j, k) / SQ(coord->J()(i, j, k) * coord->Bxy()(i, j, k)) + + coord->g_23()(i, j + 1, k) + / SQ(coord->J()(i, j + 1, k) * coord->Bxy()(i, j + 1, k))); BoutReal coef_d = 0.5 - * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) - + coord->g_23(i, j - 1, k) - / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); + * (coord->g_23()(i, j, k) / SQ(coord->J()(i, j, k) * coord->Bxy()(i, j, k)) + + coord->g_23()(i, j - 1, k) + / SQ(coord->J()(i, j - 1, k) * coord->Bxy()(i, j - 1, k))); const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; const auto zend = bout::build::use_metric_3d ? k : mesh->zend; @@ -1604,10 +1604,10 @@ Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = zstart; k <= zend; k++) { // Coefficient in front of df/dy term - BoutReal coef = - coord->g_23(i, j, k) - / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k)) - / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); + BoutReal coef = coord->g_23()(i, j, k) + / (coord->dy()(i, j + 1, k) + 2. * coord->dy()(i, j, k) + + coord->dy()(i, j - 1, k)) + / SQ(coord->J()(i, j, k) * coord->Bxy()(i, j, k)); const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; const auto zend = bout::build::use_metric_3d ? k : mesh->zend; diff --git a/src/electromagnetic.cxx b/src/electromagnetic.cxx index 6fce27005..e35c237af 100644 --- a/src/electromagnetic.cxx +++ b/src/electromagnetic.cxx @@ -244,7 +244,7 @@ void Electromagnetic::transform_impl(GuardedOptions& state) { Vector3D A; A.covariant = true; A.x = A.z = 0.0; - A.y = Apar_flutter * (coords->J * coords->Bxy); + A.y = Apar_flutter * (coords->J() * coords->Bxy()); // Perturbed magnetic field vector // Note: Contravariant components (dB_x, dB_y, dB_z) @@ -252,8 +252,8 @@ void Electromagnetic::transform_impl(GuardedOptions& state) { // Set components of the perturbed unit vector // Note: Options can't (yet) contain vectors - set(state["fields"]["deltab_flutter_x"], delta_B.x / coords->Bxy); - set(state["fields"]["deltab_flutter_z"], delta_B.z / coords->Bxy); + set(state["fields"]["deltab_flutter_x"], delta_B.x / coords->Bxy()); + set(state["fields"]["deltab_flutter_z"], delta_B.z / coords->Bxy()); #endif } } diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index 55b8719ce..48b36b70d 100644 --- a/src/evolve_density.cxx +++ b/src/evolve_density.cxx @@ -220,7 +220,7 @@ void EvolveDensity::transform_impl(GuardedOptions& state) { auto* coord = mesh->getCoordinates(); Field3D low_n_coeff = - SQ(coord->dy) * coord->g_22 + SQ(coord->dy()) * coord->g_22() * log(density_floor / clamp(N, 1e-3 * density_floor, density_floor)); low_n_coeff.applyBoundary("neumann"); set(species["low_n_coeff"], low_n_coeff); @@ -306,7 +306,7 @@ void EvolveDensity::finally(const Options& state) { if (hyper_z > 0.) { auto* coord = N.getCoordinates(); - ddt(N) -= hyper_z * SQ(SQ(coord->dz)) * D4DZ4(N); + ddt(N) -= hyper_z * SQ(SQ(coord->dz())) * D4DZ4(N); } // Collect the external source from above with all the sources from diff --git a/src/evolve_energy.cxx b/src/evolve_energy.cxx index cc5f7023c..ee75982b4 100644 --- a/src/evolve_energy.cxx +++ b/src/evolve_energy.cxx @@ -280,7 +280,7 @@ void EvolveEnergy::finally(const Options& state) { if (hyper_z > 0.) { auto* coord = N.getCoordinates(); - ddt(E) -= hyper_z * SQ(SQ(coord->dz)) * D4DZ4(E); + ddt(E) -= hyper_z * SQ(SQ(coord->dz())) * D4DZ4(E); } ////////////////////// diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index 5da1103aa..bdfa6da52 100644 --- a/src/evolve_momentum.cxx +++ b/src/evolve_momentum.cxx @@ -225,7 +225,7 @@ void EvolveMomentum::finally(const Options& state) { if (hyper_z > 0.) { auto* coord = N.getCoordinates(); - ddt(NV) -= hyper_z * SQ(SQ(coord->dz)) * D4DZ4(NV); + ddt(NV) -= hyper_z * SQ(SQ(coord->dz())) * D4DZ4(NV); } // Other sources/sinks diff --git a/src/fieldline_geometry.cxx b/src/fieldline_geometry.cxx index 399bce426..e26fa8668 100644 --- a/src/fieldline_geometry.cxx +++ b/src/fieldline_geometry.cxx @@ -26,7 +26,7 @@ Field3D calculate_Lpar() { BoutReal offset = 0; // Offset to ensure ylow domain boundary starts at 0. This is needed // because the integration starts at the first cell center, not the boundary. - auto dy = coord->dy; // Get the poloidal cell width + auto dy = coord->dy(); // Get the poloidal cell width lpar(0, 0, 0) = 0.5 * dy(0, 0, 0); // Initialize lpar at the first cell center // The parallel length at the first cell center is half of the cell width. // The local trapezoidal integration below only knows the correct starting @@ -214,9 +214,7 @@ FieldlineGeometry::FieldlineGeometry(std::string name, Options& options, Solver* // Bxy is no longer consistent with the Jacobian. // Set equal to NaN, to prevent anyone from using it. - BOUT_FOR(i, coord->Bxy.getRegion("RGN_ALL")) { - coord->Bxy[i] = std::numeric_limits::quiet_NaN(); - } + coord->setBxy(BoutNaN); // Calculate the Jacobian over the whole local field, including the guard // cells. Setting only the interior (ystart..yend) leaves the guard cells @@ -225,18 +223,16 @@ FieldlineGeometry::FieldlineGeometry(std::string name, Options& options, Solver* // inter-processor boundaries degrades the solution (and roughly doubles the // spread between different y-decompositions). effective_magnetic_field_strength // is defined over RGN_ALL, so the guard cells get a consistent value here. - BOUT_FOR(i, effective_magnetic_field_strength.getRegion("RGN_ALL")) { - // N.b. - // The Jacobian has units of [m / radian T], which is why we need an extra factor of Lnorm. - coord->J[i] = 1 / effective_magnetic_field_strength[i] / Lnorm; - } - // Fill the inter-processor guard cells with the neighbour's Jacobian so the - // metric is continuous across rank boundaries in a parallel y-decomposition. - ASSERT2(not coord->J.isFci()); - mesh->communicate(coord->J); + // N.b. + // The Jacobian has units of [m / radian T], which is why we need an extra factor of Lnorm. +#if BOUT_USE_METRIC_3D + coord->setJ(Field3D{1 / effective_magnetic_field_strength / Lnorm}); +#else + coord->setJ(DC(Field3D{1 / effective_magnetic_field_strength / Lnorm})); +#endif // Parallel length of cell - Field3D dlpar = Field3D(coord->dy) / Lnorm; + Field3D dlpar = Field3D(coord->dy()) / Lnorm; // Width of flux tube in the radial direction flux_tube_width = lambda_int * flux_expansion; // Length of the cell in the poloidal direction diff --git a/src/neutral_boundary.cxx b/src/neutral_boundary.cxx index 646eb7d4b..7d7f70ddb 100644 --- a/src/neutral_boundary.cxx +++ b/src/neutral_boundary.cxx @@ -135,17 +135,17 @@ void NeutralBoundary::transform_impl(GuardedOptions& state) { * v_th; // Thermal reflected energy // Cross-sectional area in XZ plane: - BoutReal da = (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) * 0.5 - * (coord->dx[i] + coord->dx[im]) * 0.5 - * (coord->dz[i] + coord->dz[im]); // [m^2] + BoutReal da = (coord->J()[i] + coord->J()[im]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[im])) * 0.5 + * (coord->dx()[i] + coord->dx()[im]) * 0.5 + * (coord->dz()[i] + coord->dz()[im]); // [m^2] // Multiply by area to get energy flow (power) BoutReal flow = q * da; // [W] // Divide by cell volume to get source [W/m^3] BoutReal cooling_source = - flow / (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); + flow / (coord->dx()[i] * coord->dy()[i] * coord->dz()[i] * coord->J()[i]); // Subtract from cell next to boundary energy_source[i] -= cooling_source; @@ -197,17 +197,17 @@ void NeutralBoundary::transform_impl(GuardedOptions& state) { * v_th; // Thermal reflected energy // Cross-sectional area in XZ plane: - BoutReal da = (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])) * 0.5 - * (coord->dx[i] + coord->dx[ip]) * 0.5 - * (coord->dz[i] + coord->dz[ip]); // [m^2] + BoutReal da = (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])) * 0.5 + * (coord->dx()[i] + coord->dx()[ip]) * 0.5 + * (coord->dz()[i] + coord->dz()[ip]); // [m^2] // Multiply by area to get energy flow (power) BoutReal flow = q * da; // [W] // Divide by cell volume to get source [W/m^3] BoutReal cooling_source = - flow / (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); + flow / (coord->dx()[i] * coord->dy()[i] * coord->dz()[i] * coord->J()[i]); // Subtract from cell next to boundary energy_source[i] -= cooling_source; @@ -248,12 +248,12 @@ void NeutralBoundary::transform_impl(GuardedOptions& state) { // Expanded form of the calculation for clarity // Converts dy to poloidal length: dl = dy / sqrt(g22) = dy * h_theta - BoutReal dpol = 0.5 * (coord->dy[i] + coord->dy[ig]) * 1 - / (0.5 * (sqrt(coord->g22[i]) + sqrt(coord->g22[ig]))); + BoutReal dpol = 0.5 * (coord->dy()[i] + coord->dy()[ig]) * 1 + / (0.5 * (sqrt(coord->g22()[i]) + sqrt(coord->g22()[ig]))); // Converts dz to toroidal length: = dz*sqrt(g_33) = dz * R = 2piR - BoutReal dtor = 0.5 * (coord->dz[i] + coord->dz[ig]) * 0.5 - * (sqrt(coord->g_33[i]) + sqrt(coord->g_33[ig])); + BoutReal dtor = 0.5 * (coord->dz()[i] + coord->dz()[ig]) * 0.5 + * (sqrt(coord->g_33()[i]) + sqrt(coord->g_33()[ig])); BoutReal da = dpol * dtor; // [m^2] @@ -261,9 +261,9 @@ void NeutralBoundary::transform_impl(GuardedOptions& state) { BoutReal flow = q * da; // [W] // Divide by cell volume to get source [W/m^3] - BoutReal cooling_source = - flow - / (coord->J[i] * coord->dx[i] * coord->dy[i] * coord->dz[i]); // [W m^-3] + BoutReal cooling_source = flow + / (coord->J()[i] * coord->dx()[i] * coord->dy()[i] + * coord->dz()[i]); // [W m^-3] // Subtract from cell next to boundary energy_source[i] -= cooling_source; @@ -307,12 +307,12 @@ void NeutralBoundary::transform_impl(GuardedOptions& state) { // Expanded form of the calculation for clarity // Converts dy to poloidal length: dl = dy / sqrt(g22) = dy * h_theta - BoutReal dpol = 0.5 * (coord->dy[i] + coord->dy[ig]) * 1 - / (0.5 * (sqrt(coord->g22[i]) + sqrt(coord->g22[ig]))); + BoutReal dpol = 0.5 * (coord->dy()[i] + coord->dy()[ig]) * 1 + / (0.5 * (sqrt(coord->g22()[i]) + sqrt(coord->g22()[ig]))); // Converts dz to toroidal length: = dz*sqrt(g_33) = dz * R = 2piR - BoutReal dtor = 0.5 * (coord->dz[i] + coord->dz[ig]) * 0.5 - * (sqrt(coord->g_33[i]) + sqrt(coord->g_33[ig])); + BoutReal dtor = 0.5 * (coord->dz()[i] + coord->dz()[ig]) * 0.5 + * (sqrt(coord->g_33()[i]) + sqrt(coord->g_33()[ig])); BoutReal da = dpol * dtor; // [m^2] @@ -320,9 +320,9 @@ void NeutralBoundary::transform_impl(GuardedOptions& state) { BoutReal flow = q * da; // [W] // Divide by cell volume to get source [W/m^3] - BoutReal cooling_source = - flow - / (coord->J[i] * coord->dx[i] * coord->dy[i] * coord->dz[i]); // [W m^-3] + BoutReal cooling_source = flow + / (coord->J()[i] * coord->dx()[i] * coord->dy()[i] + * coord->dz()[i]); // [W m^-3] // Subtract from cell next to boundary energy_source[i] -= cooling_source; diff --git a/src/neutral_full_velocity.cxx b/src/neutral_full_velocity.cxx index d5d285a74..6b59d054c 100644 --- a/src/neutral_full_velocity.cxx +++ b/src/neutral_full_velocity.cxx @@ -347,7 +347,7 @@ void NeutralFullVelocity::transform_impl(GuardedOptions& state) { // Vn2D is covariant and b = e_y / (JB) to write: // // V_{||n} = b dot V_n = Vn2D.y / (JB) - Vnpar = Vn2D.y / (coord->J * coord->Bxy); + Vnpar = Vn2D.y / (coord->J() * coord->Bxy()); // Set values in the state auto localstate = state["species"][objectName()]; @@ -596,19 +596,20 @@ void NeutralFullVelocity::finally(const Options& state) { if (localstate.isSet("momentum_source")) { Snv = DC(get(localstate["momentum_source"])); Field2D Fpar_mN = Snv / (AA * Nn2D_floor); - ddt(Vn2D).y += Fpar_mN * (coord->J * coord->Bxy); // Parallel flow + ddt(Vn2D).y += Fpar_mN * (coord->J() * coord->Bxy()); // Parallel flow if (toroidal_flow) { - ddt(Vn2D).z += Fpar_mN * coord->g_23 / (coord->J * coord->Bxy); // Toroidal flow + ddt(Vn2D).z += + Fpar_mN * coord->g_23() / (coord->J() * coord->Bxy()); // Toroidal flow } // NOTE: Should we add the contribution of Sn here? // Sn is introduced in the momentum equation // because we solve for Vn instead of AA*Nn*Vn // It is propably something like that: - // ddt(Vn2D).y += Vn2D.y * Sn / Nn2D_floor * (coord->J * coord->Bxy); // Parallel flow + // ddt(Vn2D).y += Vn2D.y * Sn / Nn2D_floor * (coord->J() * coord->Bxy()); // Parallel flow // if (toroidal_flow) { - // ddt(Vn2D).z += Vn2D.z * Sn / Nn2D_floor * coord->g_23 / (coord->J * coord->Bxy); + // ddt(Vn2D).z += Vn2D.z * Sn / Nn2D_floor * coord->g_23() / (coord->J() * coord->Bxy()); // // Toroidal flow // } @@ -622,7 +623,7 @@ void NeutralFullVelocity::finally(const Options& state) { // Radial flow ddt(Vn2D).x -= Vn2D.x * collision_freq; // Binormal flow - ddt(Vn2D).z -= (Vn2D.z - (coord->g_23 / coord->g_22) * Vn2D.y) * collision_freq; + ddt(Vn2D).z -= (Vn2D.z - (coord->g_23() / coord->g_22()) * Vn2D.y) * collision_freq; } ////////////////////////////////////////////////////// diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 61fd6a0ab..e4324e076 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -118,8 +118,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas heat conduction?") .withDefault(true); - conduction_method = - options["conduction_method"].withDefault(conduction_method); + conduction_method = options["conduction_method"].withDefault(conduction_method); collisionality_override = options["collisionality_override"] diff --git a/src/polarisation_drift.cxx b/src/polarisation_drift.cxx index 716c4aed7..358a5f7dd 100644 --- a/src/polarisation_drift.cxx +++ b/src/polarisation_drift.cxx @@ -23,7 +23,7 @@ PolarisationDrift::PolarisationDrift(std::string name, Options& alloptions, // Cache the B^2 value auto coord = mesh->getCoordinates(); - Bsq = SQ(coord->Bxy); + Bsq = SQ(coord->Bxy()); phiSolver = Laplacian::create(&options["laplacian"]); diff --git a/src/recalculate_metric.cxx b/src/recalculate_metric.cxx index f64dd171e..57e698bec 100644 --- a/src/recalculate_metric.cxx +++ b/src/recalculate_metric.cxx @@ -1,8 +1,8 @@ +#include #include +#include #include #include -#include -#include #include "../include/recalculate_metric.hxx" @@ -13,13 +13,13 @@ void recalculate_metric(BoutReal Lnorm, BoutReal Bnorm) { Field2D Rxy, Bpxy, Btxy, hthe, sinty; GRID_LOAD5(Rxy, Bpxy, Btxy, hthe, sinty); // Load metrics - Coordinates *coord = mesh->getCoordinates(); - + Coordinates* coord = mesh->getCoordinates(); + // Checking for dpsi and qinty used in BOUT grids 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"; @@ -28,14 +28,15 @@ void recalculate_metric(BoutReal Lnorm, BoutReal Bnorm) { 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)); // Calculate metric components - if (Options::root()["mesh"]["paralleltransform"]["type"].as() == "shifted") { + if (Options::root()["mesh"]["paralleltransform"]["type"].as() + == "shifted") { sinty = 0.0; // I disappears from metric } @@ -44,21 +45,23 @@ void recalculate_metric(BoutReal Lnorm, BoutReal Bnorm) { sbp = -1.0; } - 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 Field2D g11 = SQ(Rxy * Bpxy); + const Field2D g22 = 1.0 / SQ(hthe); + const auto g33 = SQ(sinty) * g11 + SQ(coord->Bxy()) / g11; + const Field2D g12 = 0.0; + const Field2D g13 = -sinty * g11; + const Field2D g23 = -sbp * Btxy / (hthe * Bpxy * Rxy); - coord->J = hthe / Bpxy; + const Field2D g_11 = 1.0 / g11 + SQ(sinty * Rxy); + const auto g_22 = SQ(coord->Bxy() * hthe / Bpxy); + const Field2D g_33 = Rxy * Rxy; + const Field2D g_12 = sbp * Btxy * hthe * sinty * Rxy / Bpxy; + const Field2D g_13 = sinty * Rxy * Rxy; + const Field2D 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(); + const Field2D J = hthe / Bpxy; + coord->setJ(J); } diff --git a/src/recycling.cxx b/src/recycling.cxx index 81a2a053f..14ff55210 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -197,11 +197,11 @@ void Recycling::transform_impl(GuardedOptions& state) { // Get metric tensor components Coordinates* coord = mesh->getCoordinates(); - const auto& J = coord->J; - const auto& dy = coord->dy; - const auto& dx = coord->dx; - const auto& dz = coord->dz; - const auto& g_22 = coord->g_22; + const auto& J = coord->J(); + const auto& dy = coord->dy(); + const auto& dx = coord->dx(); + const auto& dz = coord->dz(); + const auto& g_22 = coord->g_22(); for (auto& channel : channels) { const GuardedOptions species_from = state["species"][channel.from]; @@ -529,12 +529,13 @@ void Recycling::transform_impl(GuardedOptions& state) { // Calculate radial wall area in [m^2] // Calculate final cell volume [m^3] BoutReal dpolsheath = - 0.5 * (coord->dy[i] + coord->dy[ig]) * 1 - / (0.5 * (sqrt(coord->g22[i]) + sqrt(coord->g22[ig]))); - BoutReal dtorsheath = 0.5 * (coord->dz[i] + coord->dz[ig]) * 0.5 - * (sqrt(coord->g_33[i]) + sqrt(coord->g_33[ig])); + 0.5 * (coord->dy()[i] + coord->dy()[ig]) * 1 + / (0.5 * (sqrt(coord->g22()[i]) + sqrt(coord->g22()[ig]))); + BoutReal dtorsheath = 0.5 * (coord->dz()[i] + coord->dz()[ig]) * 0.5 + * (sqrt(coord->g_33()[i]) + sqrt(coord->g_33()[ig])); BoutReal dasheath = dpolsheath * dtorsheath; // [m^2] - BoutReal dv = coord->J[i] * coord->dx[i] * coord->dy[i] * coord->dz[i]; + BoutReal dv = + coord->J()[i] * coord->dx()[i] * coord->dy()[i] * coord->dz()[i]; // Calculate particle and energy fluxes of neutrals hitting the pump // Assume thermal velocity greater than perpendicular velocity and use it @@ -653,12 +654,14 @@ void Recycling::transform_impl(GuardedOptions& state) { // Calculate radial wall area in [m^2] // Calculate final cell volume [m^3] BoutReal dpolsheath = - 0.5 * (coord->dy[i] + coord->dy[ig]) * 1 - / (0.5 * (sqrt(coord->g22[i]) + sqrt(coord->g22[ig]))); - BoutReal dtorsheath = 0.5 * (coord->dz[i] + coord->dz[ig]) * 0.5 - * (sqrt(coord->g_33[i]) + sqrt(coord->g_33[ig])); + 0.5 * (coord->dy()[i] + coord->dy()[ig]) * 1 + / (0.5 * (sqrt(coord->g22()[i]) + sqrt(coord->g22()[ig]))); + BoutReal dtorsheath = + 0.5 * (coord->dz()[i] + coord->dz()[ig]) * 0.5 + * (sqrt(coord->g_33()[i]) + sqrt(coord->g_33()[ig])); BoutReal dasheath = dpolsheath * dtorsheath; // [m^2] - BoutReal dv = coord->J[i] * coord->dx[i] * coord->dy[i] * coord->dz[i]; + BoutReal dv = + coord->J()[i] * coord->dx()[i] * coord->dy()[i] * coord->dz()[i]; // Calculate particle and energy fluxes of neutrals hitting the pump // Assume thermal velocity greater than perpendicular velocity and use diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index 0fc905d5d..109b7295e 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -224,9 +224,9 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so Curlb_B.y *= SQ(Lnorm); Curlb_B.z *= SQ(Lnorm); - Curlb_B *= 2. / coord->Bxy; + Curlb_B *= 2. / coord->Bxy(); - Bsq = SQ(coord->Bxy); + Bsq = SQ(coord->Bxy()); diagnose = options["diagnose"].doc("Output additional diagnostics?").withDefault(false); @@ -720,7 +720,8 @@ void RelaxPotential::finally(const Options& state) { // Div_par(jpar) = B * Grad_par(jpar / B) // Using the approximation for small delta-B/B // b dot Grad(jpar) = Grad_par(jpar) + [jpar, Apar] - ddt(Vort) += coord->Bxy * bracket(jpar / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA); + ddt(Vort) += + coord->Bxy() * bracket(jpar / coord->Bxy(), Apar_flutter, BRACKET_ARAKAWA); } } @@ -746,7 +747,7 @@ void RelaxPotential::finally(const Options& state) { if (hyper_z > 0) { // Form of hyper-viscosity to suppress zig-zags in Z auto* coord = Vort.getCoordinates(); - ddt(Vort) -= hyper_z * SQ(SQ(coord->dz)) * D4DZ4(Vort); + ddt(Vort) -= hyper_z * SQ(SQ(coord->dz())) * D4DZ4(Vort); } if (phi_sheath_dissipation) { diff --git a/src/sheath_boundary.cxx b/src/sheath_boundary.cxx index fd79625fd..d864137a7 100644 --- a/src/sheath_boundary.cxx +++ b/src/sheath_boundary.cxx @@ -566,11 +566,11 @@ void SheathBoundary::transform_impl(GuardedOptions& state) { q = std::min(q, 0.0); // Multiply by cell area to get power - const BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + const BoutReal flux = q * (coord->J()[i] + coord->J()[im]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[im])); // Divide by volume of cell to get energy loss rate (< 0) - const BoutReal power = flux / (coord->dy[i] * coord->J[i]); + const BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); #if CHECKLEVEL >= 1 if (!std::isfinite(power)) { @@ -633,11 +633,11 @@ void SheathBoundary::transform_impl(GuardedOptions& state) { * nesheath * vesheath; q = std::max(q, 0.0); // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + BoutReal flux = q * (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])); // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); #if CHECKLEVEL >= 1 if (!std::isfinite(power)) { throw BoutException( @@ -773,11 +773,11 @@ void SheathBoundary::transform_impl(GuardedOptions& state) { q = std::min(q, 0.0); // Multiply by cell area to get power - const BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + const BoutReal flux = q * (coord->J()[i] + coord->J()[im]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[im])); // Divide by volume of cell to get energy loss rate (< 0) - const BoutReal power = flux / (coord->dy[i] * coord->J[i]); + const BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); ASSERT1(std::isfinite(power)); ASSERT2(power <= 0.0); @@ -874,11 +874,11 @@ void SheathBoundary::transform_impl(GuardedOptions& state) { q = std::max(q, 0.0); // Multiply by cell area to get power - const BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + const BoutReal flux = q * (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])); // Divide by volume of cell to get energy loss rate (> 0) - const BoutReal power = flux / (coord->dy[i] * coord->J[i]); + const BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); ASSERT1(std::isfinite(power)); ASSERT2(power >= 0.0); diff --git a/src/sheath_boundary_insulating.cxx b/src/sheath_boundary_insulating.cxx index 7d08d90e5..481f57b8a 100644 --- a/src/sheath_boundary_insulating.cxx +++ b/src/sheath_boundary_insulating.cxx @@ -333,11 +333,11 @@ void SheathBoundaryInsulating::transform_impl(GuardedOptions& state) { q = std::min(q, 0.0); // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + BoutReal flux = q * (coord->J()[i] + coord->J()[im]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[im])); // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); ASSERT1(std::isfinite(power)); ASSERT2(power <= 0.0); @@ -412,11 +412,11 @@ void SheathBoundaryInsulating::transform_impl(GuardedOptions& state) { q = std::max(q, 0.0); // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + BoutReal flux = q * (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])); // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); ASSERT1(std::isfinite(power)); ASSERT2(power >= 0.0); @@ -473,11 +473,11 @@ void SheathBoundaryInsulating::transform_impl(GuardedOptions& state) { * nesheath * vesheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + BoutReal flux = q * (coord->J()[i] + coord->J()[im]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[im])); // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); #if CHECKLEVEL >= 1 if (!std::isfinite(power)) { @@ -509,11 +509,11 @@ void SheathBoundaryInsulating::transform_impl(GuardedOptions& state) { * nesheath * vesheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + BoutReal flux = q * (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])); // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = flux / (coord->dy()[i] * coord->J()[i]); #if CHECKLEVEL >= 1 if (!std::isfinite(power)) { throw BoutException( diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index fb0828a1b..668e8dc11 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -394,11 +394,12 @@ void SheathBoundarySimple::transform_impl(GuardedOptions& state) { q -= (2.5 * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; // Cross-sectional area in XZ plane and cell volume - BoutReal da = (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) * 0.5 - * (coord->dx[i] + coord->dx[im]) * 0.5 - * (coord->dz[i] + coord->dz[im]); // [m^2] - BoutReal dv = (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] + BoutReal da = (coord->J()[i] + coord->J()[im]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[im])) * 0.5 + * (coord->dx()[i] + coord->dx()[im]) * 0.5 + * (coord->dz()[i] + coord->dz()[im]); // [m^2] + BoutReal dv = + (coord->dx()[i] * coord->dy()[i] * coord->dz()[i] * coord->J()[i]); // [m^3] // Get power and energy source BoutReal heatflow = q * da; // [W] @@ -459,11 +460,12 @@ void SheathBoundarySimple::transform_impl(GuardedOptions& state) { q -= (2.5 * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; // Cross-sectional area in XZ plane and cell volume - BoutReal da = (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])) * 0.5 - * (coord->dx[i] + coord->dx[ip]) * 0.5 - * (coord->dz[i] + coord->dz[ip]); // [m^2] - BoutReal dv = (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] + BoutReal da = (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])) * 0.5 + * (coord->dx()[i] + coord->dx()[ip]) * 0.5 + * (coord->dz()[i] + coord->dz()[ip]); // [m^2] + BoutReal dv = + (coord->dx()[i] * coord->dy()[i] * coord->dz()[i] * coord->J()[i]); // [m^3] // Get power and energy source BoutReal heatflow = q * da; // [W] @@ -607,12 +609,12 @@ void SheathBoundarySimple::transform_impl(GuardedOptions& state) { q -= (2.5 * tisheath + 0.5 * Mi * SQ(visheath)) * nisheath * visheath; // Cross-sectional area in XZ plane and cell volume - BoutReal da = (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])) * 0.5 - * (coord->dx[i] + coord->dx[ip]) * 0.5 - * (coord->dz[i] + coord->dz[ip]); // [m^2] + BoutReal da = (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])) * 0.5 + * (coord->dx()[i] + coord->dx()[ip]) * 0.5 + * (coord->dz()[i] + coord->dz()[ip]); // [m^2] BoutReal dv = - (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] + (coord->dx()[i] * coord->dy()[i] * coord->dz()[i] * coord->J()[i]); // [m^3] // Get power and energy source BoutReal heatflow = q * da; // [W] @@ -677,12 +679,12 @@ void SheathBoundarySimple::transform_impl(GuardedOptions& state) { q -= (2.5 * tisheath + 0.5 * Mi * SQ(visheath)) * nisheath * visheath; // Cross-sectional area in XZ plane and cell volume - BoutReal da = (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])) * 0.5 - * (coord->dx[i] + coord->dx[ip]) * 0.5 - * (coord->dz[i] + coord->dz[ip]); // [m^2] + BoutReal da = (coord->J()[i] + coord->J()[ip]) + / (sqrt(coord->g_22()[i]) + sqrt(coord->g_22()[ip])) * 0.5 + * (coord->dx()[i] + coord->dx()[ip]) * 0.5 + * (coord->dz()[i] + coord->dz()[ip]); // [m^2] BoutReal dv = - (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] + (coord->dx()[i] * coord->dy()[i] * coord->dz()[i] * coord->J()[i]); // [m^3] // Get power and energy source BoutReal heatflow = q * da; // [W] diff --git a/src/snb_conduction.cxx b/src/snb_conduction.cxx index 0bce5f876..71378f96e 100644 --- a/src/snb_conduction.cxx +++ b/src/snb_conduction.cxx @@ -13,14 +13,15 @@ void SNBConduction::transform_impl(GuardedOptions& state) { // SNB non-local heat flux. Also returns the Spitzer-Harm value for comparison // Note: Te in eV, Ne in Nnorm - auto dy_orig = mesh->getCoordinates()->dy; - mesh->getCoordinates()->dy *= rho_s0; // Convert distances to m + auto dy_orig = mesh->getCoordinates()->dy(); + ASSERT2(not dy_orig.isFci()); + dy_orig *= rho_s0; // Convert distances to m // Inputs in eV and m^-3 Div_Q_SNB = snb.divHeatFlux(Te, Ne, &Div_Q_SH); // Restore the metric tensor - mesh->getCoordinates()->dy = dy_orig; + mesh->getCoordinates()->setDy(dy_orig); // Normalise from eV/m^3/s Div_Q_SNB /= Tnorm * Nnorm * Omega_ci; diff --git a/src/sound_speed.cxx b/src/sound_speed.cxx index c4ce23fc9..6104c9246 100644 --- a/src/sound_speed.cxx +++ b/src/sound_speed.cxx @@ -25,7 +25,8 @@ void SoundSpeed::transform_impl(GuardedOptions& state) { auto AA = get(species["AA"]); // Atomic mass number if (species.isSet("density")) { - total_density += GET_NOBOUNDARY(Field3D, species["density"]) * get(species["AA"]); + total_density += + GET_NOBOUNDARY(Field3D, species["density"]) * get(species["AA"]); } if (species.isSet("temperature")) { @@ -45,13 +46,13 @@ void SoundSpeed::transform_impl(GuardedOptions& state) { } if (alfven_wave) { - auto *coord = fastest_wave.getCoordinates(); + auto* coord = fastest_wave.getCoordinates(); for (auto& i : fastest_wave.getRegion("RGN_NOBNDRY")) { - BoutReal alfven_speed = beta_norm * coord->Bxy[i] / sqrt(total_density[i]); + BoutReal alfven_speed = beta_norm * coord->Bxy()[i] / sqrt(total_density[i]); fastest_wave[i] = BOUTMAX(fastest_wave[i], alfven_speed); } } set(state["sound_speed"], sound_speed); - set(state["fastest_wave"], fastest_wave*fastest_wave_factor); + set(state["fastest_wave"], fastest_wave * fastest_wave_factor); } diff --git a/src/vorticity.cxx b/src/vorticity.cxx index c723cc4af..abcbe3439 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -64,8 +64,8 @@ Vector3D Grad_perp_XZ(const Field3D& f) { result.y = 0.0; result.z = emptyFrom(f); - auto dx = metric->dx; - auto dz = metric->dz; + auto dx = metric->dx(); + auto dz = metric->dz(); BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) { auto xp = i.xp(); @@ -190,11 +190,11 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) // Create an XY solver for n=0 component laplacexy = LaplaceXY::create(mesh); // Set coefficients for Boussinesq solve - laplacexy->setCoefs(average_atomic_mass / SQ(DC(coord->Bxy)), 0.0); + laplacexy->setCoefs(average_atomic_mass / SQ(DC(coord->Bxy())), 0.0); } phiSolver = Laplacian::create(&options["laplacian"]); // Set coefficients for Boussinesq solve - phiSolver->setCoefC(average_atomic_mass / SQ(coord->Bxy)); + phiSolver->setCoefC(average_atomic_mass / SQ(coord->Bxy())); if (phi_boundary_relax) { // Set the last update time to -1, so it will reset @@ -249,9 +249,9 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) Curlb_B.y *= SQ(Lnorm); Curlb_B.z *= SQ(Lnorm); - Curlb_B *= 2. / coord->Bxy; + Curlb_B *= 2. / coord->Bxy(); - Bsq = SQ(coord->Bxy); + Bsq = SQ(coord->Bxy()); diagnose = options["diagnose"].doc("Output additional diagnostics?").withDefault(false); @@ -862,7 +862,8 @@ void Vorticity::finally(const Options& state) { // Div_par(jpar) = B * Grad_par(jpar / B) // Using the approximation for small delta-B/B // b dot Grad(jpar) = Grad_par(jpar) + [jpar, Apar] - ddt(Vort) += coord->Bxy * bracket(jpar / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA); + ddt(Vort) += + coord->Bxy() * bracket(jpar / coord->Bxy(), Apar_flutter, BRACKET_ARAKAWA); } } @@ -882,7 +883,7 @@ void Vorticity::finally(const Options& state) { if (hyper_z > 0) { // Form of hyper-viscosity to suppress zig-zags in Z const auto* coord = Vort.getCoordinates(); - ddt(Vort) -= hyper_z * SQ(SQ(coord->dz)) * D4DZ4(Vort); + ddt(Vort) -= hyper_z * SQ(SQ(coord->dz())) * D4DZ4(Vort); } if (phi_sheath_dissipation) { diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 1be1cdd43..fb86eb539 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -90,7 +91,7 @@ public: void setCoordinates(std::shared_ptr coords, CELL_LOC location = CELL_CENTRE) { - coords_map[location] = coords; + coords_map[location] = std::move(coords); } void setGridDataSource(GridDataSource* source_in) { source = source_in; } @@ -366,3 +367,54 @@ public: private: Options values; ///< Store values to be returned by get() }; + +// A mock ParallelTransform to test transform_from_field_aligned +// property of FieldFactory. For now, the transform just returns the +// negative of the input. Ideally, this will get moved to GoogleMock +// when we start using it. +// +// Can turn off the ability to do the transform. Should still be valid +class MockParallelTransform : public ParallelTransform { +public: + MockParallelTransform(Mesh& mesh, bool allow_transform_) + : ParallelTransform(mesh), allow_transform(allow_transform_) {} + ~MockParallelTransform() = default; + + void calcParallelSlices(Field3D&) override {} + + bool canToFromFieldAligned() const override { return allow_transform; } + + bool requiresTwistShift(bool, YDirectionType) override { return false; } + + void checkInputGrid() override {} + + Field3D fromFieldAligned(const Field3D& f, const std::string&) override { + if (f.getDirectionY() != YDirectionType::Aligned) { + throw BoutException("Unaligned field passed to fromFieldAligned"); + } + return -f; + } + + FieldPerp fromFieldAligned(const FieldPerp& f, const std::string&) override { + if (f.getDirectionY() != YDirectionType::Aligned) { + throw BoutException("Unaligned field passed to fromFieldAligned"); + } + return -f; + } + + Field3D toFieldAligned(const Field3D& f, const std::string&) override { + if (f.getDirectionY() != YDirectionType::Standard) { + throw BoutException("Aligned field passed to toFieldAligned"); + } + return -f; + } + FieldPerp toFieldAligned(const FieldPerp& f, const std::string&) override { + if (f.getDirectionY() != YDirectionType::Standard) { + throw BoutException("Aligned field passed to toFieldAligned"); + } + return -f; + } + +private: + const bool allow_transform; +}; diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 2758dbe41..021f37a24 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -3,11 +3,14 @@ #include #include +#include "bout/build_config.hxx" +#include "bout/metric_tensor.hxx" #include #include #include #include #include +#include #include #include #include @@ -31,7 +34,7 @@ /// Use this template class directly to use different sized grid: /// /// using MyTest = FakeMeshFixture_tmpl<7, 9, 11>; -template +template class FakeMeshFixture_tmpl : public ::testing::Test { public: FakeMeshFixture_tmpl() @@ -47,22 +50,22 @@ public: Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}); - - // Set some auxilliary variables - // Usually set in geometry() - // Note: For testing these are set to non-zero values - test_coords->G1 = test_coords->G2 = test_coords->G3 = 0.1; + mesh_m.setCoordinates(test_coords); // Set nonuniform corrections - test_coords->non_uniform = true; - test_coords->d1_dx = test_coords->d1_dy = 0.2; - test_coords->d1_dz = 0.0; -#if BOUT_USE_METRIC_3D - test_coords->Bxy.splitParallelSlices(); - test_coords->Bxy.yup() = test_coords->Bxy.ydown() = test_coords->Bxy; -#endif - - // No call to Coordinates::geometry() needed here + test_coords->setNon_uniform(true); + test_coords->setD1_dx(0.2); + test_coords->setD1_dy(0.2); + test_coords->setD1_dz(0.0); + + if (bout::build::use_metric_3d) { + bout::FieldMetric mutable_Bxy = test_coords->Bxy(); + mutable_Bxy.splitParallelSlices(); + mutable_Bxy.yup() = test_coords->Bxy(); + mutable_Bxy.ydown() = test_coords->Bxy(); + test_coords->setBxy(mutable_Bxy); + } + mesh_m.setCoordinates(test_coords); mesh_m.setGridDataSource(new FakeGridDataSource()); // May need a ParallelTransform to create fields, because create3D calls @@ -89,22 +92,25 @@ public: Field2D{1.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}); - - // Set some auxilliary variables - test_coords_staggered->G1 = test_coords_staggered->G2 = test_coords_staggered->G3 = - 0.1; + mesh_staggered_m.setCoordinates(test_coords_staggered); // Set nonuniform corrections - test_coords_staggered->non_uniform = true; - test_coords_staggered->d1_dx = test_coords_staggered->d1_dy = 0.2; - test_coords_staggered->d1_dz = 0.0; -#if BOUT_USE_METRIC_3D - test_coords_staggered->Bxy.splitParallelSlices(); - test_coords_staggered->Bxy.yup() = test_coords_staggered->Bxy.ydown() = - test_coords_staggered->Bxy; -#endif - - // No call to Coordinates::geometry() needed here + test_coords_staggered->setNon_uniform(true); + test_coords_staggered->setD1_dx(0.2); + test_coords_staggered->setD1_dy(0.2); + test_coords_staggered->setD1_dz(0.0); + + if (bout::build::use_metric_3d) { + bout::FieldMetric mutable_Bxy = test_coords_staggered->Bxy(); + mutable_Bxy.splitParallelSlices(); + test_coords_staggered->setBxy(mutable_Bxy); + + mutable_Bxy = test_coords_staggered->Bxy(); + mutable_Bxy.yup() = test_coords_staggered->Bxy(); + mutable_Bxy.ydown() = test_coords_staggered->Bxy(); + test_coords_staggered->setBxy(mutable_Bxy); + } + test_coords_staggered->setParallelTransform( bout::utils::make_unique(mesh_staggered_m)); @@ -113,6 +119,11 @@ public: mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_XLOW); mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_YLOW); mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_ZLOW); + + if constexpr (FCI) { + mesh_m.getCoordinates()->setParallelTransform( + bout::utils::make_unique(mesh_m, false)); + } } FakeMeshFixture_tmpl(const FakeMeshFixture_tmpl&) = delete; @@ -148,3 +159,4 @@ public: }; using FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>; +using FakeMeshFixtureFCI = FakeMeshFixture_tmpl<3, 5, 7, true>; diff --git a/tests/unit/test_anomalous_diffusion.cxx b/tests/unit/test_anomalous_diffusion.cxx index 561757953..1858d2a0b 100644 --- a/tests/unit/test_anomalous_diffusion.cxx +++ b/tests/unit/test_anomalous_diffusion.cxx @@ -45,7 +45,7 @@ TEST_F(AnomalousDiffusionTest, NoDiffusion) { TEST_F(AnomalousDiffusionTest, ParticleDiffusion) { Coordinates* coords = mesh->getCoordinates(); - coords->Bxy = 1.0; // Note: This is non-finite or zero? + coords->setBxy(1.0); // Note: This is non-finite or zero? Options options; options["units"]["meters"] = 1.0; @@ -74,7 +74,7 @@ TEST_F(AnomalousDiffusionTest, ParticleDiffusion) { "RGN_NOBNDRY")); // Expect the sum over all cells of density source to be zero - auto dV = coords->J * coords->dx * coords->dy * coords->dz; // Cell volume + auto dV = coords->cell_volume(); // Cell volume Field3D source = get(state["species"]["h"]["density_source"]); BoutReal integral = 0.0; diff --git a/tests/unit/test_fieldline_geometry.cxx b/tests/unit/test_fieldline_geometry.cxx index c7d5c42ba..3325124e1 100644 --- a/tests/unit/test_fieldline_geometry.cxx +++ b/tests/unit/test_fieldline_geometry.cxx @@ -132,8 +132,8 @@ TEST_F(FieldlineGeometryTest, SetsCoordinatesJacobianAndBxy) { // Bxy is no longer consistent with the new Jacobian, so should be NaN // everywhere to stop anyone accidentally using it. - for (auto i : coord->Bxy.getRegion("RGN_ALL")) { - ASSERT_TRUE(std::isnan(coord->Bxy[i])); + for (auto i : coord->Bxy().getRegion("RGN_ALL")) { + ASSERT_TRUE(std::isnan(coord->Bxy()[i])); } // J = 1 / (Beff) / Lnorm. With constant inputs, transport_broadening = 1 @@ -236,7 +236,7 @@ TEST_F(FieldlineGeometryTest, GeometryFactorsAreSelfConsistentForNonTrivialProfi component.outputVars(outputs); Coordinates* coord = mesh->getCoordinates(); - auto dy = coord->dy; + auto dy = coord->dy(); auto pitch_angle = get(outputs["fieldline_geometry_magnetic_pitch"]); auto Rxy = get(outputs["fieldline_geometry_Rxy"]); auto lambda_int = get(outputs["fieldline_geometry_lambda_int"]);