diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index d100e2d797..d8ffcc65fe 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -1,17 +1,11 @@ /************************************************************************** * Describes coordinate systems * - * ChangeLog - * ========= - * - * 2014-11-10 Ben Dudson - * * Created by separating metric from Mesh - * * ************************************************************************** - * Copyright 2014 B.D.Dudson + * Copyright 2014 - 2024 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -239,6 +233,8 @@ public: FieldMetric& J() const; ///< Magnitude of B = nabla z times nabla x + ///< Note: This should always be positive + ///< for both right- and left-handed coordinates. const FieldMetric& Bxy() const { return Bxy_; } void setJ(const FieldMetric& J); @@ -389,6 +385,19 @@ public: const FieldMetric& Grad2_par2_DDY_invSg(CELL_LOC outloc, const std::string& method) const; + /// Calculate 1 / (J |B|) + /// This is cached as it is used frequently in parallel operators. + /// + /// In a Clebsch coordinate system + /// B = (1 / J) e_y + /// so the unit vector along the magnetic field is: + /// b = B / |B| = (1 / ( J |B| )) e_y + /// so e.g. + /// Grad_par = b dot Grad = (1 / J |B|) d / dy. + /// + /// Note: In a right-handed Clebsch coordinate system + /// this is 1 / sqrt(g_22) i.e. positive. + /// In a left-handed coordinate system it is negative. const FieldMetric& invSg() const; ChristoffelSymbols& christoffel_symbols(); @@ -466,6 +475,8 @@ private: FieldMetric getUnaligned(const std::string& name, BoutReal default_value); + /// Recalculate magnetic field magnitude Bxy = |B| + /// Note: Always positive FieldMetric recalculateBxy() const; /// Non-uniform meshes. Need to use DDX, DDY diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d4df92b275..5d49b6443f 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -905,7 +905,8 @@ MetricTensor::FieldMetric Coordinates::recalculateJacobian() const { } MetricTensor::FieldMetric Coordinates::recalculateBxy() const { - return sqrt(g_22()) / J(); + // Note: J may be negative, by return is always positive + return sqrt(g_22()) / abs(J()); } void Coordinates::setParallelTransform(Options* mesh_options) { @@ -1452,7 +1453,7 @@ GValues& Coordinates::g_values() const { const Coordinates::FieldMetric& Coordinates::invSg() const { if (invSgCache == nullptr) { auto ptr = std::make_unique(); - (*ptr) = 1.0 / sqrt(g_22()); + (*ptr) = 1.0 / (J() * Bxy()); invSgCache = std::move(ptr); } return *invSgCache; @@ -1502,6 +1503,7 @@ void Coordinates::setJ(const FieldMetric& J) { void Coordinates::setBxy(FieldMetric Bxy) { //TODO: Calculate Bxy and check value is close + // Also check that Bxy is positive Bxy_ = std::move(Bxy); localmesh->communicate(Bxy_); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 78c1ef35a1..e793c85fcf 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -2,9 +2,9 @@ * Various differential operators defined on BOUT grid * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2024 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -104,7 +104,9 @@ Field3D Grad_parP(const Field3D& apar, const Field3D& f) { for (int x = 1; x <= mesh->LocalNx - 2; x++) { for (int y = mesh->ystart; y <= mesh->yend; y++) { for (int z = 0; z < ncz; z++) { - BoutReal by = 1. / sqrt(metric->g_22(x, y, z)); + + // Note: by is negative in a left-handed coordinate system + BoutReal by = 1. / (metric->J(x, y, z) * metric->Bxy(x, y, z)); // Z indices zm and zp int const zm = (z - 1 + ncz) % ncz; int const zp = (z + 1) % ncz; @@ -258,14 +260,13 @@ Field3D Div_par(const Field3D& f, const Field3D& v) { BoutReal const vR = 0.5 * (v(i, j, k) + v.yup()(i, j + 1, k)); // Calculate flux at right boundary (y+1/2) + // Note: Magnitude of B at the midpoint used rather than J/sqrt(g_22) BoutReal const fluxRight = - fR * vR * (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + fR * vR * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j + 1, k)); // Calculate at left boundary (y-1/2) BoutReal const fluxLeft = - fL * vL * (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + fL * vL * 2 / (coord->Bxy(i, j, k) + coord->Bxy(i, j - 1, k)); result(i, j, k) = (fluxRight - fluxLeft) / (coord->dy(i, j, k) * coord->J(i, j, k)); @@ -285,7 +286,7 @@ Field3D Div_par_flux(const Field3D& v, const Field3D& f, CELL_LOC outloc, auto Bxy_floc = f.getCoordinates()->Bxy(); if (!f.hasParallelSlices()) { - return metric->Bxy() * FDDY(v, f / Bxy_floc, outloc, method) / sqrt(metric->g_22()); + return FDDY(v, f / Bxy_floc, outloc, method) / metric->J(); } // Need to modify yup and ydown fields @@ -294,7 +295,7 @@ Field3D Div_par_flux(const Field3D& v, const Field3D& f, CELL_LOC outloc, f_B.splitParallelSlices(); f_B.yup() = f.yup() / Bxy_floc; f_B.ydown() = f.ydown() / Bxy_floc; - return metric->Bxy() * FDDY(v, f_B, outloc, method) / sqrt(metric->g_22()); + return FDDY(v, f_B, outloc, method) / metric->J(); } Field3D Div_par_flux(const Field3D& v, const Field3D& f, const std::string& method, @@ -478,7 +479,8 @@ Coordinates::FieldMetric b0xGrad_dot_Grad(const Field2D& phi, const Field2D& A, // Upwind A using these velocities Coordinates::FieldMetric result = VDDX(vx, A, outloc) + VDDY(vy, A, outloc); - result /= metric->J() * sqrt(metric->g_22()); + result /= SQ(metric->J()) + * metric->Bxy(); // J^2 B same sign for left & right handed coordinates ASSERT1(result.getLocation() == outloc); @@ -519,7 +521,7 @@ Field3D b0xGrad_dot_Grad(const Field2D& phi, const Field3D& A, CELL_LOC outloc) Field3D result = VDDX(vx, A, outloc) + VDDY(vy, A, outloc) + VDDZ(vz, A, outloc); - result /= (metric->J() * sqrt(metric->g_22())); + result /= SQ(metric->J()) * metric->Bxy(); // J^2 B #if BOUT_USE_TRACK result.name = "b0xGrad_dot_Grad(" + phi.name + "," + A.name + ")"; @@ -554,7 +556,7 @@ Field3D b0xGrad_dot_Grad(const Field3D& p, const Field2D& A, CELL_LOC outloc) { Field3D result = VDDX(vx, A, outloc) + VDDY(vy, A, outloc); - result /= (metric->J() * sqrt(metric->g_22())); + result /= SQ(metric->J()) * metric->Bxy(); // J^2 B #if BOUT_USE_TRACK result.name = "b0xGrad_dot_Grad(" + p.name + "," + A.name + ")"; @@ -595,7 +597,7 @@ Field3D b0xGrad_dot_Grad(const Field3D& phi, const Field3D& A, CELL_LOC outloc) Field3D result = VDDX(vx, A, outloc) + VDDY(vy, A, outloc) + VDDZ(vz, A, outloc); - result /= (metric->J() * sqrt(metric->g_22())); + result /= SQ(metric->J()) * metric->Bxy(); // J^2 B #if BOUT_USE_TRACK result.name = "b0xGrad_dot_Grad(" + phi.name + "," + A.name + ")";