diff --git a/.pip_install_for_tests.sh b/.pip_install_for_tests.sh index 443fdf7db..c9df6fe2a 100755 --- a/.pip_install_for_tests.sh +++ b/.pip_install_for_tests.sh @@ -19,9 +19,8 @@ do fi done -# Install Hermes-3 -git clone https://github.com/boutproject/xhermes.git -cd xhermes -pip3 install --user . -cd .. +# Install xBOUT to make sure we get latest master version +pip3 install --user git+https://github.com/boutproject/xBOUT.git@v0.3.7 +# Install xHermes for Hermes-3 python +pip3 install --user git+https://github.com/boutproject/xhermes.git diff --git a/CMakeLists.txt b/CMakeLists.txt index 907f3fef5..666b98a31 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -357,6 +357,28 @@ if(HERMES_TESTS) # Have hermes_unit_tests look for dynamic libs in ./ before searching elsewhere set_target_properties(hermes_unit_tests PROPERTIES BUILD_RPATH "$ORIGIN") endif() + +# Method of manufactured solutions (MMS) tests of operators + option(HERMES_DIVOPS_MMS_TESTS "Build the differential operator MMS tests" ON) + if(HERMES_DIVOPS_MMS_TESTS) + #TEST_SOURCES= + add_executable(hermes_mms_operator_tests + tests/mms_operator/main.cxx) + #${TEST_SOURCES}) + target_link_libraries(hermes_mms_operator_tests PRIVATE gtest bout++::bout++ hermes-3-lib) + add_test(NAME differential_operators_mms_orthogonal + WORKING_DIRECTORY tests/mms_operator/ + COMMAND python3 orthogonal_test.py) + set_tests_properties(differential_operators_mms_orthogonal PROPERTIES LABELS integration) + add_test(NAME differential_operators_mms_nonorthogonal + WORKING_DIRECTORY tests/mms_operator/ + COMMAND python3 nonorthogonal_test.py) + set_tests_properties(differential_operators_mms_nonorthogonal PROPERTIES LABELS integration) + add_test(NAME rhs_mms_neutral_mixed + WORKING_DIRECTORY tests/mms_operator/ + COMMAND python3 neutral_mixed_ddt_test.py) + set_tests_properties(rhs_mms_neutral_mixed PROPERTIES LABELS integration) + endif() endif() # Compile-time options diff --git a/include/classical_diffusion.hxx b/include/classical_diffusion.hxx index cb738cc90..a87817a94 100644 --- a/include/classical_diffusion.hxx +++ b/include/classical_diffusion.hxx @@ -15,6 +15,11 @@ private: Field3D Dn; ///< Particle diffusion coefficient BoutReal custom_D; ///< User-set particle diffusion coefficient override + // Flow diagnostics + Field3D cls_pf_perp_xlow, cls_pf_perp_ylow; + Field3D cls_mf_perp_xlow, cls_mf_perp_ylow; + Field3D cls_nef_perp_xlow, cls_nef_perp_ylow; + Field3D cls_tef_perp_xlow, cls_tef_perp_ylow; void transform_impl(GuardedOptions& state) override; }; diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 5ed2841ef..ce2b58f7a 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -76,6 +76,16 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, const Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow, bool bndry_flux = true); +/*! + * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity + * + * This version includes corrections for non-orthogonal meshes + * in which the g12 and g13 components can be non-zero + * i.e. X-Y, X-Z and Y-Z coordinates can all be non-orthogonal. + */ +Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& x, Field3D& flux_xlow, + Field3D& flux_ylow); + namespace FV { /// Superbee limiter diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 37979b9a2..4242f9f4e 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -18,7 +18,7 @@ struct NeutralMixed : public Component { /// @param options Top-level options. Settings will be taken from options[name] /// @param solver Time-integration solver to be used NeutralMixed(const std::string& name, Options& options, Solver *solver); - + /// Use the final simulation state to update internal state /// (e.g. time derivatives) void finally(const Options &state) override; @@ -30,7 +30,7 @@ struct NeutralMixed : public Component { void precon(const Options &state, BoutReal gamma) override; private: std::string name; ///< Species name - + Field3D Nn, Pn, NVn; // Density, pressure and parallel momentum Field3D Vn; ///< Neutral parallel velocity Field3D Tn; ///< Neutral temperature @@ -45,6 +45,7 @@ private: Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators BoutReal flux_limit; ///< Diffusive flux limit BoutReal diffusion_limit; ///< Maximum diffusion coefficient + BoutReal neutral_lmax; bool sheath_ydown, sheath_yup; @@ -53,12 +54,15 @@ private: BoutReal pressure_floor; ///< Minimum Pn used when dividing Pn by Nn to get Tn. bool freeze_low_density; ///< Freeze evolution in low density regions? - + BoutReal collisionality_override; ///< Rnn input for testing + BoutReal density_norm, pressure_norm; ///< Normalisations + BoutReal momentum_norm; ///< Normalisations bool neutral_viscosity; ///< include viscosity? bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? - + bool normalise_sources; ///< Normalise input sources? + Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity bool precondition {true}; ///< Enable preconditioner? @@ -66,7 +70,7 @@ private: bool lax_flux; ///< Use Lax flux for advection terms std::unique_ptr inv; ///< Laplacian inversion used for preconditioning - Field3D density_source, pressure_source; ///< External input source + Field3D density_source, pressure_source, momentum_source; ///< External input source Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source Field3D sound_speed; ///< Sound speed for use with Lax flux diff --git a/src/classical_diffusion.cxx b/src/classical_diffusion.cxx index 5761c7f98..267618df0 100644 --- a/src/classical_diffusion.cxx +++ b/src/classical_diffusion.cxx @@ -1,5 +1,6 @@ #include "classical_diffusion.hxx" +#include "../include/div_ops.hxx" #include ClassicalDiffusion::ClassicalDiffusion(std::string name, Options& alloptions, Solver*) @@ -36,29 +37,30 @@ void ClassicalDiffusion::transform_impl(GuardedOptions& state) { // The only term here comes from the resistive drift Field3D Ptotal = 0.0; - for (auto& kv : allspecies.getChildren()) { - const auto species = kv.second; + if (!(custom_D > 0)) { + for (auto& kv : allspecies.getChildren()) { + const auto species = kv.second; - if (!(species.isSet("charge") and species.isSet("pressure"))) { - continue; // Skip, go to next species - } - auto q = get(species["charge"]); - if (fabs(q) < 1e-5) { - continue; + if (!(species.isSet("charge") and species.isSet("pressure"))) { + continue; // Skip, go to next species + } + auto q = get(species["charge"]); + if (fabs(q) < 1e-5) { + continue; + } + Ptotal += GET_VALUE(Field3D, species["pressure"]); } - Ptotal += GET_VALUE(Field3D, species["pressure"]); } - auto electrons = allspecies["e"]; - const auto me = get(electrons["AA"]); - const Field3D Ne = GET_VALUE(Field3D, electrons["density"]); - // Particle diffusion coefficient. Applied to all charged species // so that net transport is ambipolar if (custom_D > 0) { // User-set Dn = custom_D; } else { // Calculated from collisions + auto electrons = allspecies["e"]; + const auto me = get(electrons["AA"]); + const Field3D Ne = GET_VALUE(Field3D, electrons["density"]); const Field3D nu_e = floor(GET_VALUE(Field3D, electrons["collision_frequency"]), 1e-10); Dn = floor(Ptotal, 1e-5) * me * nu_e / (floor(Ne, 1e-5) * Bsq); } @@ -70,40 +72,49 @@ void ClassicalDiffusion::transform_impl(GuardedOptions& state) { for (auto kv : allspecies.getChildren()) { GuardedOptions species = allspecies[kv.first]; // Note: Need non-const - - if (!(species.isSet("charge") and species.isSet("density"))) { - continue; // Skip, go to next species - } - auto q = get(species["charge"]); - if (fabs(q) < 1e-5) { - continue; + if (!(custom_D > 0)) { + if (!(species.isSet("charge") and species.isSet("density"))) { + continue; // Skip, go to next species + } + auto q = get(species["charge"]); + if (fabs(q) < 1e-5) { + continue; + } } - const auto N = GET_VALUE(Field3D, species["density"]); - add(species["density_source"], FV::Div_a_Grad_perp(Dn, N)); + // add(species["density_source"], FV::Div_a_Grad_perp(Dn, N)); + add(species["density_source"], + Div_a_Grad_perp_nonorthog(Dn, N, cls_pf_perp_xlow, cls_pf_perp_ylow)); if (IS_SET(species["velocity"])) { const auto V = GET_VALUE(Field3D, species["velocity"]); const auto AA = GET_VALUE(BoutReal, species["AA"]); - add(species["momentum_source"], FV::Div_a_Grad_perp(Dn * AA * V, N)); + // add(species["momentum_source"], FV::Div_a_Grad_perp(Dn * AA * V, N)); + add(species["momentum_source"], + Div_a_Grad_perp_nonorthog(Dn * AA * V, N, cls_mf_perp_xlow, cls_mf_perp_ylow)); } if (IS_SET(species["temperature"])) { const auto T = GET_VALUE(Field3D, species["temperature"]); - add(species["energy_source"], FV::Div_a_Grad_perp(Dn * (3. / 2) * T, N)); - - // Cross-field heat conduction - // kappa_perp = 2 * n * nu_ii * rho_i^2 - - const auto P = GET_VALUE(Field3D, species["pressure"]); - const auto AA = GET_VALUE(BoutReal, species["AA"]); + // add(species["energy_source"], FV::Div_a_Grad_perp(Dn * (3. / 2) * T, N)); + add(species["energy_source"], + Div_a_Grad_perp_nonorthog(Dn * (3. / 2) * T, N, cls_nef_perp_xlow, + cls_nef_perp_ylow)); // TODO: Figure out what to do with the below if(custom_D < 0) { + // Cross-field heat conduction + // kappa_perp = 2 * n * nu_ii * rho_i^2 + const auto P = GET_VALUE(Field3D, species["pressure"]); + const auto AA = GET_VALUE(BoutReal, species["AA"]); const Field3D nu = floor(GET_VALUE(Field3D, species["collision_frequency"]), 1e-10); - add(species["energy_source"], FV::Div_a_Grad_perp(2. * floor(P, 1e-5) * nu * AA / Bsq, T)); + // add(species["energy_source"], FV::Div_a_Grad_perp(2. * floor(P, 1e-5) * nu * AA + // / Bsq, T)); + add(species["energy_source"], + Div_a_Grad_perp_nonorthog(2. * floor(P, 1e-5) * nu * AA / Bsq, T, + cls_tef_perp_xlow, cls_tef_perp_ylow)); } } } diff --git a/src/div_ops.cxx b/src/div_ops.cxx index ed9a61dce..adcaf8572 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -36,12 +36,12 @@ using bout::globals::mesh; -const Field3D Div_par_diffusion_index(const Field3D &f, bool bndry_flux) { +const Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux) { Field3D result; result = 0.0; - Coordinates *coord = mesh->getCoordinates(); - + Coordinates* coord = mesh->getCoordinates(); + for (int i = mesh->xstart; i <= mesh->xend; i++) for (int j = mesh->ystart - 1; j <= mesh->yend; j++) for (int k = 0; k < mesh->LocalNz; k++) { @@ -54,8 +54,7 @@ const Field3D Div_par_diffusion_index(const Field3D &f, bool bndry_flux) { if ((j == mesh->ystart - 1) && mesh->firstY(i)) continue; } - BoutReal J = - 0.5 * (coord->J(i, j) + coord->J(i, j + 1)); // Jacobian at boundary + BoutReal J = 0.5 * (coord->J(i, j) + coord->J(i, j + 1)); // Jacobian at boundary BoutReal gradient = f(i, j + 1, k) - f(i, j, k); @@ -105,13 +104,12 @@ void MC(Stencil1D& n) { * poloidal - If true, includes X-Y flows * positive - If true, limit advected quantity (n_in) to be positive */ -const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, - bool bndry_flux, bool poloidal, - bool positive) { +const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flux, + bool poloidal, bool positive) { Field3D result{0.0}; - Coordinates *coord = mesh->getCoordinates(); - + Coordinates* coord = mesh->getCoordinates(); + ////////////////////////////////////////// // X-Z advection. // @@ -138,24 +136,25 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // 1) Interpolate stream function f onto corners fmp, fpp, fpm - BoutReal fmm = 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + - f(i - 1, j, km)); - BoutReal fmp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) + - f(i - 1, j, kp)); // 2nd order accurate - BoutReal fpp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + - f(i + 1, j, kp)); - BoutReal fpm = 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + - f(i + 1, j, km)); + BoutReal fmm = + 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + f(i - 1, j, km)); + BoutReal fmp = 0.25 + * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) + + f(i - 1, j, kp)); // 2nd order accurate + BoutReal fpp = + 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + f(i + 1, j, kp)); + BoutReal fpm = + 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + f(i + 1, j, km)); // 2) Calculate velocities on cell faces BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx - BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) / - coord->dz(i, j); // J*df/dz - BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) / - coord->dz(i, j); // J*df/dz + BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) + / coord->dz(i, j); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) + / coord->dz(i, j); // J*df/dz // 3) Calculate n on the cell faces. The sign of the // velocity determines which side is used. @@ -184,8 +183,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)); } result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); } } else { // Not at a boundary @@ -193,8 +191,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Flux out into next cell BoutReal flux = vR * s.R; result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); // if(i==mesh->xend) // output.write("Setting flux (%d,%d) : %e\n", @@ -219,8 +216,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)); } result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += - flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); } } else { // Not at a boundary @@ -228,8 +224,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, if (vL < 0.0) { BoutReal flux = vL * s.L; result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += - flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); } } @@ -299,11 +294,11 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Average dfdy to right X boundary BoutReal f_R = - 0.5 * ((coord->g11(i + 1, j) * coord->g23(i + 1, j) / - SQ(coord->Bxy(i + 1, j))) * - dfdy(i + 1, j, k) + - (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) * - dfdy(i, j, k)); + 0.5 + * ((coord->g11(i + 1, j) * coord->g23(i + 1, j) / SQ(coord->Bxy(i + 1, j))) + * dfdy(i + 1, j, k) + + (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) + * dfdy(i, j, k)); // Advection velocity across cell face BoutReal Vx = 0.5 * (coord->J(i + 1, j) + coord->J(i, j)) * f_R; @@ -312,8 +307,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, BoutReal flux = Vx; if (Vx > 0) { // Right boundary of cell (i,j,k) - BoutReal nval = - n(i, j, k) + 0.25 * (n(i + 1, j, k) - n(i - 1, j, k)); + BoutReal nval = n(i, j, k) + 0.25 * (n(i + 1, j, k) - n(i - 1, j, k)); if (positive && (nval < 0.0)) { // Limit value to positive nval = 0.0; @@ -321,8 +315,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, flux *= nval; } else { // Left boundary of cell (i+1,j,k) - BoutReal nval = - n(i + 1, j, k) - 0.25 * (n(i + 2, j, k) - n(i, j, k)); + BoutReal nval = n(i + 1, j, k) - 0.25 * (n(i + 2, j, k) - n(i, j, k)); if (positive && (nval < 0.0)) { nval = 0.0; } @@ -330,14 +323,13 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, } result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); } } if (poloidal) { // Y flux - + Field3D dfdx = DDX(f); mesh->communicate(dfdx); dfdx.applyBoundary("neumann"); @@ -345,9 +337,9 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // This calculation is in field aligned coordinates dfdx = toFieldAligned(dfdx); Field3D n_fa = toFieldAligned(n); - + Field3D yresult{zeroFrom(n_fa)}; - + for (int i = mesh->xstart; i <= mesh->xend; i++) { int ys = mesh->ystart - 1; int ye = mesh->yend; @@ -370,16 +362,15 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Average dfdx to upper Y boundary BoutReal f_U = - 0.5 * ((coord->g11(i, j + 1) * coord->g23(i, j + 1) / - SQ(coord->Bxy(i, j + 1))) * - dfdx(i, j + 1, k) + - (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) * - dfdx(i, j, k)); + 0.5 + * ((coord->g11(i, j + 1) * coord->g23(i, j + 1) / SQ(coord->Bxy(i, j + 1))) + * dfdx(i, j + 1, k) + + (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) + * dfdx(i, j, k)); BoutReal Vy = -0.5 * (coord->J(i, j + 1) + coord->J(i, j)) * f_U; - if (mesh->firstY(i) && !mesh->periodicY(i) && - (j == mesh->ystart - 1)) { + if (mesh->firstY(i) && !mesh->periodicY(i) && (j == mesh->ystart - 1)) { // Lower y boundary. Allow flows out of the domain only Vy = std::min(Vy, 0.0); } @@ -415,7 +406,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, } result += fromFieldAligned(yresult); } - + return result; } @@ -438,8 +429,8 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { // o --- gD --- o // - Coordinates *coord = mesh->getCoordinates(); - + Coordinates* coord = mesh->getCoordinates(); + for (int i = mesh->xstart; i <= mesh->xend; i++) for (int j = mesh->ystart; j <= mesh->yend; j++) for (int k = 0; k < mesh->LocalNz; k++) { @@ -457,16 +448,16 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { BoutReal gU = fs(i, j, kp) - fs(i, j, k); // Flow right - BoutReal flux = gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) * - (coord->dx(i + 1, j) + coord->dx(i, j)) * - (as(i + 1, j, k) + as(i, j, k)); + BoutReal flux = gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) + * (coord->dx(i + 1, j) + coord->dx(i, j)) + * (as(i + 1, j, k) + as(i, j, k)); result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) * - (coord->dx(i - 1, j) + coord->dx(i, j)) * - (as(i - 1, j, k) + as(i, j, k)); + flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) + * (coord->dx(i - 1, j) + coord->dx(i, j)) + * (as(i - 1, j, k) + as(i, j, k)); result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); @@ -478,7 +469,7 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { flux = gD * 0.5 * (as(i, j, k) + as(i, j, km)); result(i, j, k) -= flux; } - + return result; } @@ -508,22 +499,22 @@ const Field3D Div_Z_FV_Index(const Field3D &as, const Field3D &fs) { } // *** USED *** -const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { +const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { Field3D result = 0.0; - Coordinates *coord = mesh->getCoordinates(); - + Coordinates* coord = mesh->getCoordinates(); + for (int i = mesh->xstart; i <= mesh->xend; i++) for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { // 3rd derivative at right boundary - BoutReal d3fdx3 = (f(i + 2, j, k) - 3. * f(i + 1, j, k) + - 3. * f(i, j, k) - f(i - 1, j, k)); + BoutReal d3fdx3 = + (f(i + 2, j, k) - 3. * f(i + 1, j, k) + 3. * f(i, j, k) - f(i - 1, j, k)); - BoutReal flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + BoutReal flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) + * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; if (mesh->lastX() && (i == mesh->xend)) { // Boundary @@ -531,15 +522,15 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { if (bndry_flux) { // Use a one-sided difference formula - d3fdx3 = -((16. / 5) * 0.5 * - (f(i + 1, j, k) + f(i, j, k)) // Boundary value f_b - - 6. * f(i, j, k) // f_0 - + 4. * f(i - 1, j, k) // f_1 - - (6. / 5) * f(i - 2, j, k) // f_2 - ); + d3fdx3 = + -((16. / 5) * 0.5 * (f(i + 1, j, k) + f(i, j, k)) // Boundary value f_b + - 6. * f(i, j, k) // f_0 + + 4. * f(i - 1, j, k) // f_1 + - (6. / 5) * f(i - 2, j, k) // f_2 + ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) + * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; } else { // No fluxes through boundary @@ -557,32 +548,30 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { // On an X boundary if (bndry_flux) { - d3fdx3 = -(-(16. / 5) * 0.5 * - (f(i - 1, j, k) + f(i, j, k)) // Boundary value f_b - + 6. * f(i, j, k) // f_0 - - 4. * f(i + 1, j, k) // f_1 - + (6. / 5) * f(i + 2, j, k) // f_2 - ); + d3fdx3 = + -(-(16. / 5) * 0.5 * (f(i - 1, j, k) + f(i, j, k)) // Boundary value f_b + + 6. * f(i, j, k) // f_0 + - 4. * f(i + 1, j, k) // f_1 + + (6. / 5) * f(i + 2, j, k) // f_2 + ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) + * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); - result(i - 1, j, k) += - flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + result(i - 1, j, k) += flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); } } else { // Not on a boundary - d3fdx3 = (f(i + 1, j, k) - 3. * f(i, j, k) + 3. * f(i - 1, j, k) - - f(i - 2, j, k)); + d3fdx3 = + (f(i + 1, j, k) - 3. * f(i, j, k) + 3. * f(i - 1, j, k) - f(i - 2, j, k)); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) + * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); - result(i - 1, j, k) += - flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + result(i - 1, j, k) += flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); } } } @@ -607,56 +596,49 @@ const Field3D D4DZ4_Index(const Field3D& f) { * we would need the corner cell values to take Y derivatives along X edges * */ -const Field2D Laplace_FV(const Field2D &k, const Field2D &f) { +const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { Field2D result; result.allocate(); - Coordinates *coord = mesh->getCoordinates(); - + Coordinates* coord = mesh->getCoordinates(); + for (int i = mesh->xstart; i <= mesh->xend; i++) for (int j = mesh->ystart; j <= mesh->yend; j++) { // Calculate gradients on cell faces - BoutReal gR = (coord->g11(i, j) + coord->g11(i + 1, j)) * - (f(i + 1, j) - f(i, j)) / - (coord->dx(i + 1, j) + coord->dx(i, j)); + BoutReal gR = (coord->g11(i, j) + coord->g11(i + 1, j)) * (f(i + 1, j) - f(i, j)) + / (coord->dx(i + 1, j) + coord->dx(i, j)); - BoutReal gL = (coord->g11(i - 1, j) + coord->g11(i, j)) * - (f(i, j) - f(i - 1, j)) / - (coord->dx(i - 1, j) + coord->dx(i, j)); + BoutReal gL = (coord->g11(i - 1, j) + coord->g11(i, j)) * (f(i, j) - f(i - 1, j)) + / (coord->dx(i - 1, j) + coord->dx(i, j)); - BoutReal gU = (coord->g22(i, j) + coord->g22(i, j + 1)) * - (f(i, j + 1) - f(i, j)) / - (coord->dy(i, j + 1) + coord->dy(i, j)); + BoutReal gU = (coord->g22(i, j) + coord->g22(i, j + 1)) * (f(i, j + 1) - f(i, j)) + / (coord->dy(i, j + 1) + coord->dy(i, j)); - BoutReal gD = (coord->g22(i, j - 1) + coord->g22(i, j)) * - (f(i, j) - f(i, j - 1)) / - (coord->dy(i, j) + coord->dy(i, j - 1)); + BoutReal gD = (coord->g22(i, j - 1) + coord->g22(i, j)) * (f(i, j) - f(i, j - 1)) + / (coord->dy(i, j) + coord->dy(i, j - 1)); // Flow right - BoutReal flux = gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) * - (k(i + 1, j) + k(i, j)); + BoutReal flux = + gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) * (k(i + 1, j) + k(i, j)); result(i, j) = flux / (coord->dx(i, j) * coord->J(i, j)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) * - (k(i - 1, j) + k(i, j)); + flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) * (k(i - 1, j) + k(i, j)); result(i, j) -= flux / (coord->dx(i, j) * coord->J(i, j)); // Flow up - flux = gU * 0.25 * (coord->J(i, j + 1) + coord->J(i, j)) * - (k(i, j + 1) + k(i, j)); + flux = gU * 0.25 * (coord->J(i, j + 1) + coord->J(i, j)) * (k(i, j + 1) + k(i, j)); result(i, j) += flux / (coord->dy(i, j) * coord->J(i, j)); // Flow down - flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j)) * - (k(i, j - 1) + k(i, j)); + flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j)) * (k(i, j - 1) + k(i, j)); result(i, j) -= flux / (coord->dy(i, j) * coord->J(i, j)); } return result; @@ -914,7 +896,7 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { // Calculate flux from i to i+1 const BoutReal gradient = (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) + + coord->J(i + 1, j) * coord->g11(i + 1, j)) * (f(i + 1, j, k) - f(i, j, k)) / (coord->dx(i, j) + coord->dx(i + 1, j)); @@ -1052,6 +1034,363 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { return result; } +// Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity +// Includes nonorthogonal X-Y and X-Z metric corrections +// +Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& flow_xlow, + Field3D& flow_ylow) { + ASSERT2(a.getLocation() == f.getLocation()); + + Mesh* mesh = a.getMesh(); + + Coordinates* coord = f.getCoordinates(); + + // Zero all flows + flow_xlow = 0.0; + flow_ylow = 0.0; + + // Y and Z fluxes require Y derivatives + + // Fields containing values along the magnetic field + Field3D fup(mesh), fdown(mesh); + Field3D aup(mesh), adown(mesh); + + Field3D g23up(mesh), g23down(mesh); + Field3D g_23up(mesh), g_23down(mesh); + Field3D g12up(mesh), g12down(mesh); + Field3D g_12up(mesh), g_12down(mesh); + Field3D Jup(mesh), Jdown(mesh); + Field3D dyup(mesh), dydown(mesh); + Field3D dzup(mesh), dzdown(mesh); + Field3D Bxyup(mesh), Bxydown(mesh); + + // Values on this y slice (centre). + // This is needed because toFieldAligned may modify the field + 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; + + // 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 + Field3D fddx_xhigh(mesh); + fddx_xhigh.allocate(); + for (int i = mesh->xstart - 1; i <= mesh->xend; i++) { + for (int j = mesh->ystart - 1; j <= mesh->yend + 1; + j++) { // Note: Including one guard cell + for (int k = 0; k < mesh->LocalNz; k++) { + fddx_xhigh(i, j, k) = 2. * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + } + } + } + Field3D fddx_xhigh_up(mesh), fddx_xhigh_down(mesh); + + // Result of the Y and Z fluxes + Field3D yzresult(mesh); + yzresult.allocate(); + + if (f.hasParallelSlices() && a.hasParallelSlices()) { + // Both inputs have yup and ydown + + fup = f.yup(); + fdown = f.ydown(); + + aup = a.yup(); + adown = a.ydown(); + + mesh->communicate(fddx_xhigh); + fddx_xhigh_up = fddx_xhigh.yup(); + fddx_xhigh_down = fddx_xhigh.ydown(); + } else { + // At least one input doesn't have yup/ydown fields. + // Need to shift to/from field aligned coordinates + + fup = fdown = fc = toFieldAligned(f); + aup = adown = ac = toFieldAligned(a); + + fddx_xhigh_up = fddx_xhigh_down = toFieldAligned(fddx_xhigh); + + yzresult.setDirectionY(YDirectionType::Aligned); + flow_ylow.setDirectionY(YDirectionType::Aligned); + } + + if (bout::build::use_metric_3d) { + // 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()) { + throw BoutException("metrics have no yup/down"); + } + + g23up = coord->g23.yup(); + g23down = coord->g23.ydown(); + + g_23up = coord->g_23.yup(); + g_23down = coord->g_23.ydown(); + + g12up = coord->g12.yup(); + g12down = coord->g12.ydown(); + + g_12up = coord->g_12.yup(); + g_12down = coord->g_12.ydown(); + + Jup = coord->J.yup(); + Jdown = coord->J.ydown(); + + dyup = coord->dy.yup(); + dydown = coord->dy.ydown(); + + dzup = coord->dz.yup(); + dzdown = coord->dz.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); + } + + // Y flux + // Includes fluxes due to Z derivatives (non-zero g23 metric) + // and due to X derivatives if grid is nonorthogonal (non-zero g12 metric) + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + BoutReal coef_yz = + 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k))); + + BoutReal coef_xy = + 0.5 + * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_12up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k))); + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.5 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / (dzc(i, j, k) + dzup(i, j + 1, k)); + + // Y derivative + BoutReal dfdy = + 2. * (fup(i, j + 1, k) - fc(i, j, k)) / (dyup(i, j + 1, k) + dyc(i, j, k)); + + // X derivative at Y boundary + BoutReal dfdx = 0.25 + * (fddx_xhigh(i - 1, j, k) + fddx_xhigh(i, j, k) + + fddx_xhigh_up(i - 1, j + 1, k) + fddx_xhigh_up(i, j + 1, k)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * ( + // Component of flux from d/dz and d/dy + (Jc(i, j, k) * g23c(i, j, k) + Jup(i, j + 1, k) * g23up(i, j + 1, k)) + * (dfdz - coef_yz * dfdy) + // Non-orthogonal mesh correction with g12 metric + + (Jc(i, j, k) * g12c(i, j, k) + Jup(i, j + 1, k) * g12up(i, j + 1, k)) + * (dfdx - coef_xy * dfdy)); + + yzresult(i, j, k) = fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Calculate flux between j and j-1 + coef_yz = + 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23down(i, j - 1, k) / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k))); + + coef_xy = + 0.5 + * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_12down(i, j - 1, k) / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k))); + + dfdz = 0.5 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) + / (dzc(i, j, k) + dzdown(i, j - 1, k)); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (dyc(i, j, k) + dydown(i, j - 1, k)); + + dfdx = 0.25 + * (fddx_xhigh(i - 1, j, k) + fddx_xhigh(i, j, k) + + fddx_xhigh_down(i - 1, j - 1, k) + fddx_xhigh_down(i, j - 1, k)); + + fout = + 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * ( + // Component of flux from d/dz and d/dy + (Jc(i, j, k) * g23c(i, j, k) + Jdown(i, j - 1, k) * g23down(i, j - 1, k)) + * (dfdz - coef_yz * dfdy) + // Non-orthogonal mesh correction with g12 metric + + (Jc(i, j, k) * g12c(i, j, k) + + Jdown(i, j - 1, k) * g12down(i, j - 1, k)) + * (dfdx - coef_xy * dfdy)); + + yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + } + } + } + + // Calculate the Y derivative, including in X guard cells + // This is used to calculate X flux contribution from g12 * d/dy + Field3D fddy(mesh); + fddy.allocate(); + fddy.setDirectionY(yzresult.getDirectionY()); + for (int i = mesh->xstart - 1; i <= mesh->xend + 1; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + fddy(i, j, k) = + (fup(i, j + 1, k) - fdown(i, j - 1, k)) + / (0.5 * dyup(i, j + 1, k) + dyc(i, j, k) + 0.5 * dydown(i, j - 1, k)); + } + } + } + + // Z flux + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; + + BoutReal ddx = + 0.5 + * ((fc(i + 1, j, k) - fc(i - 1, j, k)) + / (0.5 * dxc(i + 1, j, k) + dxc(i, j, k) + 0.5 * dxc(i - 1, j, k)) + + (fc(i + 1, j, kp) - fc(i - 1, j, kp)) + / (0.5 * dxc(i + 1, j, kp) + dxc(i, j, kp) + + 0.5 * dxc(i - 1, j, kp))); + + BoutReal ddy = + 0.5 + * ((fup(i, j + 1, k) - fdown(i, j - 1, k)) + / (0.5 * dyup(i, j + 1, k) + dyc(i, j, k) + 0.5 * dydown(i, j - 1, k)) + + (fup(i, j + 1, kp) - fdown(i, j - 1, kp)) + / (0.5 * dyup(i, j + 1, kp) + dyc(i, j, kp) + + 0.5 * dydown(i, j - 1, kp))); + + BoutReal ddz = 2. * (fc(i, j, kp) - fc(i, j, k)) / (dzc(i, j, k) + dzc(i, j, kp)); + + BoutReal fout = + 0.5 * (ac(i, j, k) + ac(i, j, kp)) + * ( + // Jg^xz (d/dx - g_xy / g_yy d/dy) + 0.5 + * (Jc(i, j, k) * coord->g13(i, j, k) + + Jc(i, j, kp) * coord->g13(i, j, kp)) + * (ddx + - ddy * 0.5 + * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_12c(i, j, kp) / SQ(Jc(i, j, kp) * Bxyc(i, j, kp)))) + // Jg^zz (d/dz - g_yz / g_yy d/dy) + + 0.5 + * (Jc(i, j, k) * coord->g33(i, j, k) + + Jc(i, j, kp) * coord->g33(i, j, kp)) + * (ddz + - ddy * 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23c(i, j, kp) + / SQ(Jc(i, j, kp) * Bxyc(i, j, kp))))); + + yzresult(i, j, k) += fout / (Jc(i, j, k) * dzc(i, j, k)); + yzresult(i, j, kp) -= fout / (Jc(i, j, kp) * dzc(i, j, kp)); + } + } + } + + // Check if we need to transform back + Field3D result; + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result = yzresult; + } else { + result = fromFieldAligned(yzresult); + fddy = fromFieldAligned(fddy); + } + + // Flux in X + + for (int i = mesh->xstart - 1; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux from i to i+1 + + const int kp = (k + 1) % mesh->LocalNz; + const int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + BoutReal ddx = 2 * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + BoutReal ddy = 0.5 * (fddy(i, j, k) + fddy(i + 1, j, k)); + + BoutReal ddz = 0.5 + * ((f(i, j, kp) - f(i, j, km)) + / (0.5 * coord->dz(i, j, kp) + coord->dz(i, j, k) + + 0.5 * coord->dz(i, j, km)) + + (f(i + 1, j, kp) - f(i + 1, j, km)) + / (0.5 * coord->dz(i + 1, j, kp) + coord->dz(i + 1, j, k) + + 0.5 * coord->dz(i + 1, j, km))); + + BoutReal fout = + 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * ( + // Jg^xx (d/dx - g_xy / g_yy d/dy) + 0.5 + * (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + * (ddx + - ddy * 0.5 + * (coord->g_12(i, j, k) / coord->g_22(i, j, k) + + coord->g_12(i + 1, j, k) / coord->g_22(i + 1, j, k))) + // Jg^xz (d/dz - g_yz / g_yy d/dy) + + 0.5 + * ((coord->J(i, j, k) * coord->g13(i, j, k) + + coord->J(i + 1, j, k) * coord->g13(i + 1, j, k)) + * (ddz + - ddy * 0.5 + * (coord->g_23(i, j, k) / coord->g_22(i, j, k) + + coord->g_23(i + 1, j, k) + / coord->g_22(i + 1, j, k))))); + + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); + } + } + } + + return result; +} + /// Div ( a Grad_perp(f) ) -- diffusion /// /// Returns the flows in the final arguments @@ -1211,7 +1550,6 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, BoutReal coef = coord->g_23(i, j) / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) / SQ(coord->J(i, j) * coord->Bxy(i, j)); - for (int k = 0; k < mesh->LocalNz; k++) { // Calculate flux between k and k+1 int kp = (k + 1) % mesh->LocalNz; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index a05e0bef3..a88b8c9d5 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -27,6 +27,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* const BoutReal Nnorm = units["inv_meters_cubed"]; const BoutReal Tnorm = units["eV"]; const BoutReal Omega_ci = 1. / units["seconds"].as(); + const BoutReal Cs0 = sqrt(SI::qe * Tnorm / SI::Mp); // Need to take derivatives in X for cross-field diffusion terms ASSERT0(mesh->xstart > 0); @@ -80,6 +81,8 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Enable stabilising lax flux?") .withDefault(true); + neutral_lmax = 0.1 / meters; // Normalised length + flux_limit = options["flux_limit"] .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") @@ -98,6 +101,15 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas heat conduction?") .withDefault(true); + collisionality_override = + options["collisionality_override"] + .doc( + "Paramter for overriding the neutral collision frequency in Dn for testing") + .withDefault(-1.0); + + normalise_sources = options["normalise_sources"] + .doc("Normalise input sources?") + .withDefault(true); diffusion_collisions_mode = options["diffusion_collisions_mode"] .doc("Can be multispecies: all enabled collisions excl. IZ, or afn: CX, IZ and NN collisions") .withDefault("multispecies"); @@ -118,6 +130,15 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* AA = options["AA"].doc("Particle atomic mass. Proton = 1").withDefault(1.0); + if (normalise_sources) { + density_norm = Nnorm * Omega_ci; + pressure_norm = SI::qe * Nnorm * Tnorm * Omega_ci; + momentum_norm = SI::Mp * Nnorm * Cs0 * Omega_ci; + } else { + density_norm = 1.0; + pressure_norm = 1.0; + momentum_norm = 1.0; + } // Try to read the density source from the mesh // Units of particles per cubic meter per second density_source = 0.0; @@ -127,7 +148,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* alloptions[std::string("N") + name]["source"] .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) .withDefault(density_source) - / (Nnorm * Omega_ci); + / density_norm; // Try to read the pressure source from the mesh // Units of Pascals per second @@ -138,7 +159,17 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc(std::string("Source term in ddt(P") + name + std::string("). Units [N/m^2/s]")) .withDefault(pressure_source) - / (SI::qe * Nnorm * Tnorm * Omega_ci); + / pressure_norm; + // Try to read the momentum source from the mesh + momentum_source = 0.0; + mesh->get(momentum_source, fmt::format("NV{}_src", name)); + // Allow the user to override the source + momentum_source = + alloptions[fmt::format("NV{}", name)]["source"] + .doc(fmt::format("Source term in ddt(NV{}). Units [kg m^-2 s^-2]", name)) + .withDefault(momentum_source) + / momentum_norm; + // need some normalisation convention here // Set boundary condition defaults: Neumann for all but the diffusivity. // The dirichlet on diffusivity ensures no radial flux. @@ -193,9 +224,6 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { Vn = NVn / (AA * Nnlim); Vn.applyBoundary("neumann"); - Pnlim = softFloor(Pn, pressure_floor); - Pnlim.applyBoundary(); - ///////////////////////////////////////////////////// // Parallel boundary conditions TRACE("Neutral boundary conditions"); @@ -221,7 +249,6 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { // Zero-gradient pressure Pn(r.ind, mesh->ystart - 1, jz) = Pn(r.ind, mesh->ystart, jz); - Pnlim(r.ind, mesh->ystart - 1, jz) = Pnlim(r.ind, mesh->ystart, jz); // No flow into wall Vn(r.ind, mesh->ystart - 1, jz) = -Vn(r.ind, mesh->ystart, jz); @@ -246,7 +273,6 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { // Zero-gradient pressure Pn(r.ind, mesh->yend + 1, jz) = Pn(r.ind, mesh->yend, jz); - Pnlim(r.ind, mesh->yend + 1, jz) = Pnlim(r.ind, mesh->yend, jz); // No flow into wall Vn(r.ind, mesh->yend + 1, jz) = -Vn(r.ind, mesh->yend, jz); @@ -268,6 +294,19 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { void NeutralMixed::finally(const Options& state) { auto& localstate = state["species"][name]; + // extract auxiliary variables derived from + // Nn, Pn, NVn, from the local state + // and set boundary conditions on evolved quantities + Tn = get(localstate["temperature"]); + Vn = get(localstate["velocity"]); + Pn.setBoundaryTo(get(localstate["pressure"])); + Nn.setBoundaryTo(get(localstate["density"])); + if (!evolve_momentum) { + // momentum is not evolved, so need to get the value from the localstate + NVn = get(localstate["momentum"]); + } else { + NVn.setBoundaryTo(get(localstate["momentum"])); + } // Logarithms used to calculate perpendicular velocity // V_perp = -Dnn * ( Grad_perp(Nn)/Nn + Grad_perp(Tn)/Tn ) // @@ -276,91 +315,93 @@ void NeutralMixed::finally(const Options& state) { // Field3D logNn = log(Nn); // Field3D logTn = log(Tn); + // Nnlim Used where division by neutral density is needed + Nnlim = softFloor(Nn, density_floor); + // Tnlim used where positivity of Tn is required + Field3D Tnlim = softFloor(Tn, temperature_floor); + // Pnlim used where positivity of Pn is required + Pnlim = softFloor(Pn, pressure_floor); logPnlim = log(Pnlim); logPnlim.applyBoundary(); - /////////////////////////////////////////////////////// // Calculate cross-field diffusion from collision frequency // // - Field3D Tnlim = softFloor(Tn, temperature_floor); - - BoutReal neutral_lmax = - 0.1 / get(state["units"]["meters"]); // Normalised length - Field3D Rnn = sqrt(Tnlim / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] - - if (localstate.isSet("collision_frequency")) { - - // Collisionality - // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ - if (collision_names.empty()) { // Calculate only once - at the beginning - - if (diffusion_collisions_mode == "afn") { - for (const auto& collision : localstate["collision_frequencies"].getChildren()) { - - std::string collision_name = collision.second.name(); - - if (// Charge exchange - (collisionSpeciesMatch( - collision_name, name, "+", "cx", "partial")) or - // Ionisation - (collisionSpeciesMatch( - collision_name, name, "+", "iz", "partial")) or - // Neutral-neutral collisions - (collisionSpeciesMatch( - collision_name, name, name, "coll", "exact"))) { - collision_names.push_back(collision_name); - } + if (collisionality_override > 0.0) { + // user has set an override for collision frequency + Dnn = (Tn / AA) / collisionality_override; + } else { + if (localstate.isSet("collision_frequency")) { + // Collisionality + // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ + if (collision_names.empty()) { // Calculate only once - at the beginning + + if (diffusion_collisions_mode == "afn") { + for (const auto& collision : + localstate["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if ( // Charge exchange + (collisionSpeciesMatch(collision_name, name, "+", "cx", "partial")) or + // Ionisation + (collisionSpeciesMatch(collision_name, name, "+", "iz", "partial")) or + // Neutral-neutral collisions + (collisionSpeciesMatch(collision_name, name, name, "coll", "exact"))) { + collision_names.push_back(collision_name); + } + } + // Multispecies mode: all collisions and CX are included + } else if (diffusion_collisions_mode == "multispecies") { + for (const auto& collision : + localstate["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if ( // Charge exchange + (collisionSpeciesMatch(collision_name, name, "", "cx", "partial")) or + // Any collision (en, in, ee, ii, nn) + (collisionSpeciesMatch(collision_name, name, "", "coll", "partial"))) { + collision_names.push_back(collision_name); + } + } + + } else { + throw BoutException("\ndiffusion_collisions_mode for {:s} must be either " + "multispecies or braginskii", + name); } - // Multispecies mode: all collisions and CX are included - } else if (diffusion_collisions_mode == "multispecies") { - for (const auto& collision : localstate["collision_frequencies"].getChildren()) { - - std::string collision_name = collision.second.name(); - - if (// Charge exchange - (collisionSpeciesMatch( - collision_name, name, "", "cx", "partial")) or - // Any collision (en, in, ee, ii, nn) - (collisionSpeciesMatch( - collision_name, name, "", "coll", "partial"))) { - collision_names.push_back(collision_name); - } + + if (collision_names.empty()) { + throw BoutException("\tNo collisions found for {:s} in neutral_mixed for " + "selected collisions mode", + name); } - - } else { - throw BoutException("\ndiffusion_collisions_mode for {:s} must be either multispecies or braginskii", name); - } - if (collision_names.empty()) { - throw BoutException("\tNo collisions found for {:s} in neutral_mixed for selected collisions mode", name); + // Write chosen collisions to log file + output_info.write("\t{:s} neutral collisionality mode: '{:s}' using ", name, + diffusion_collisions_mode); + for (const auto& collision : collision_names) { + output_info.write("{:s} ", collision); + } + output_info.write("\n"); } - // Write chosen collisions to log file - output_info.write("\t{:s} neutral collisionality mode: '{:s}' using ", - name, diffusion_collisions_mode); - for (const auto& collision : collision_names) { - output_info.write("{:s} ", collision); - } - output_info.write("\n"); + // Collect the collisionalities based on list of names + nu = 0; + for (const auto& collision_name : collision_names) { + nu += GET_VALUE(Field3D, localstate["collision_frequencies"][collision_name]); } - // Collect the collisionalities based on list of names - nu = 0; - for (const auto& collision_name : collision_names) { - nu += GET_VALUE(Field3D, localstate["collision_frequencies"][collision_name]); + // Dnn = Vth^2 / sigma + Dnn = (Tnlim / AA) / (nu + Rnn); + } else { + Dnn = (Tnlim / AA) / Rnn; } - - - // Dnn = Vth^2 / sigma - Dnn = (Tnlim / AA) / (nu + Rnn); - } else { - Dnn = (Tnlim / AA) / Rnn; } - if (flux_limit > 0.0) { // Apply flux limit to diffusion, // using the local thermal speed and pressure gradient magnitude @@ -422,7 +463,7 @@ void NeutralMixed::finally(const Options& state) { } } - // Heat conductivity + // Heat conductivity // Note: This is kappa_n = (5/2) * Pn / (m * nu) // where nu is the collision frequency used in Dnn kappa_n = (5. / 2) * DnnNn; @@ -439,13 +480,11 @@ void NeutralMixed::finally(const Options& state) { ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - - ddt(Nn) = - - FV::Div_par_mod(Nn, Vn, sound_speed, pf_adv_par_ylow); // Parallel advection - - ddt(Nn) += Div_a_Grad_perp_flows(DnnNn, logPnlim, - pf_adv_perp_xlow, - pf_adv_perp_ylow); // Perpendicular advection + ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed, + pf_adv_par_ylow) // Parallel advection + + Div_a_Grad_perp_nonorthog(DnnNn, logPnlim, pf_adv_perp_xlow, + pf_adv_perp_ylow); // Perpendicular diffusion + ; Sn = density_source; // Save for possible output if (localstate.isSet("density_source")) { @@ -457,28 +496,30 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = - (5. / 3) * FV::Div_par_mod( // Parallel advection + ddt(Pn) = -(5. / 3) + * FV::Div_par_mod( // Parallel advection Pn, Vn, sound_speed, ef_adv_par_ylow) - + (2. / 3) * Vn * Grad_par(Pn) // Work done - + (5. / 3) * Div_a_Grad_perp_flows( // Perpendicular advection - DnnPn, logPnlim, - ef_adv_perp_xlow, ef_adv_perp_ylow) - ; + + (2. / 3) * Vn * Grad_par(Pn) // Work done + + (5. / 3) + * Div_a_Grad_perp_nonorthog( // Perpendicular advection + DnnPn, logPnlim, ef_adv_perp_xlow, ef_adv_perp_ylow); // The factor here is 5/2 as we're advecting internal energy and pressure. ef_adv_par_ylow *= 5./2; - ef_adv_perp_xlow *= 5./2; + ef_adv_perp_xlow *= 5. / 2; ef_adv_perp_ylow *= 5./2; if (neutral_conduction) { - ddt(Pn) += (2. / 3) * Div_a_Grad_perp_flows( - kappa_n, Tn, // Perpendicular conduction - ef_cond_perp_xlow, ef_cond_perp_ylow) - - + (2. / 3) * Div_par_K_Grad_par_mod(kappa_n, Tn, // Parallel conduction - ef_cond_par_ylow, - false) // No conduction through target boundary - ; + ddt(Pn) += + (2. / 3) + * Div_a_Grad_perp_nonorthog(kappa_n, Tn, // Perpendicular conduction + ef_cond_perp_xlow, ef_cond_perp_ylow) + + + (2. / 3) + * Div_par_K_Grad_par_mod(kappa_n, Tn, // Parallel conduction + ef_cond_par_ylow, + false) // No conduction through target boundary + ; // The factor here is likely 3/2 as this is pure energy flow, but needs checking. ef_cond_perp_xlow *= 3./2; ef_cond_perp_ylow *= 3./2; @@ -497,48 +538,46 @@ void NeutralMixed::finally(const Options& state) { // Neutral momentum TRACE("Neutral momentum"); - ddt(NVn) = - -AA * FV::Div_par_fvv( // Momentum flow - Nnlim, Vn, sound_speed) + ddt(NVn) = -AA + * FV::Div_par_fvv( // Momentum flow + Nnlim, Vn, sound_speed) - - Grad_par(Pn) // Pressure gradient - - + Div_a_Grad_perp_flows(DnnNVn, logPnlim, - mf_adv_perp_xlow, - mf_adv_perp_ylow) // Perpendicular advection - ; + - Grad_par(Pn) // Pressure gradient + + + Div_a_Grad_perp_nonorthog(DnnNVn, logPnlim, mf_adv_perp_xlow, + mf_adv_perp_ylow) // Perpendicular diffusion + // + Div_a_Grad_perp_flows(DnnNVn, logPnlim, + // mf_adv_perp_xlow, + // mf_adv_perp_ylow) // Perpendicular advection + ; if (neutral_viscosity) { // NOTE: The following viscosity terms are not (yet) balanced // by a viscous heating term - // Relationship between heat conduction and viscosity for neutral // gas Chapman, Cowling "The Mathematical Theory of Non-Uniform // Gases", CUP 1952 Ferziger, Kaper "Mathematical Theory of // Transport Processes in Gases", 1972 // eta_n = (2. / 5) * kappa_n; - Field3D viscosity_source = Div_a_Grad_perp_flows( - eta_n, Vn, // Perpendicular viscosity - mf_visc_perp_xlow, - mf_visc_perp_ylow) - - + Div_par_K_Grad_par_mod( // Parallel viscosity - eta_n, Vn, - mf_visc_par_ylow, - false) // No viscosity through target boundary - ; + Field3D viscosity_source = + Div_a_Grad_perp_nonorthog(eta_n, Vn, // Perpendicular viscosity + mf_visc_perp_xlow, + mf_visc_perp_ylow) + + + Div_par_K_Grad_par_mod( // Parallel viscosity + eta_n, Vn, mf_visc_par_ylow, + false) // No viscosity through target boundary + ; ddt(NVn) += viscosity_source; ddt(Pn) += -(2. /3) * Vn * viscosity_source; } - + Snv = momentum_source; if (localstate.isSet("momentum_source")) { - Snv = get(localstate["momentum_source"]); - ddt(NVn) += Snv; - } else { - Snv = 0; + Snv += get(localstate["momentum_source"]); } + ddt(NVn) += Snv; } else { ddt(NVn) = 0; diff --git a/tests/integrated/drift-wave/analysis.py b/tests/integrated/drift-wave/analysis.py index 2fefe281b..30afdfd23 100644 --- a/tests/integrated/drift-wave/analysis.py +++ b/tests/integrated/drift-wave/analysis.py @@ -89,6 +89,9 @@ def extract_data(path): dt = float(bd["t"][1] - bd["t"][0]) # Seconds gamma = np.mean(np.gradient(np.log(nrms))[-10:]) / dt # Growth rate [1/s] + # close dataset after last use + bd.close() + print(f"gamma = {gamma}, gamma / w* = {gamma / wstar}") # Track the motion of the peak to infer phase velocity diff --git a/tests/mms_operator/BOUT.inp.neutral_mixed.template b/tests/mms_operator/BOUT.inp.neutral_mixed.template new file mode 100644 index 000000000..ab767eed4 --- /dev/null +++ b/tests/mms_operator/BOUT.inp.neutral_mixed.template @@ -0,0 +1,46 @@ +nout = 0 +timestep = 1.0 +MYG=1 +MXG=2 +NXPE=1 + +[solver] +type = pvode +use_precon = true +mxstep = 1e5 + +atol = 1e-12 +rtol = 1e-10 + +[hermes] +components = (d+, d, e) + +Nnorm = 1 # Reference density [m^-3] +Bnorm = 1 # Reference magnetic field [T] +Tnorm = 1 # Reference temperature [eV] +recalculate_metric = false +normalise_metric = false +################################################################ +# Ions + +[d+] +type = fixed_density, fixed_velocity, fixed_temperature + +AA = 2 +charge = 1 +density = 1.0 +velocity = 0.0 +temperature = 1.0 +diagnose = true + +################################################################# +## Electrons +[e] +type = fixed_density, fixed_velocity, fixed_temperature +AA = 1/1836 +charge = -1 +density = 1.0 +velocity = 0.0 +temperature = 1.0 +diagnose = true + diff --git a/tests/mms_operator/BOUT.inp.template b/tests/mms_operator/BOUT.inp.template new file mode 100644 index 000000000..8a8d80791 --- /dev/null +++ b/tests/mms_operator/BOUT.inp.template @@ -0,0 +1,34 @@ +TwistShift = false # use twist-shift condition? +ballooning = false # ballooning transformation? +shiftinitial = false +ShiftWithoutTwist=false + +MYG=1 +MXG=2 +NXPE=1 + +dump_format = nc + +[input] +error_on_unused_options=false + +################################################## +# derivative methods + +[ddx] + +first = C2 +second = C2 +upwind = W3 + +[ddy] + +first = C2 +second = C2 +upwind = W3 + +[ddz] + +first = C2 +second = C2 +upwind = W3 \ No newline at end of file diff --git a/tests/mms_operator/README.md b/tests/mms_operator/README.md new file mode 100644 index 000000000..58c0acec3 --- /dev/null +++ b/tests/mms_operator/README.md @@ -0,0 +1,29 @@ +MMS tests of differential operators used in Hermes-3 +==================================================== + +These tests are designed to test that for any given +differential operator of two arguments `L(a,f)` +returns the expected result for a known `a`, `f`, and +contravariant metric coefficients, i.e., we compute +a symbolic result + + S = L(a,f) + +from prescribed symbolic expressions using `sympy`, +and we check that the numerical result for `S` is equal to the +symbolic one up to the numerical error. + +We carry out two tests, one on operators which are designed +for "orthogonal" metrics (`orthogonal_test.py`), and another for operators which +are generalised to fully nonorthogonal metrics (`nonorthogonal_test.py`). + +For each operator, we plot the numerical error, showing +how it declines compared to the expected convergence order +of 2 for the conservative finite difference methods. +These plots can be made interactively, or are otherwise plotted +to `fig_{i}.png` for the `ith` operator. + +The tests can be extended to include new operators by including +new operators in `const auto differential_operators` in `main.cxx` +and in `"differential_operator_list":` from the `test_input` dictionary +in the python test scripts `orthogonal_test.py` or `nonorthogonal_test.py`. diff --git a/tests/mms_operator/job_functions.py b/tests/mms_operator/job_functions.py new file mode 100644 index 000000000..5155ff631 --- /dev/null +++ b/tests/mms_operator/job_functions.py @@ -0,0 +1,535 @@ +# module to provide a test function for running the MMS tests of differential operators +import os +from xbout import open_boutdataset +import numpy as np +from scipy.optimize import curve_fit +from boututils.run_wrapper import launch_safe + +def lin_func(x,b,a): + return b*x + a + +def collectvar(datasets, var, mesh=0): + return datasets[mesh][var] + +def volume_integral(var,dx,dy,dz,J): + intV = np.sum(var[-1,:,:,0]*dx*dy*dz*J) + return intV + +def run_manufactured_solutions_test(test_input): + # expand inputs from input dictionary + # the number of grids generated + ntest = test_input["ntest"] + # the minimum number of points in each of the x, y, z grids in the test + # the number of points in the ith test is ngrid*i + nnbase = test_input["ngrid"] + # list of [name, symbolic string, expected convergence order] + differential_operator_test_list = test_input["differential_operator_list"] + astr = test_input["a_string"] + fstr = test_input["f_string"] + g11_str = test_input["g11_string"] + g22_str = test_input["g22_string"] + g33_str = test_input["g33_string"] + g12_str = test_input["g12_string"] + g13_str = test_input["g13_string"] + g23_str = test_input["g23_string"] + base_test_dir = test_input["test_dir"] + interactive_plots = test_input["interactive_plots"] + + # create directory + if not os.path.isdir(base_test_dir): + os.system("mkdir "+base_test_dir) + workdirs = [] + # make test for each resolution based on template file + for i in range(0,ntest): + workdir = f"{base_test_dir}/slab-mms-test-{i}" + workdirs.append(workdir) + # create directory + if not os.path.isdir(workdir): + os.system("mkdir "+workdir) + # copy template + file = workdir+"/BOUT.inp" + os.system(f"cp BOUT.inp.template "+file) + # update with mesh values for test + nn = nnbase*(i+1) # number of points in each grid + dd = 2.0*np.pi/nn # y z grid spacing, z, y on [0, 2pi] + ddx = 1.0/(nn-4) # x grid spacing to account for guard cells, x on [0,1] + # symmetricGlobalX = true ensures that x = 0 and x = 1 sits + # halfway between the last true grid point and the first guard point. + with open(file,"a") as file: + mesh_string = f""" + [mesh] + symmetricGlobalX = true + extrapolate_y = false + extrapolate_x= false + extrapolate_z= false + + nx = {nn} + dx = {ddx} + ny = {nn} + dy = {dd} + nz = {nn} + dz = {dd} + + g11 = {g11_str} + g22 = {g22_str} + g33 = {g33_str} + g12 = {g12_str} + g23 = {g23_str} + g13 = {g13_str} + x_input = x + y_input = y + z_input = z + a = {astr} + f = {fstr} + """ + n_operators = len(differential_operator_test_list) + mesh_string += f""" + # information about differential operators + n_operators = {n_operators} + """ + for i in range(0,n_operators): + mesh_string += f""" + differential_operator_name_{i} = {differential_operator_test_list[i][0]} + expected_result_{i} = {differential_operator_test_list[i][1]} + """ + file.write(mesh_string.replace("**","^")) + + # run job on this input + cmd = "../.././hermes_mms_operator_tests -d "+workdir+" > "+workdir+"/output.txt" + print(cmd) + # Launch using MPI + s, out = launch_safe(cmd, nproc=1, mthread=1, pipe=True) + if s != 0: + print(f"Command exited with status {s}. STDOUT printed below:") + with open(workdir + "/output.txt") as f: + print(f.read()) + + # now analyse the results of the test + # this slice avoids including guard cells in the test + # need guard cells in x (assume 2 here) and guard cells in y (assume 1) + # no guard cells in z + s = slice(2, -2), slice(1, -1), slice(None) + + # a dictionary of plot data, filled later on + plot_data = dict() + + # open the series of "workdir/BOUT.0.nc" files, + # saving them in a list `datasets` + datasets = [] + for workdir in workdirs: + boutmeshpath = workdir+"/"+f'BOUT.0.nc' + boutinppath = workdir+"/"+'BOUT.inp' + datasets.append(open_boutdataset(boutmeshpath, inputfilepath=boutinppath, keep_yboundaries=False)) + + # make a easy scan over the two operators, generalisation to N operators possible + for i, values in enumerate(differential_operator_test_list): + label = values[0] + expected_slope = values[2] + l2norm = [] + nylist = [] + dylist = [] + for m in range(0,ntest): + numerical = collectvar(datasets, f"result_{i}", m) + expected = collectvar(datasets, f"expected_result_{i}", m) + xx = collectvar(datasets, "x_input", m) + yy = collectvar(datasets, "y_input", m) + zz = collectvar(datasets, "z_input", m) + ff = collectvar(datasets, "f", m) + aa = collectvar(datasets, "a", m) + + error_values = (numerical - expected)[s] + thisl2 = np.sqrt(np.mean(error_values**2)) + l2norm.append(thisl2) + nylist.append(yy.shape[1]) + # proxy for grid spacing + dylist.append(1.0/yy.shape[1]) + + # cast lists as numpy arrays for further manipulation + l2norm = np.array(l2norm) + #print("test error: ",l2norm) + nylist = np.array(nylist) + dylist = np.array(dylist) + + # find linear fit coefficients to test convergence rate + # and construct fit function for plotting + try: + logl2norm = np.log(l2norm) + logdylist = np.log(dylist) + outvals = curve_fit(lin_func,logdylist,logl2norm) + coeffs = outvals[0] + slope = coeffs[0] # the convergence order + offset = coeffs[1] + logfit = slope*logdylist + offset + fitfunc = np.exp(logfit) + except ValueError: + print("Infs/Nans encountered in fit, skipping") + slope = None + offset = None + fitfunc = None + + # record results in dictionary and plot + #label = attrs["operator"] + " : f = " + attrs["inp"] + #label = "FV::Div_a_Grad_perp(a, f)" + plot_data[label] = [dylist, l2norm, fitfunc, slope, offset, expected_slope] + + # close the datasets + for dataset in datasets: + dataset.close() + + # plot the results + try: + import matplotlib.pyplot as plt + ifig = 0 + for key, variable_set in plot_data.items(): + (xaxis, yaxis, fit, slope, offset, expected_slope) = variable_set + plt.figure() + plt.plot(xaxis, yaxis, "x-", label="$\\epsilon(\\mathcal{L}\\ast f)$: "+key) + plt.plot(xaxis, yaxis[0]*(xaxis/xaxis[0])**expected_slope, "x-", label="$\\propto \\Delta^{{{:.2f}}}$".format(expected_slope)) + if not fit is None: + plt.plot(xaxis, fit, "x-", label="$e^{{{:.2f}}}\\Delta^{{{:.2f}}}$".format(offset,slope)) + plt.xlabel("$\\Delta = 1/N_y$") + plt.title(key) + plt.legend() + if not fit is None: + plt.gca().set_yscale("log") + plt.gca().set_xscale("log") + else: + print("l2 error: ",yaxis) + plt.savefig(f"{base_test_dir}/fig_{ifig}.png") + if interactive_plots: + plt.show() + plt.close() + ifig+=1 + except: + # Plotting could fail for any number of reasons, and the actual + # error raised may depend on, among other things, the current + # matplotlib backend, so catch everything + pass + + # test the convergence rates + success = True + output_message = "" + for key, variable_set in plot_data.items(): + this_test_success = True + (xaxis, yaxis, fit, slope, offset, expected_slope) = variable_set + # check slope of fit ~= 2 + slope_min = 0.975*expected_slope + if not slope is None: + if slope < slope_min: + this_test_success = False + else: # or permit near-zero errors, but nothing larger + for error in yaxis: + if error > 1.0e-10: + this_test_success = False + # append test message and set global success variable + if this_test_success: + output_message += f"{key} convergence order {slope:.2f} > {slope_min:.2f} => Test passed \n" + else: + output_message += f"{key} convergence order {slope:.2f} < {slope_min:.2f} => Test failed \n" + success = False + + return success, output_message + +def run_neutral_mixed_manufactured_solutions_test(test_input): + # expand inputs from input dictionary + # the number of grids generated + ntest = test_input["ntest"] + # the minimum number of points in each of the x, y, z grids in the test + # the number of points in the ith test is ngrid*i + nnbase = test_input["ngrid"] + # list of [name, symbolic string, expected convergence order] + nu_cfreq = test_input["collision_frequency"] + Nd_string = test_input["Nd_string"] + Pd_string = test_input["Pd_string"] + NVd_string = test_input["NVd_string"] + source_Nd_string = test_input["source_Nd_string"] + source_Pd_string = test_input["source_Pd_string"] + source_NVd_string = test_input["source_NVd_string"] + g11_str = test_input["g11_string"] + g22_str = test_input["g22_string"] + g33_str = test_input["g33_string"] + g12_str = test_input["g12_string"] + g13_str = test_input["g13_string"] + g23_str = test_input["g23_string"] + g_11_str = test_input["g_11_string"] + g_22_str = test_input["g_22_string"] + g_33_str = test_input["g_33_string"] + g_12_str = test_input["g_12_string"] + g_13_str = test_input["g_13_string"] + g_23_str = test_input["g_23_string"] + J_str = test_input["J_string"] + mass = test_input["mass"] + neutral_conduction = test_input["neutral_conduction"] + neutral_viscosity = test_input["neutral_viscosity"] + evolve_momentum = test_input["evolve_momentum"] + base_test_dir = test_input["test_dir"] + sub_test_dir = test_input["sub_test_dir"] + interactive_plots = test_input["interactive_plots"] + conservation_test = test_input["conservation_test"] + expected_convergence_order = test_input["expected_convergence_order"] + + # create directory + if not os.path.isdir(base_test_dir): + os.system("mkdir "+base_test_dir) + # create sub-directory, if required + if not sub_test_dir is None: + base_test_dir = base_test_dir + "/" + sub_test_dir + if not os.path.isdir(base_test_dir): + os.system("mkdir "+base_test_dir) + + workdirs = [] + # make test for each resolution based on template file + for i in range(0,ntest): + workdir = f"{base_test_dir}/slab-mms-test-{i}" + workdirs.append(workdir) + # create directory + if not os.path.isdir(workdir): + os.system("mkdir "+workdir) + # copy template + file = workdir+"/BOUT.inp" + os.system(f"cp BOUT.inp.neutral_mixed.template "+file) + # update with mesh values for test + nn = nnbase*(i+1) # number of points in each grid + dd = 2.0*np.pi/nn # y z grid spacing, z, y on [0, 2pi] + ddx = 1.0/(nn-4) # x grid spacing to account for guard cells, x on [0,1] + # symmetricGlobalX = true ensures that x = 0 and x = 1 sits + # halfway between the last true grid point and the first guard point. + with open(file,"a") as file: + mesh_string = f""" + [mesh] + symmetricGlobalX = true + extrapolate_y = false + extrapolate_x= false + + nx = {nn} + dx = {ddx} + ny = {nn} + dy = {dd} + nz = {1} + dz = {1.0} + + g11 = {g11_str} + g22 = {g22_str} + g33 = {g33_str} + g12 = {g12_str} + g23 = {g23_str} + g13 = {g13_str} + g_11 = {g_11_str} + g_22 = {g_22_str} + g_33 = {g_33_str} + g_12 = {g_12_str} + g_23 = {g_23_str} + g_13 = {g_13_str} + J = {J_str} + + Nd_src = {source_Nd_string} + Pd_src = {source_Pd_string} + NVd_src = {source_NVd_string} + + ################################################################# + # Neutrals + + [d] + type = neutral_mixed + + AA = {mass} + + diagnose = true + output_ddt = true + evolve_momentum = {evolve_momentum} + precondition = true + diagnose = true + collisionality_override = {nu_cfreq} + flux_limit = -1.0 + diffusion_limit = -1.0 + lax_flux = false + density_floor = 1.0e-15 + neutral_conduction = {neutral_conduction} + neutral_viscosity = {neutral_viscosity} + normalise_sources = false + [Nd] + + function = {Nd_string} + bndry_core = neumann + bndry_all = neumann + + [Pd] + function = {Pd_string} + bndry_core = neumann + bndry_all = neumann + """ + if evolve_momentum: + mesh_string = mesh_string + f""" + [NVd] + function = {NVd_string} + bndry_core = neumann + bndry_all = neumann + """ + file.write(mesh_string.replace("**","^")) + + # run job on this input + cmd = "../.././hermes-3 -d "+workdir+" > "+workdir+"/output.txt" + print(cmd) + # Launch using MPI + s, out = launch_safe(cmd, nproc=1, mthread=1, pipe=True) + if s != 0: + print(f"Command exited with status {s}. STDOUT printed below:") + with open(workdir + "/output.txt") as f: + print(f.read()) + + # now analyse the results of the test + # this slice avoids including guard cells in the test + # need to select last time point + # need guard cells in x (assume 2 here) and guard cells in y (assume 1) + # no guard cells in z + s = slice(None), slice(2, -2), slice(1, -1), slice(None) + #s = slice(1, 2), slice(2, -2), slice(1, -1), slice(None) + sxyz = slice(2, -2), slice(1, -1), slice(None) + sxy = slice(2, -2), slice(1, -1) + sx = slice(2, -2) + sy = slice(1, -1) + sz = slice(None) + # a dictionary of plot data, filled later on + plot_data = dict() + + # open the series of "workdir/BOUT.0.nc" files, + # saving them in a list `datasets` + datasets = [] + for workdir in workdirs: + boutmeshpath = workdir+"/"+f'BOUT.dmp.0.nc' + boutinppath = workdir+"/"+'BOUT.inp' + datasets.append(open_boutdataset(boutmeshpath, inputfilepath=boutinppath, keep_yboundaries=False)) + # By default xarray loads NetCDF files lazily. For some reason + # this causes the CI to fail (and not even in the same way each + # time!). Therefore, we load it eagerly here. + datasets[-1].load() + #for dataset in datasets: + # keys = dataset.keys() + # for key in keys: + # print(key) + # make a easy scan over the two operators, generalisation to N operators possible + function_list = ["ddt(Nd)", "ddt(Pd)"] + if evolve_momentum and not conservation_test: + function_list.append("ddt(NVd)") + + for varstring in function_list: + expected_slope = expected_convergence_order + l2norm = [] + nylist = [] + dylist = [] + for m in range(0,ntest): + if varstring == "ddt(Pd)" and evolve_momentum and conservation_test: + NVd = collectvar(datasets, "NVd", m) + Nd = collectvar(datasets, "Nd", m) + Vd = NVd/(Nd*mass) + numerical = (3.0/2.0)*collectvar(datasets, "ddt(Pd)", m) + Vd*collectvar(datasets, "ddt(NVd)", m) + label = "ddt(Ed) = 3/2 ddt(Pd) + Vd * ddt(NVd)" + else: + label = varstring + numerical = collectvar(datasets, varstring, m) + normvar = (collectvar(datasets, varstring[4:-1], m))[s] + dx = (collectvar(datasets, "dx", m))[sxy] + dy = (collectvar(datasets, "dy", m))[sxy] + dz = (collectvar(datasets, "dz", m))[sxy] + J = (collectvar(datasets, "J", m))[sxy] + error_values = (numerical)[s] + + if conservation_test: + thisl2 = np.abs(volume_integral(error_values,dx,dy,dz,J)/ + volume_integral(normvar,dx,dy,dz,J)) + else: + thisl2 = np.sqrt(volume_integral(error_values**2,dx,dy,dz,J)/ + volume_integral(normvar**2,dx,dy,dz,J)) + #print("thisl2",thisl2) + l2norm.append(thisl2) + nylist.append(numerical.shape[1]) + # proxy for grid spacing + dylist.append(1.0/numerical.shape[1]) + + # cast lists as numpy arrays for further manipulation + print(varstring) + l2norm = np.array(l2norm) + print("test error: ",l2norm) + nylist = np.array(nylist) + dylist = np.array(dylist) + print("dylist: ",dylist) + + # find linear fit coefficients to test convergence rate + # and construct fit function for plotting + try: + logl2norm = np.log(l2norm) + logdylist = np.log(dylist) + outvals = curve_fit(lin_func,logdylist,logl2norm) + coeffs = outvals[0] + slope = coeffs[0] # the convergence order + offset = coeffs[1] + logfit = slope*logdylist + offset + fitfunc = np.exp(logfit) + except ValueError: + print("Infs/Nans encountered in fit, skipping") + slope = None + offset = None + fitfunc = None + + # record results in dictionary and plot + #label = attrs["operator"] + " : f = " + attrs["inp"] + #label = "FV::Div_a_Grad_perp(a, f)" + plot_data[label] = [dylist, l2norm, fitfunc, slope, offset, expected_slope] + + # close the datasets + for dataset in datasets: + dataset.close() + + # plot the results + try: + import matplotlib.pyplot as plt + ifig = 0 + for key, variable_set in plot_data.items(): + (xaxis, yaxis, fit, slope, offset, expected_slope) = variable_set + plt.figure() + plt.plot(xaxis, yaxis, "x-", label="$\\epsilon$ "+key) + plt.plot(xaxis, yaxis[0]*(xaxis/xaxis[0])**expected_slope, "x-", label="$\\propto \\Delta^{{{:.2f}}}$".format(expected_slope)) + if not fit is None: + plt.plot(xaxis, fit, "x-", label="$e^{{{:.2f}}}\\Delta^{{{:.2f}}}$".format(offset,slope)) + plt.xlabel("$\\Delta = 1/N_y$") + plt.title(key) + plt.legend() + if not fit is None: + plt.gca().set_yscale("log") + plt.gca().set_xscale("log") + else: + print("l2 error: ",yaxis) + plt.savefig(f"{base_test_dir}/fig_{ifig}.png") + if interactive_plots: + plt.show() + plt.close() + ifig+=1 + except: + # Plotting could fail for any number of reasons, and the actual + # error raised may depend on, among other things, the current + # matplotlib backend, so catch everything + pass + + # test the convergence rates + success = True + output_message = "" + for key, variable_set in plot_data.items(): + this_test_success = True + (xaxis, yaxis, fit, slope, offset, expected_slope) = variable_set + # check slope of fit ~= 2 + slope_min = 0.975*expected_slope + if not slope is None: + if slope < slope_min: + this_test_success = False + else: # or permit near-zero errors, but nothing larger + for error in yaxis: + if error > 1.0e-10: + this_test_success = False + # append test message and set global success variable + if this_test_success: + output_message += f"{base_test_dir}: {key} convergence order {slope:.2f} > {slope_min:.2f} => Test passed \n" + else: + output_message += f"{base_test_dir}: {key} convergence order {slope:.2f} < {slope_min:.2f} => Test failed \n" + success = False + + return success, output_message diff --git a/tests/mms_operator/main.cxx b/tests/mms_operator/main.cxx new file mode 100644 index 000000000..eab40edba --- /dev/null +++ b/tests/mms_operator/main.cxx @@ -0,0 +1,264 @@ +#include "bout/field_factory.hxx" + +// #include "hermes-2.hxx" +#include "../include/div_ops.hxx" +#include "div_ops.hxx" +#include "bout/difops.hxx" +#include "bout/fv_ops.hxx" +#include "bout/version.hxx" +using ParLimiter = FV::Upwind; + +// Class used to store function of one argument +// in operator list below +class nameandfunction1 { +public: + std::string name; + std::function func; +}; + +// Class used to store function of two arguments +// in operator list below +class nameandfunction2 { +public: + std::string name; + std::function func; +}; + +class nameandfunction3 { +public: + std::string name; + std::function func; +}; + +class nameandfunction3_mod { +public: + std::string name; + std::function func; +}; + +class nameandfunction3_mod_diagnose { +public: + std::string name; + std::function func; +}; + +// Class used to store function of four arguments +// in operator list below. Second pair of arguments are +// internal diagnostics to the function +class nameandfunction4 { +public: + std::string name; + std::function func; +}; + +// List of tested operators of 1 arguments +const auto differential_operators_1_arg = { + nameandfunction1{"Div_par(f)", [](const Field3D& f) { return Div_par(f); }}, + nameandfunction1{"Grad_par(f)", [](const Field3D& f) { return Grad_par(f); }}, +}; +// List of tested operators of 2 arguments +const auto differential_operators_2_arg = { + nameandfunction2{ + "FV::Div_a_Grad_perp(a, f)", + [](const Field3D& a, const Field3D& f) { return FV::Div_a_Grad_perp(a, f); }}, +}; +// List of tested operators of 3 arguments +const auto differential_operators_3_arg = { + nameandfunction3{"Div_par_K_Grad_par_mod(a, f)", + [](const Field3D& a, const Field3D& f, Field3D& flow_ylow) { + return Div_par_K_Grad_par_mod(a, f, flow_ylow, true); + }}, +}; +// List of tested operators of 3 arguments +const auto differential_operators_3_arg_mod = { + nameandfunction3_mod{ + "FV::Div_par_fvv(f, v, wave_speed)", + [](const Field3D& f, const Field3D& v, const Field3D& wave_speed) { + return FV::Div_par_fvv(f, v, wave_speed, true); + }}, +}; +// List of tested operators of 3 arguments with a diagnostic +const auto differential_operators_3_arg_mod_diagnose = { + nameandfunction3_mod_diagnose{"FV::Div_par_mod(f, v, wave_speed)", + [](const Field3D& f, const Field3D& v, + const Field3D& wave_speed, Field3D& flow_ylow) { + return FV::Div_par_mod(f, v, wave_speed, + flow_ylow, true); + }}, +}; +// List of tested operators of 4 arguments +const auto differential_operators_4_arg = { + nameandfunction4{ + "Div_a_Grad_perp_nonorthog(a, f)", + [](const Field3D& a, const Field3D& f, Field3D& flow_xlow, Field3D& flow_ylow) { + return Div_a_Grad_perp_nonorthog(a, f, flow_xlow, flow_ylow); + }}, + nameandfunction4{ + "Div_a_Grad_perp_flows(a, f)", + [](const Field3D& a, const Field3D& f, Field3D& flow_xlow, Field3D& flow_ylow) { + return Div_a_Grad_perp_flows(a, f, flow_xlow, flow_ylow); + }}, +}; + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + int n_operators = Options::root()["mesh"]["n_operators"].withDefault(1); + + Mesh* mesh = Mesh::create(&Options::root()["mesh"]); + mesh->load(); + + Options dump; + + // Note that sticking all the inputs in the [mesh] section is a bit of a hack, but makes + // for a more compact input file than using the 'standard' variable-initialisation + // routines, which would require a separate section for each variable. + + // Load and save coordinate variables + Field3D xx{mesh}, yy{mesh}, zz{mesh}; + mesh->get(xx, "x_input", 0.0, false); + mesh->get(yy, "y_input", 0.0, false); + mesh->get(zz, "z_input", 0.0, false); + dump["x_input"] = xx; + dump["y_input"] = yy; + dump["z_input"] = zz; + + // Get coefficient from input file + Field3D a{mesh}; + mesh->get(a, "a", 1.0, false); + mesh->communicate(a); + dump["a"] = a; + + // Get test variable from input file + Field3D f{mesh}; + mesh->get(f, "f", 0.0, false); + mesh->communicate(f); + dump["f"] = f; + + // diagnostic variables + Field3D flow_xlow{mesh}; + Field3D flow_ylow{mesh}; + // auxiliary variables + Field3D ones{1.0, mesh}; + Field3D zeros{0.0, mesh}; + + for (int i = 0; i < n_operators; i++) { + std::string inputname = "differential_operator_name_" + std::to_string(i); + std::string expectedname = "expected_result_" + std::to_string(i); + std::string outname = "result_" + std::to_string(i); + std::string outname_flow_xlow = "result_flow_xlow_" + std::to_string(i); + std::string outname_flow_ylow = "result_flow_ylow_" + std::to_string(i); + std::string differential_operator_name = + Options::root()["mesh"][inputname].withDefault("FV::Div_a_Grad_perp(a, f)"); + // the for loop and if statement below should be replaced + // by a neater indexing syntax below if possible + for (const auto& difop : differential_operators_1_arg) { + if (difop.name.compare(differential_operator_name) == 0) { + // Get result of applying the named differential operator + Field3D result = difop.func(f); + dump[outname] = result; + dump[outname].setAttributes({ + {"operator", difop.name}, + }); + // Get expected result from input file + Field3D expected_result{mesh}; + mesh->get(expected_result, expectedname, 0.0, false); + dump[expectedname] = expected_result; + } + } + for (const auto& difop : differential_operators_2_arg) { + if (difop.name.compare(differential_operator_name) == 0) { + // Get result of applying the named differential operator + Field3D result = difop.func(a, f); + dump[outname] = result; + dump[outname].setAttributes({ + {"operator", difop.name}, + }); + // Get expected result from input file + Field3D expected_result{mesh}; + mesh->get(expected_result, expectedname, 0.0, false); + dump[expectedname] = expected_result; + } + } + for (const auto& difop : differential_operators_3_arg) { + if (difop.name.compare(differential_operator_name) == 0) { + // Get result of applying the named differential operator + Field3D result = difop.func(a, f, flow_ylow); + dump[outname] = result; + dump[outname].setAttributes({ + {"operator", difop.name}, + }); + // Get expected result from input file + Field3D expected_result{mesh}; + mesh->get(expected_result, expectedname, 0.0, false); + dump[expectedname] = expected_result; + // dump diagnostics + dump[outname_flow_ylow] = flow_ylow; + } + } + for (const auto& difop : differential_operators_3_arg_mod) { + if (difop.name.compare(differential_operator_name) == 0) { + // Get result of applying the named differential operator + Field3D result = difop.func(f, ones, zeros); + dump[outname] = result; + dump[outname].setAttributes({ + {"operator", difop.name}, + }); + // Get expected result from input file + Field3D expected_result{mesh}; + mesh->get(expected_result, expectedname, 0.0, false); + dump[expectedname] = expected_result; + } + } + for (const auto& difop : differential_operators_3_arg_mod_diagnose) { + if (difop.name.compare(differential_operator_name) == 0) { + // Get result of applying the named differential operator + Field3D result = difop.func(f, ones, zeros, flow_ylow); + dump[outname] = result; + dump[outname].setAttributes({ + {"operator", difop.name}, + }); + // Get expected result from input file + Field3D expected_result{mesh}; + mesh->get(expected_result, expectedname, 0.0, false); + dump[expectedname] = expected_result; + // dump diagnostics + dump[outname_flow_ylow] = flow_ylow; + } + } + for (const auto& difop : differential_operators_4_arg) { + if (difop.name.compare(differential_operator_name) == 0) { + // Get result of applying the named differential operator + Field3D result = difop.func(a, f, flow_xlow, flow_ylow); + dump[outname] = result; + dump[outname].setAttributes({ + {"operator", difop.name}, + }); + // Get expected result from input file + Field3D expected_result{mesh}; + mesh->get(expected_result, expectedname, 0.0, false); + dump[expectedname] = expected_result; + // dump diagnostics + dump[outname_flow_xlow] = flow_xlow; + dump[outname_flow_ylow] = flow_ylow; + } + } + } + // Field3D result = FV::Div_a_Grad_perp(a, f); + // dump["result"] = result; + + // Field3D result_nonorthog = Div_a_Grad_perp_nonorthog(a, f); + // dump["result_nonorthog"] = result_nonorthog; + + mesh->outputVars(dump); + + std::string outname = fmt::format( + "{}/BOUT.{}.nc", Options::root()["datadir"].withDefault("data"), + BoutComm::rank()); + + bout::OptionsIO::create(outname)->write(dump); + + BoutFinalise(); + + return 0; +} diff --git a/tests/mms_operator/neutral_mixed_ddt_test.py b/tests/mms_operator/neutral_mixed_ddt_test.py new file mode 100644 index 000000000..cdc993a53 --- /dev/null +++ b/tests/mms_operator/neutral_mixed_ddt_test.py @@ -0,0 +1,196 @@ +#!/usr/bin/env python3 +from job_functions import run_neutral_mixed_manufactured_solutions_test +from boutdata.mms import x, y, z, Div_par, Grad_par +from perpendicular_laplacian import Div_a_Grad_perp_f, Div_par_k_Grad_par_f, GeneralMetric +from sympy import sin, cos, log + +# specify symbolic inputs +# contravariant metric coeffs +# g11 = 1.1 + 0.16*x*cos(y) +# g22 = 0.9 + 0.09*x*cos(y) +# g33 = 1.2 + 0.2*x*cos(y) +# g12 = 0.0 +# g23 = 0.5 + 0.15*x*cos(y) +# g13 = 0.0 +# we can safely specify the simplest metric +# here knowing that the symbolic operators +# are used in orthogonal_test.py and nonorthogonal_test.py +# to test more complex metrics. +g11 = 1.0 +g22 = 1.0 +g33 = 1.0 +g12 = 0.0 +g23 = 0.0 +g13 = 0.0 +# the metric object +metric = GeneralMetric( + g11=g11, g12=g12, g13=g13, + g22=g22, g23=g23, g33=g33) +# extract covariant metric coeffs and Jacobian +g_11 = metric.g_11 +g_12 = metric.g_12 +g_13 = metric.g_13 +g_22 = metric.g_22 +g_23 = metric.g_23 +g_33 = metric.g_33 +J = metric.J + +def neutral_mixed_equations(evolve_momentum=False, + conservation_test=False, + # expected N in discretisation error ~ (grid_spacing)^N + expected_convergence_order=2.0): + pfac = 0.2 + nfac = 0.2 + # mass + AA = 2.0 + # collision frequency + nu_cfreq = 1.0 + # Manufactured solutions + Nn = 0.5 + nfac*(0.5 + x**3 - 1.5*x**2)*sin(y) + Pn = 1.0 + pfac*(0.5 + x**3 - 1.5*x**2)*sin(2*y) + if evolve_momentum: + # choose function that is both zero on boundary in x + # and has no derivative + NVn = (x**2 - 2*x**3 + x**4)*sin(y) + #NVn = 0.0 + x*(x-1)*cos(y) + else: + NVn = 0.0 + + # Derived quantities + Vn = NVn/(Nn*AA) + Tn = Pn/Nn + logPn = log(Pn) + Dn = (Tn/AA)/nu_cfreq + Kn = (5.0/2.0)*Nn*Dn + Etan = (2.0/5.0)*AA*Kn + neutral_conduction = True + neutral_viscosity = True + + # time evolution equations without sources + # N.B. + # Div . ( a Grad_perp f ) + #div_a_grad_perp_f = div_a_grad_perp_f_symbolic(g11, g12, g13, g22, g23, g33, a, f) + # Div . ( a Grad_par f) + #div_par_k_grad_par_f = div_par_k_grad_par_f_symbolic(g11, g12, g13, g22, g23, g33, a, f) + # Div . ( vec(b) f) + #div_par_f = div_par_f_symbolic(g11, g12, g13, g22, g23, g33, f) + # vec(b) . Grad f + #grad_par_f = grad_par_f_symbolic(g11, g12, g13, g22, g23, g33, f) + + ddt_Nn = ( -Div_par(Nn*Vn, metric=metric) + -Div_a_Grad_perp_f(-Dn*Nn, logPn, metric=metric) + ) + + ddt_Pn = ( -Div_par(Pn*Vn, metric=metric) + -(5.0/3.0)*Div_a_Grad_perp_f(-Dn*Pn, logPn, metric=metric) + -(2.0/3.0)*Pn*Div_par(Vn, metric=metric) + ) + + ddt_NVn = ( - AA*Div_par(Nn*Vn*Vn, metric=metric) + + AA*Div_a_Grad_perp_f(Nn*Vn*Dn, logPn, metric=metric) + - Grad_par(Pn, metric=metric) + ) + # if neutral conduction included in test + if neutral_conduction: + ddt_Pn = (ddt_Pn + + +(2.0/3.0)*Div_a_Grad_perp_f(Kn, Tn, metric=metric) + +(2.0/3.0)*Div_par_k_Grad_par_f(Kn, Tn, metric=metric) + ) + # if neutral viscosity included in test + if neutral_viscosity: + viscosity_source = (Div_a_Grad_perp_f(Etan, Vn, metric=metric) + + Div_par_k_Grad_par_f(Etan, Vn, metric=metric)) + ddt_NVn = (ddt_NVn + + viscosity_source + ) + ddt_Pn = (ddt_Pn - + (2.0/3.0)*Vn*viscosity_source + ) + + # Sources for steady state, set to exactly cancel + # the RHS operators if aiming to test the definition + # of the RHS, or set to zero if we just want to test + # that the RHS is conservative. + if conservation_test: + # we just want to make ddt and test conservation properties + # so no sources are required here + source_Nn = 0.0 + source_Pn = 0.0 + source_NVn = 0.0 + else: + # we want to cancel the numerically calculated ddt with an + # analytical source, so set the sources + source_Nn = -ddt_Nn + source_Pn = -ddt_Pn + source_NVn = -ddt_NVn + + test_input = { + "ntest" : 3, + "ngrid" : 10, + "collision_frequency" : nu_cfreq, + "g11_string": str(g11), + "g22_string": str(g22), + "g33_string": str(g33), + "g12_string": str(g12), + "g13_string": str(g13), + "g23_string": str(g23), + "g_11_string": str(g_11), + "g_22_string": str(g_22), + "g_33_string": str(g_33), + "g_12_string": str(g_12), + "g_13_string": str(g_13), + "g_23_string": str(g_23), + "J_string": str(J), + "mass" : AA, + "Nd_string" : str(Nn), + "Pd_string" : str(Pn), + "NVd_string" : str(NVn), + "source_Nd_string" : str(source_Nn), + "source_Pd_string" : str(source_Pn), + "source_NVd_string" : str(source_NVn), + "neutral_conduction" : neutral_conduction, + "neutral_viscosity" : neutral_viscosity, + "evolve_momentum" : evolve_momentum, + "test_dir" : "neutral_mixed", + "sub_test_dir" : "evolve_momentum_"+str(evolve_momentum)+"_conservation_test_"+str(conservation_test), + "interactive_plots" : False, + "conservation_test" : conservation_test, + "expected_convergence_order" : expected_convergence_order, + } + return test_input + +global_success = True + +test_input = neutral_mixed_equations(evolve_momentum=False, + conservation_test=False, + expected_convergence_order=2.0) +success, output_message_0 = run_neutral_mixed_manufactured_solutions_test(test_input) +global_success = success and global_success + +test_input = neutral_mixed_equations(evolve_momentum=False, + conservation_test=True, + expected_convergence_order=2.0) +success, output_message_1 = run_neutral_mixed_manufactured_solutions_test(test_input) +global_success = success and global_success + +test_input = neutral_mixed_equations(evolve_momentum=True, + conservation_test=False, + expected_convergence_order=1.0) +success, output_message_2 = run_neutral_mixed_manufactured_solutions_test(test_input) +global_success = success and global_success + +test_input = neutral_mixed_equations(evolve_momentum=True, + conservation_test=True, + expected_convergence_order=0.8) +success, output_message_3 = run_neutral_mixed_manufactured_solutions_test(test_input) +global_success = success and global_success + +print(output_message_0) +print(output_message_1) +print(output_message_2) +print(output_message_3) + +if global_success: + exit(0) +else: + exit(1) diff --git a/tests/mms_operator/nonorthogonal_test.py b/tests/mms_operator/nonorthogonal_test.py new file mode 100644 index 000000000..bbd521923 --- /dev/null +++ b/tests/mms_operator/nonorthogonal_test.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python3 +from job_functions import run_manufactured_solutions_test +from boutdata.mms import x, y, z, Div_par, Grad_par +from perpendicular_laplacian import Div_a_Grad_perp_f, Div_par_k_Grad_par_f, GeneralMetric +from sympy import sin, cos + +# specify symbolic inputs +# contravariant metric coeffs +g11 = 1.1 + 0.16*x*cos(y) +g22 = 0.9 + 0.09*x*cos(y) +g33 = 1.2 + 0.2*x*cos(y) +g12 = 0.2 + 0.11*x*cos(y) +g23 = 0.5 + 0.15*x*cos(y) +g13 = 0.3 + 0.13*x*cos(y) +# f and a +f = (x**2)*sin(y)*sin(z) +a = 1.0 + 0.1*x**3*sin(2*y)*sin(2*z) +# the metric object +metric = GeneralMetric( + g11=g11, g12=g12, g13=g13, + g22=g22, g23=g23, g33=g33) +# Div . ( a Grad_perp f ) +div_a_grad_perp_f = Div_a_Grad_perp_f(a, f, metric=metric) +# Div . ( a Grad_par f) +div_par_k_grad_par_f = Div_par_k_Grad_par_f(a, f, metric=metric) +# Div . ( vec(b) f) +div_par_f = Div_par(f, metric=metric) +# vec(b) . Grad f +grad_par_f = Grad_par(f, metric=metric) + +test_input = { + "ntest" : 3, + "ngrid" : 20, + # list of list of ["name", symbolic function, expected convergence order] + "differential_operator_list": [["Div_a_Grad_perp_nonorthog(a, f)", str(div_a_grad_perp_f), 2], + ["Div_par_K_Grad_par_mod(a, f)", str(div_par_k_grad_par_f), 2], + ["Div_par(f)", str(div_par_f), 2], + ["FV::Div_par_fvv(f, v, wave_speed)", str(div_par_f), 1], + ["FV::Div_par_mod(f, v, wave_speed)", str(div_par_f), 1], + ["Grad_par(f)", str(grad_par_f), 2], + ], + "a_string": str(a), + "f_string": str(f), + "g11_string": str(g11), + "g22_string": str(g22), + "g33_string": str(g33), + "g12_string": str(g12), + "g13_string": str(g13), + "g23_string": str(g23), + "test_dir" : "nonorthogonal", + "interactive_plots" : False +} + +success, output_message = run_manufactured_solutions_test(test_input) + +print(output_message) +if success: + exit(0) +else: + exit(1) diff --git a/tests/mms_operator/orthogonal_test.py b/tests/mms_operator/orthogonal_test.py new file mode 100644 index 000000000..714441db9 --- /dev/null +++ b/tests/mms_operator/orthogonal_test.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +from job_functions import run_manufactured_solutions_test +from boutdata.mms import x, y, z, Div_par, Grad_par +from perpendicular_laplacian import Div_a_Grad_perp_f, Div_par_k_Grad_par_f, GeneralMetric +from sympy import sin, cos + +# specify symbolic inputs +# contravariant metric coeffs +g11 = 1.1 + 0.16*x*cos(y) +g22 = 0.9 + 0.09*x*cos(y) +g33 = 1.2 + 0.2*x*cos(y) +g12 = 0.0 +g23 = 0.5 + 0.15*x*cos(y) +g13 = 0.0 +# f and a +f = (x**2)*sin(y)*sin(z) +a = 1.0 + 0.1*x**3*sin(2*y)*sin(2*z) +# the metric object +metric = GeneralMetric( + g11=g11, g12=g12, g13=g13, + g22=g22, g23=g23, g33=g33) +# Div . ( a Grad_perp f ) +div_a_grad_perp_f = Div_a_Grad_perp_f(a, f, metric=metric) +# Div . ( a Grad_par f) +div_par_k_grad_par_f = Div_par_k_Grad_par_f(a, f, metric=metric) +# Div . ( vec(b) f) +div_par_f = Div_par(f, metric=metric) +# vec(b) . Grad f +grad_par_f = Grad_par(f, metric=metric) + +test_input = { + "ntest" : 3, + "ngrid" : 20, + # list of list of ["name", symbolic function, expected convergence order] + "differential_operator_list": [["FV::Div_a_Grad_perp(a, f)", str(div_a_grad_perp_f), 2], + ["Div_a_Grad_perp_nonorthog(a, f)", str(div_a_grad_perp_f), 2], + ["Div_a_Grad_perp_flows(a, f)", str(div_a_grad_perp_f), 2], + ["Div_par_K_Grad_par_mod(a, f)", str(div_par_k_grad_par_f), 2], + ["Div_par(f)", str(div_par_f), 2], + ["FV::Div_par_fvv(f, v, wave_speed)", str(div_par_f), 1], + ["FV::Div_par_mod(f, v, wave_speed)", str(div_par_f), 1], + ["Grad_par(f)", str(grad_par_f), 2], + ], + "a_string": str(a), + "f_string": str(f), + "g11_string": str(g11), + "g22_string": str(g22), + "g33_string": str(g33), + "g12_string": str(g12), + "g13_string": str(g13), + "g23_string": str(g23), + "test_dir" : "orthogonal", + "interactive_plots" : False +} + +success, output_message = run_manufactured_solutions_test(test_input) + +print(output_message) +if success: + exit(0) +else: + exit(1) diff --git a/tests/mms_operator/perpendicular_laplacian.py b/tests/mms_operator/perpendicular_laplacian.py new file mode 100644 index 000000000..b2f028d91 --- /dev/null +++ b/tests/mms_operator/perpendicular_laplacian.py @@ -0,0 +1,86 @@ +from sympy import diff, sqrt +from sympy.matrices import Matrix + +from boutdata.mms import identity, x, y, z + +# A class very similar to boutdata.mms.Metric +# to permit construction of similar functions +# and use for Div_par and Grad_par from boutdata.mms +# A key difference is that the symbol variables +# x, y, z are the same as the global Cartesian variables x, y, z +# rather than the transformed x', y', z'. +class GeneralMetric: + def __init__(self, + g11=1.0, + g12=0.0, + g13=0.0, + g22=1.0, + g23=0.0, + g33=1.0): + # obtain g_{ij} from g^{ij} + gup = Matrix(3,3,[g11,g12,g13,g12,g22,g23,g13,g23,g33]) + gdown = gup.inv() + g_11 = gdown[0,0] + g_12 = gdown[0,1] + g_13 = gdown[0,2] + g_22 = gdown[1,1] + g_23 = gdown[1,2] + g_33 = gdown[2,2] + detgup = gup.det() + J = 1/sqrt(detgup) + #detgup = g11*g22*g33 - g11*g23*g23 - g12*g12*g33 + g12*g13*g23 - g13*g13*g22 + g13*g12*g23 + + self.g11 = g11 + self.g12 = g12 + self.g13 = g13 + self.g22 = g22 + self.g23 = g23 + self.g33 = g33 + self.g_11 = g_11 + self.g_12 = g_12 + self.g_13 = g_13 + self.g_22 = g_22 + self.g_23 = g_23 + self.g_33 = g_33 + self.J = J + self.B = sqrt(g_22)/J + # for the GeneralMetric, use the cartesian x, y, z + # global variables from boutdata as the internal variables + self.x = x + self.y = y + self.z = z + +# Div . ( a Grad_perp f ) +# see notes for formulae +# https://bout-dev.readthedocs.io/en/stable/user_docs/coordinates.html#the-perpendicular-laplacian-in-divergence-form +def Div_a_Grad_perp_f(a, f, metric=identity): + g11, g12, g13, g22, g23, g33 = metric.g11, metric.g12, metric.g13, metric.g22, metric.g23, metric.g33 + g_11, g_12, g_13, g_22, g_23, g_33 = metric.g_11, metric.g_12, metric.g_13, metric.g_22, metric.g_23, metric.g_33 + J = metric.J + # use the metric internal variables for derivatives + x, y, z = metric.x, metric.y, metric.z + + dfdx = diff(f,x) + dfdy = diff(f,y) + dfdz = diff(f,z) + + df1 = dfdx - (g_12/g_22)*dfdy + df3 = dfdz - (g_23/g_22)*dfdy + + a_grad_perp_f_x = a*J*(g11*df1 + g13*df3) + a_grad_perp_f_y = a*J*(g12*df1 + g23*df3) + a_grad_perp_f_z = a*J*(g13*df1 + g33*df3) + + div_a_grad_perp_f = (1/J)*(diff(a_grad_perp_f_x,x)+diff(a_grad_perp_f_y,y)+diff(a_grad_perp_f_z,z)) + return div_a_grad_perp_f + +def Div_par_k_Grad_par_f(a, f, metric=identity): + g11, g12, g13, g22, g23, g33 = metric.g11, metric.g12, metric.g13, metric.g22, metric.g23, metric.g33 + g_11, g_12, g_13, g_22, g_23, g_33 = metric.g_11, metric.g_12, metric.g_13, metric.g_22, metric.g_23, metric.g_33 + J = metric.J + # use the metric internal variables for derivatives + x, y, z = metric.x, metric.y, metric.z + dfdy = diff(f,y) + k_grad_par_f = a*J*dfdy/g_22 + div_par_k_grad_par_f = (1/J)*(diff(k_grad_par_f,y)) + return div_par_k_grad_par_f diff --git a/tests/unit/test_neutral_mixed.cxx b/tests/unit/test_neutral_mixed.cxx index fff532889..6c96521b1 100644 --- a/tests/unit/test_neutral_mixed.cxx +++ b/tests/unit/test_neutral_mixed.cxx @@ -19,36 +19,56 @@ using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using NeutralMixedTest = FakeMeshFixture; - +// Component test checking only that the state +// includes Nd, Pd, NVd as evolved variables. TEST_F(NeutralMixedTest, CreateComponent) { FakeSolver solver; Options options{ {"units", {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}, {"meters", 1.0}}}, - {"n", {{"AA", 2.0}}}}; - NeutralMixed component("n", options, &solver); + {"d", {{"type", "neutral_mixed"}, {"AA", 2.0}, {"evolve_momentum", true}}}}; + NeutralMixed component("d", options, &solver); Options state = solver.getState(); - EXPECT_TRUE(state.isSet("Nn")); - EXPECT_TRUE(state.isSet("Pn")); - EXPECT_TRUE(state.isSet("NVn")); + EXPECT_TRUE(state.isSet("Nd")); + EXPECT_TRUE(state.isSet("Pd")); + EXPECT_TRUE(state.isSet("NVd")); } +// Component test checking only that the state +// includes Nd, Pd as evolved variables when evolve_momentum = false. +TEST_F(NeutralMixedTest, CreateComponentEvolveMomentumFalse) { + FakeSolver solver; + Options options{ + {"units", + {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}, {"meters", 1.0}}}, + {"d", {{"type", "neutral_mixed"}, {"AA", 2.0}, {"evolve_momentum", false}}}}; + NeutralMixed component("d", options, &solver); + + Options state = solver.getState(); + + EXPECT_TRUE(state.isSet("Nd")); + EXPECT_TRUE(state.isSet("Pd")); + EXPECT_FALSE(state.isSet("NVd")); +} +// Transform test checking only that the state has data set in +// the the required auxiliary variables. Note that boundary conditions +// are also set in the evolved variables. TEST_F(NeutralMixedTest, Transform) { FakeSolver solver; Options options{ {"units", {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}, {"meters", 1.0}}}, - {"n", {{"AA", 2.0}}}}; - NeutralMixed component("n", options, &solver); + {"d", {{"type", "neutral_mixed"}, {"AA", 2.0}, {"evolve_momentum", true}}}}; + NeutralMixed component("d", options, &solver); Options state; component.transform(state); - auto& species = state["species"]["n"]; + Options& species = state["species"]["d"]; EXPECT_TRUE(species.isSet("density")); EXPECT_TRUE(species.isSet("AA")); EXPECT_TRUE(species.isSet("pressure")); @@ -56,3 +76,148 @@ TEST_F(NeutralMixedTest, Transform) { EXPECT_TRUE(species.isSet("velocity")); EXPECT_TRUE(species.isSet("temperature")); } +// Test of finally() for this component following +// tests of evolve_density and evolve_pressure. +// We provide a state vector where all fields are constants. +// We use finally to compute the ddt() of the evolved +// state variables. For these inputs, the only non-zero +// contribution comes from the sources, which may be tested +// by checking that the ddt() variables match expected constants. +// Note that this test does not check the definition of the +// differential operators in the right-hand-side of the equations. +TEST_F(NeutralMixedTest, Finally) { + FakeSolver solver; + + Options options{ + {"units", + {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}, {"meters", 1.0}}}, + {"d", {{"type", "neutral_mixed"}, {"AA", 2.0}, {"evolve_momentum", true}}}}; + NeutralMixed component("d", options, &solver); + + // Call the finally() method with a density, energy, and momentum source + const Options state = {{"species", + {{"d", + {{"density", 1.0}, + {"density_source", 0.5}, + {"pressure", 1.0}, + {"energy_source", 1.5}, + {"momentum", 1.0}, + {"momentum_source", 0.75}, + {"temperature", 1.0}, + {"velocity", 1.0}}}}}}; + component.finally(state); + + Options ddt = solver.getTimeDerivs(); + + EXPECT_TRUE(ddt.isSet("Nd")); + Field3D ddt_Nd = ddt["Nd"].as(); + + BOUT_FOR_SERIAL(i, ddt_Nd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(0.5, ddt_Nd[i]); + } + + EXPECT_TRUE(ddt.isSet("Pd")); + Field3D ddt_Pd = ddt["Pd"].as(); + + BOUT_FOR_SERIAL(i, ddt_Pd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(1.0, ddt_Pd[i]); + } + + EXPECT_TRUE(ddt.isSet("NVd")); + Field3D ddt_NVd = ddt["NVd"].as(); + + BOUT_FOR_SERIAL(i, ddt_NVd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(0.75, ddt_NVd[i]); + } +} +// Identical to the test above, but using the collisionality_override variable. +TEST_F(NeutralMixedTest, FinallyCollisionalityOverride) { + FakeSolver solver; + + Options options{ + {"units", + {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}, {"meters", 1.0}}}, + {"d", + {{"type", "neutral_mixed"}, + {"AA", 2.0}, + {"evolve_momentum", true}, + {"collisionality_override", 1.0}}}}; + NeutralMixed component("d", options, &solver); + + // Call the finally() method with a density, energy, and momentum source + const Options state = {{"species", + {{"d", + {{"density", 1.0}, + {"density_source", 0.5}, + {"pressure", 1.0}, + {"energy_source", 1.5}, + {"momentum", 1.0}, + {"momentum_source", 0.75}, + {"temperature", 1.0}, + {"velocity", 1.0}}}}}}; + component.finally(state); + + Options ddt = solver.getTimeDerivs(); + + EXPECT_TRUE(ddt.isSet("Nd")); + Field3D ddt_Nd = ddt["Nd"].as(); + + BOUT_FOR_SERIAL(i, ddt_Nd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(0.5, ddt_Nd[i]); + } + + EXPECT_TRUE(ddt.isSet("Pd")); + Field3D ddt_Pd = ddt["Pd"].as(); + + BOUT_FOR_SERIAL(i, ddt_Pd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(1.0, ddt_Pd[i]); + } + + EXPECT_TRUE(ddt.isSet("NVd")); + Field3D ddt_NVd = ddt["NVd"].as(); + + BOUT_FOR_SERIAL(i, ddt_NVd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(0.75, ddt_NVd[i]); + } +} +// Identical to the test above, but using evolve_momentum = false. +TEST_F(NeutralMixedTest, FinallyEvolveMomentumFalse) { + FakeSolver solver; + + Options options{ + {"units", + {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}, {"meters", 1.0}}}, + {"d", {{"type", "neutral_mixed"}, {"AA", 2.0}, {"evolve_momentum", false}}}}; + NeutralMixed component("d", options, &solver); + + // Call the finally() method with a density, energy, and momentum source + const Options state = {{"species", + {{"d", + {{"density", 1.0}, + {"density_source", 0.5}, + {"pressure", 1.0}, + {"energy_source", 1.5}, + {"momentum", 1.0}, + {"momentum_source", 0.75}, + {"temperature", 1.0}, + {"velocity", 1.0}}}}}}; + component.finally(state); + + Options ddt = solver.getTimeDerivs(); + + EXPECT_TRUE(ddt.isSet("Nd")); + Field3D ddt_Nd = ddt["Nd"].as(); + + BOUT_FOR_SERIAL(i, ddt_Nd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(0.5, ddt_Nd[i]); + } + + EXPECT_TRUE(ddt.isSet("Pd")); + Field3D ddt_Pd = ddt["Pd"].as(); + + BOUT_FOR_SERIAL(i, ddt_Pd.getRegion("RGN_NOBNDRY")) { + ASSERT_DOUBLE_EQ(1.0, ddt_Pd[i]); + } + + EXPECT_FALSE(ddt.isSet("NVd")); +}