diff --git a/external/BOUT-dev b/external/BOUT-dev index 7d28d67c3..5a18d876d 160000 --- a/external/BOUT-dev +++ b/external/BOUT-dev @@ -1 +1 @@ -Subproject commit 7d28d67c3f12c24ec281c0982e870f5369c65a6f +Subproject commit 5a18d876d3d69b5a02010f8be317e9e16f9d1c22 diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index 6efd3a8f3..0bde97a36 100644 --- a/src/evolve_density.cxx +++ b/src/evolve_density.cxx @@ -1,16 +1,17 @@ #include -#include -#include -#include #include #include +#include +#include #include +#include +#include #include "../include/div_ops.hxx" #include "../include/evolve_density.hxx" -#include "../include/hermes_utils.hxx" #include "../include/hermes_build_config.hxx" +#include "../include/hermes_utils.hxx" using bout::globals::mesh; @@ -28,9 +29,11 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7); - BoutReal temperature_floor = options["temperature_floor"].doc("Low temperature scale for low_T_diffuse_perp") - .withDefault(0.1) / get(alloptions["units"]["eV"]); - + BoutReal temperature_floor = options["temperature_floor"] + .doc("Low temperature scale for low_T_diffuse_perp") + .withDefault(0.1) + / get(alloptions["units"]["eV"]); + low_n_diffuse = options["low_n_diffuse"] .doc("Parallel diffusion at low density") .withDefault(false); @@ -84,15 +87,15 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv auto& n_options = alloptions[std::string("N") + name]; source_time_dependent = n_options["source_time_dependent"] - .doc("Use a time-dependent source?") - .withDefault(false); + .doc("Use a time-dependent source?") + .withDefault(false); source_only_in_core = n_options["source_only_in_core"] - .doc("Zero the source outside the closed field-line region?") - .withDefault(false); + .doc("Zero the source outside the closed field-line region?") + .withDefault(false); source_normalisation = Nnorm * Omega_ci; - time_normalisation = 1./Omega_ci; + time_normalisation = 1. / Omega_ci; // Try to read the density source from the mesh // Units of particles per cubic meter per second @@ -100,16 +103,17 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv mesh->get(source, std::string("N") + name + "_src"); // Allow the user to override the source from input file source = n_options["source"] - .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) - .withDefault(source) - / source_normalisation; + .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) + .withDefault(source) + / source_normalisation; // If time dependent, parse the function with respect to time from the input file if (source_time_dependent) { auto str = n_options["source_prefactor"] - .doc("Time-dependent function of multiplier on ddt(N" + name + std::string(") source.")) - .as(); - source_prefactor_function = FieldFactory::get()->parse(str, &n_options); + .doc("Time-dependent function of multiplier on ddt(N" + name + + std::string(") source.")) + .as(); + source_prefactor_function = FieldFactory::get()->parse(str, &n_options); } // Putting source at first X index would put it in both PFR in core, this ensures only core @@ -126,9 +130,10 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv } } - neumann_boundary_average_z = alloptions[std::string("N") + name]["neumann_boundary_average_z"] - .doc("Apply neumann boundary with Z average?") - .withDefault(false); + neumann_boundary_average_z = + alloptions[std::string("N") + name]["neumann_boundary_average_z"] + .doc("Apply neumann boundary with Z average?") + .withDefault(false); std::vector outputs = {"AA", "density", "density_source"}; if (charge != 0.) { @@ -213,7 +218,9 @@ void EvolveDensity::transform_impl(GuardedOptions& state) { if (source_time_dependent) { // Evaluate the source_prefactor function at the current time in seconds and scale source with it BoutReal time = get(state["time"]); - BoutReal source_prefactor = source_prefactor_function ->generate(bout::generator::Context().set("x",0,"y",0,"z",0,"t",time*time_normalisation)); + BoutReal source_prefactor = + source_prefactor_function->generate(bout::generator::Context().set( + "x", 0, "y", 0, "z", 0, "t", time * time_normalisation)); final_source = source * source_prefactor; } else { final_source = source; @@ -230,7 +237,8 @@ void EvolveDensity::finally(const Options& state) { // but retain densities which fall below zero N.setBoundaryTo(get(species["density"])); - if ((fabs(charge) > 1e-5) and state.isSection("fields") and state["fields"].isSet("phi")) { + if ((fabs(charge) > 1e-5) and state.isSection("fields") + and state["fields"].isSet("phi")) { // Electrostatic potential set and species is charged -> include ExB flow Field3D phi = get(state["fields"]["phi"]); @@ -373,25 +381,26 @@ void EvolveDensity::outputVars(Options& state) { auto rho_s0 = get(state["rho_s0"]); if (flow_xlow.isAllocated()) { - set_with_attrs(state[fmt::format("pf{}_tot_xlow", name)], flow_xlow, - {{"time_dimension", "t"}, - {"units", "s^-1"}, - {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, - {"standard_name", "particle flow"}, - {"long_name", name + " particle flow in X. Note: May be incomplete."}, - {"species", name}, - {"source", "evolve_density"}}); + set_with_attrs( + state[fmt::format("pf{}_tot_xlow", name)], flow_xlow, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, + {"standard_name", "particle flow"}, + {"long_name", name + " particle flow in X. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_density"}}); } if (flow_ylow.isAllocated()) { - set_with_attrs(state[fmt::format("pf{}_tot_ylow", name)], flow_ylow, - {{"time_dimension", "t"}, - {"units", "s^-1"}, - {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, - {"standard_name", "particle flow"}, - {"long_name", name + " particle flow in Y. Note: May be incomplete."}, - {"species", name}, - {"source", "evolve_density"}}); + set_with_attrs( + state[fmt::format("pf{}_tot_ylow", name)], flow_ylow, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, + {"standard_name", "particle flow"}, + {"long_name", name + " particle flow in Y. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_density"}}); } } } - diff --git a/src/evolve_energy.cxx b/src/evolve_energy.cxx index 6aef6bad8..7b182f361 100644 --- a/src/evolve_energy.cxx +++ b/src/evolve_energy.cxx @@ -6,11 +6,12 @@ #include #include #include +#include #include "../include/div_ops.hxx" #include "../include/evolve_energy.hxx" -#include "../include/hermes_utils.hxx" #include "../include/hermes_build_config.hxx" +#include "../include/hermes_utils.hxx" #include @@ -226,8 +227,8 @@ void EvolveEnergy::finally(const Options& state) { Field3D Pfloor = P; - if (species.isSet("charge") and (fabs(get(species["charge"])) > 1e-5) and - state.isSection("fields") and state["fields"].isSet("phi")) { + if (species.isSet("charge") and (fabs(get(species["charge"])) > 1e-5) + and state.isSection("fields") and state["fields"].isSet("phi")) { // Electrostatic potential set -> include ExB flow Field3D phi = get(state["fields"]["phi"]); @@ -279,7 +280,8 @@ void EvolveEnergy::finally(const Options& state) { } #if CHECKLEVEL >= 1 if (species.isSet("pressure_source")) { - throw BoutException("Components must evolve `energy_source` rather then `pressure_source`"); + throw BoutException( + "Components must evolve `energy_source` rather then `pressure_source`"); } #endif if (species.isSet("momentum_source")) { @@ -384,24 +386,26 @@ void EvolveEnergy::outputVars(Options& state) { {"source", "evolve_energy"}}); if (flow_xlow.isAllocated()) { - set_with_attrs(state[fmt::format("ef{}_tot_xlow", name)], flow_xlow, - {{"time_dimension", "t"}, - {"units", "W"}, - {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, - {"standard_name", "power"}, - {"long_name", name + " power through X cell face. Note: May be incomplete."}, - {"species", name}, - {"source", "evolve_energy"}}); + set_with_attrs( + state[fmt::format("ef{}_tot_xlow", name)], flow_xlow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " power through X cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_energy"}}); } if (flow_ylow.isAllocated()) { - set_with_attrs(state[fmt::format("ef{}_tot_ylow", name)], flow_ylow, - {{"time_dimension", "t"}, - {"units", "W"}, - {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, - {"standard_name", "power"}, - {"long_name", name + " power through Y cell face. Note: May be incomplete."}, - {"species", name}, - {"source", "evolve_energy"}}); + set_with_attrs( + state[fmt::format("ef{}_tot_ylow", name)], flow_ylow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " power through Y cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_energy"}}); } } } diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index 688632121..cacad8e4c 100644 --- a/src/evolve_momentum.cxx +++ b/src/evolve_momentum.cxx @@ -1,12 +1,13 @@ +#include #include #include -#include #include #include +#include -#include "../include/evolve_momentum.hxx" #include "../include/div_ops.hxx" +#include "../include/evolve_momentum.hxx" #include "../include/hermes_build_config.hxx" #include "../include/hermes_utils.hxx" @@ -25,8 +26,10 @@ EvolveMomentum::EvolveMomentum(std::string name, Options& alloptions, Solver* so density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7); - BoutReal temperature_floor = options["temperature_floor"].doc("Low temperature scale for low_T_diffuse_perp") - .withDefault(0.1) / get(alloptions["units"]["eV"]); + BoutReal temperature_floor = options["temperature_floor"] + .doc("Low temperature scale for low_T_diffuse_perp") + .withDefault(0.1) + / get(alloptions["units"]["eV"]); low_n_diffuse_perp = options["low_n_diffuse_perp"] .doc("Perpendicular diffusion at low density") @@ -39,24 +42,23 @@ EvolveMomentum::EvolveMomentum(std::string name, Options& alloptions, Solver* so .withDefault(false); bndry_flux = options["bndry_flux"] - .doc("Allow flows through radial boundaries") - .withDefault(true); + .doc("Allow flows through radial boundaries") + .withDefault(true); - poloidal_flows = options["poloidal_flows"] - .doc("Include poloidal ExB flow") - .withDefault(true); + poloidal_flows = + options["poloidal_flows"].doc("Include poloidal ExB flow").withDefault(true); hyper_z = options["hyper_z"].doc("Hyper-diffusion in Z").withDefault(-1.0); V.setBoundary(std::string("V") + name); - diagnose = options["diagnose"] - .doc("Output additional diagnostics?") - .withDefault(false); + diagnose = + options["diagnose"].doc("Output additional diagnostics?").withDefault(false); - fix_momentum_boundary_flux = options["fix_momentum_boundary_flux"] - .doc("Fix Y boundary momentum flux to boundary midpoint value?") - .withDefault(false); + fix_momentum_boundary_flux = + options["fix_momentum_boundary_flux"] + .doc("Fix Y boundary momentum flux to boundary midpoint value?") + .withDefault(false); // Set to zero so set for output momentum_source = 0.0; @@ -80,7 +82,7 @@ void EvolveMomentum::transform_impl(GuardedOptions& state) { V.applyBoundary(); set(species["velocity"], V); - NV_solver = NV; // Save the momentum as calculated by the solver + NV_solver = NV; // Save the momentum as calculated by the solver NV = AA * N * V; // Re-calculate consistent with V and N // Note: Now NV and NV_solver will differ when N < density_floor NV_err = NV - NV_solver; // This is used in the finally() function @@ -134,17 +136,16 @@ void EvolveMomentum::finally(const Options& state) { // Include a correction term for electromagnetic simulations const Field3D Apar = get(state["fields"]["Apar"]); - const Field3D density_source = species.isSet("density_source") ? - get(species["density_source"]) - : zeroFrom(Apar); + const Field3D density_source = species.isSet("density_source") + ? get(species["density_source"]) + : zeroFrom(Apar); Field3D dummy; - + // This is Z * Apar * dn/dt, keeping just leading order terms Field3D dndt = density_source - - FV::Div_par_mod(N, V, fastest_wave, dummy) - - Div_n_bxGrad_f_B_XPPM(N, phi, bndry_flux, poloidal_flows, true) - ; + - FV::Div_par_mod(N, V, fastest_wave, dummy) + - Div_n_bxGrad_f_B_XPPM(N, phi, bndry_flux, poloidal_flows, true); if (low_n_diffuse_perp) { dndt += Div_Perp_Lap_FV_Index( density_floor / softFloor(N, 1e-3 * density_floor), N); @@ -170,8 +171,10 @@ void EvolveMomentum::finally(const Options& state) { // - Density floor should be consistent with calculation of V // otherwise energy conservation is affected // - using the same operator as in density and pressure equations doesn't work - ddt(NV) -= AA * FV::Div_par_fvv(Nlim, V, fastest_wave, fix_momentum_boundary_flux); - + ddt(NV) -= AA + * FV::Div_par_fvv(Nlim, V, fastest_wave, + fix_momentum_boundary_flux); + // Parallel pressure gradient if (species.isSet("pressure")) { Field3D P = get(species["pressure"]); @@ -192,7 +195,8 @@ void EvolveMomentum::finally(const Options& state) { if (species.isSet("low_n_coeff")) { // Low density parallel diffusion Field3D low_n_coeff = get(species["low_n_coeff"]); - ddt(NV) += FV::Div_par_K_Grad_par(low_n_coeff * V, N) + FV::Div_par_K_Grad_par(low_n_coeff * Nlim, V); + ddt(NV) += FV::Div_par_K_Grad_par(low_n_coeff * V, N) + + FV::Div_par_K_Grad_par(low_n_coeff * Nlim, V); } if (low_n_diffuse_perp) { @@ -298,24 +302,26 @@ void EvolveMomentum::outputVars(Options& state) { auto rho_s0 = get(state["rho_s0"]); if (flow_xlow.isAllocated()) { - set_with_attrs(state[fmt::format("mf{}_tot_xlow", name)], flow_xlow, - {{"time_dimension", "t"}, - {"units", "N"}, - {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"standard_name", "momentum flow"}, - {"long_name", name + " momentum flow in X. Note: May be incomplete."}, - {"species", name}, - {"source", "evolve_momentum"}}); + set_with_attrs( + state[fmt::format("mf{}_tot_xlow", name)], flow_xlow, + {{"time_dimension", "t"}, + {"units", "N"}, + {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum flow"}, + {"long_name", name + " momentum flow in X. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_momentum"}}); } if (flow_ylow.isAllocated()) { - set_with_attrs(state[fmt::format("mf{}_tot_ylow", name)], flow_ylow, - {{"time_dimension", "t"}, - {"units", "N"}, - {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"standard_name", "momentum flow"}, - {"long_name", name + " momentum flow in Y. Note: May be incomplete."}, - {"species", name}, - {"source", "evolve_momentum"}}); + set_with_attrs( + state[fmt::format("mf{}_tot_ylow", name)], flow_ylow, + {{"time_dimension", "t"}, + {"units", "N"}, + {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum flow"}, + {"long_name", name + " momentum flow in Y. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_momentum"}}); } } } diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index ebfb606f8..90666cfe6 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -7,6 +7,7 @@ #include #include #include +#include #include "../include/div_ops.hxx" #include "../include/evolve_pressure.hxx" diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 6f632c8f9..19b74ef7c 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -4,6 +4,7 @@ #include #include #include +#include #include "../include/div_ops.hxx" #include "../include/hermes_build_config.hxx" diff --git a/src/recycling.cxx b/src/recycling.cxx index 8cf6cf5f6..603c75f4c 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -8,6 +8,7 @@ #include // for trim, strsplit #include +#include using bout::globals::mesh; @@ -54,8 +55,9 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) for (const auto& species : species_list) { std::string from = trim(species, " \t\r()"); // The species name in the list - if (from.empty()) + if (from.empty()) { continue; // Missing + } // Get the options for this species Options& from_options = alloptions[from]; @@ -272,7 +274,7 @@ void Recycling::transform_impl(GuardedOptions& state) { BoutReal tisheath = (T[i] + T[ig]) * 0.5; BoutReal visheath = (V[i] + V[ig]) * 0.5; BoutReal sheath_ion_heat_flow = - abs(gamma_i * nisheath * tisheath * visheath * daparsheath / volume); + std::abs(gamma_i * nisheath * tisheath * visheath * daparsheath / volume); // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] @@ -366,7 +368,7 @@ void Recycling::transform_impl(GuardedOptions& state) { BoutReal tisheath = (T[i] + T[ig]) * 0.5; BoutReal visheath = (V[i] + V[ig]) * 0.5; BoutReal sheath_ion_heat_flow = - abs(gamma_i * nisheath * tisheath * visheath * daparsheath / volume); + std::abs(gamma_i * nisheath * tisheath * visheath * daparsheath / volume); // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index 425c1f6f6..cde39db88 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -7,6 +7,7 @@ #include #include #include +#include using bout::globals::mesh; diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 116fb4c22..2a04cde9d 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -4,11 +4,13 @@ #include "../include/hermes_utils.hxx" #include -#include -#include #include #include +#include +#include #include +#include +#include using bout::globals::mesh; @@ -55,16 +57,16 @@ Vector3D Grad_perp_XZ(const Field3D& f) { auto xp = i.xp(); auto xm = i.xm(); // DDX(f) - result.x[i] = ((4. * f[xp] + f[xp.zp()] + f[xp.zm()]) - - (4. * f[xm] + f[xm.zp()] + f[xm.zm()])) - / (12. * dx[i]); + result.x[i] = + ((4. * f[xp] + f[xp.zp()] + f[xp.zm()]) - (4. * f[xm] + f[xm.zp()] + f[xm.zm()])) + / (12. * dx[i]); // DDZ(f) auto zp = i.zp(); auto zm = i.zm(); - result.z[i] = ((4. * f[zp] + f[zp.xp()] + f[zp.xm()]) - - (4. * f[zm] + f[zm.xp()] + f[zm.xm()])) - / (12. * dz[i]); + result.z[i] = + ((4. * f[zp] + f[zp.xp()] + f[zp.xm()]) - (4. * f[zm] + f[zm.xp()] + f[zm.xm()])) + / (12. * dz[i]); } result.setLocation(f.getLocation()); @@ -91,15 +93,16 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) .withDefault(true); exb_advection_simplified = options["exb_advection_simplified"] - .doc("Simplify nonlinear ExB advection form?") - .withDefault(true); + .doc("Simplify nonlinear ExB advection form?") + .withDefault(true); diamagnetic = options["diamagnetic"].doc("Include diamagnetic current?").withDefault(true); - sheath_boundary = options["sheath_boundary"] - .doc("Set potential to j=0 sheath at radial boundaries? (default = 0)") - .withDefault(false); + sheath_boundary = + options["sheath_boundary"] + .doc("Set potential to j=0 sheath at radial boundaries? (default = 0)") + .withDefault(false); diamagnetic_polarisation = options["diamagnetic_polarisation"] @@ -128,10 +131,9 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) .withDefault(false); include_viscosity = options.isSet("viscosity"); // Only include if set - viscosity = options["viscosity"] - .doc("Kinematic viscosity [m^2/s]") - .withDefault(0.0) - / (Lnorm * Lnorm * Omega_ci); + viscosity = + options["viscosity"].doc("Kinematic viscosity [m^2/s]").withDefault(0.0) + / (Lnorm * Lnorm * Omega_ci); viscosity.applyBoundary("dirichlet"); hyper_z = options["hyper_z"].doc("Hyper-viscosity in Z. < 0 -> off").withDefault(-1.0); @@ -154,12 +156,12 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) .withDefault(false); phi_sheath_dissipation = options["phi_sheath_dissipation"] - .doc("Add dissipation when phi < 0.0 at the sheath") - .withDefault(false); + .doc("Add dissipation when phi < 0.0 at the sheath") + .withDefault(false); damp_core_vorticity = options["damp_core_vorticity"] - .doc("Damp vorticity at the core boundary?") - .withDefault(false); + .doc("Damp vorticity at the core boundary?") + .withDefault(false); // Add phi to restart files so that the value in the boundaries // is restored on restart. This is done even when phi is not evolving, @@ -186,8 +188,9 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) phi_boundary_last_update = -1.; phi_core_averagey = options["phi_core_averagey"] - .doc("Average phi core boundary in Y?") - .withDefault(false) and mesh->periodicY(mesh->xstart); + .doc("Average phi core boundary in Y?") + .withDefault(false) + and mesh->periodicY(mesh->xstart); phi_boundary_timescale = options["phi_boundary_timescale"] .doc("Timescale for phi boundary relaxation [seconds]") @@ -237,9 +240,8 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) Bsq = SQ(coord->Bxy); - diagnose = options["diagnose"] - .doc("Output additional diagnostics?") - .withDefault(false); + diagnose = + options["diagnose"].doc("Output additional diagnostics?").withDefault(false); if (diamagnetic or diamagnetic_polarisation) { // FIXME: These will only be read if BOTH charge and pressure (and possibly AA) are @@ -329,16 +331,18 @@ Field3D Vorticity::calculateDivJdia(Field3D phi, GuardedOptions allspecies) { // Note: We need boundary conditions on P, so apply the same // free boundary condition as sheath_boundary. if (P.hasParallelSlices()) { - Field3D &P_ydown = P.ydown(); - Field3D &P_yup = P.yup(); + Field3D& P_ydown = P.ydown(); + Field3D& P_yup = P.yup(); for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { - P_ydown(r.ind, mesh->ystart - 1, jz) = 2 * P(r.ind, mesh->ystart, jz) - P_yup(r.ind, mesh->ystart + 1, jz); + P_ydown(r.ind, mesh->ystart - 1, jz) = + 2 * P(r.ind, mesh->ystart, jz) - P_yup(r.ind, mesh->ystart + 1, jz); } } for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { - P_yup(r.ind, mesh->yend + 1, jz) = 2 * P(r.ind, mesh->yend, jz) - P_ydown(r.ind, mesh->yend - 1, jz); + P_yup(r.ind, mesh->yend + 1, jz) = + 2 * P(r.ind, mesh->yend, jz) - P_ydown(r.ind, mesh->yend - 1, jz); } } } else { @@ -415,10 +419,7 @@ void Vorticity::transform_impl(GuardedOptions& state) { MPI_Comm comm_inner = mesh->getYcomm(0); int np; MPI_Comm_size(comm_inner, &np); - MPI_Allreduce(&philocal, - &phivalue, - 1, MPI_DOUBLE, - MPI_SUM, comm_inner); + MPI_Allreduce(&philocal, &phivalue, 1, MPI_DOUBLE, MPI_SUM, comm_inner); phivalue /= (np * mesh->LocalNz * mesh->LocalNy); } @@ -534,7 +535,7 @@ void Vorticity::transform_impl(GuardedOptions& state) { } if (mesh->xend < mesh->LocalNx - 2) { for (int k = 0; k < mesh->LocalNz; k++) { - + // Note: This seems to make a difference, but don't know why. // Without this, get convergence failures with no apparent instability // (all fields apparently smooth, well behaved) @@ -620,28 +621,32 @@ void Vorticity::transform_impl(GuardedOptions& state) { // Note: The below calculation requires phi derivatives at the Y boundaries // Setting to free boundaries if (phi.hasParallelSlices()) { - Field3D &phi_ydown = phi.ydown(); - Field3D &phi_yup = phi.yup(); + Field3D& phi_ydown = phi.ydown(); + Field3D& phi_yup = phi.yup(); for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_ydown(r.ind, mesh->ystart - 1, jz) = 2 * phi(r.ind, mesh->ystart, jz) - phi_yup(r.ind, mesh->ystart + 1, jz); + phi_ydown(r.ind, mesh->ystart - 1, jz) = + 2 * phi(r.ind, mesh->ystart, jz) - phi_yup(r.ind, mesh->ystart + 1, jz); } } for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_yup(r.ind, mesh->yend + 1, jz) = 2 * phi(r.ind, mesh->yend, jz) - phi_ydown(r.ind, mesh->yend - 1, jz); + phi_yup(r.ind, mesh->yend + 1, jz) = + 2 * phi(r.ind, mesh->yend, jz) - phi_ydown(r.ind, mesh->yend - 1, jz); } } } else { Field3D phi_fa = toFieldAligned(phi); for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_fa(r.ind, mesh->ystart - 1, jz) = 2 * phi_fa(r.ind, mesh->ystart, jz) - phi_fa(r.ind, mesh->ystart + 1, jz); + phi_fa(r.ind, mesh->ystart - 1, jz) = + 2 * phi_fa(r.ind, mesh->ystart, jz) - phi_fa(r.ind, mesh->ystart + 1, jz); } } for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_fa(r.ind, mesh->yend + 1, jz) = 2 * phi_fa(r.ind, mesh->yend, jz) - phi_fa(r.ind, mesh->yend - 1, jz); + phi_fa(r.ind, mesh->yend + 1, jz) = + 2 * phi_fa(r.ind, mesh->yend, jz) - phi_fa(r.ind, mesh->yend - 1, jz); } } phi = fromFieldAligned(phi_fa); @@ -723,7 +728,7 @@ void Vorticity::transform_impl(GuardedOptions& state) { //phi_hat.applyBoundary("neumann"); // Note: Grad_perp_XZ ignores the DDY terms - viscous_heating = - viscosity * (Grad_perp_XZ(Vort) * Grad_perp_XZ(phi_hat)); + viscous_heating = -viscosity * (Grad_perp_XZ(Vort) * Grad_perp_XZ(phi_hat)); // Weight the heating by mass density of charged species // This is an approximation motivated because the perpendicular @@ -734,7 +739,8 @@ void Vorticity::transform_impl(GuardedOptions& state) { for (const auto& kv : allspecies.getChildren()) { const GuardedOptions species = kv.second; - if (!(species.isSet("charge") and species.isSet("AA") and species.isSet("density"))) { + if (!(species.isSet("charge") and species.isSet("AA") + and species.isSet("density"))) { continue; // No charge or mass -> no current } if (fabs(get(species["charge"])) < 1e-5) { @@ -749,7 +755,8 @@ void Vorticity::transform_impl(GuardedOptions& state) { for (const auto& kv : allspecies.getChildren()) { GuardedOptions species = allspecies[kv.first]; // Note: need non-const - if (!(species.isSet("charge") and species.isSet("AA") and species.isSet("density"))) { + if (!(species.isSet("charge") and species.isSet("AA") + and species.isSet("density"))) { continue; // No charge or mass -> no current } if (fabs(get(species["charge"])) < 1e-5) { @@ -786,8 +793,7 @@ void Vorticity::finally(const Options& state) { // Because this is implemented in terms of an operation on the result // of an operation, we need to communicate and the resulting stencil is // wider than the simple form. - ddt(Vort) -= - Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, bndry_flux, poloidal_flows); + ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, bndry_flux, poloidal_flows); // V_ExB dot Grad(Pi) Field3D vEdotGradPi = bracket(phi, Pi_hat, BRACKET_ARAKAWA); @@ -800,8 +806,8 @@ void Vorticity::finally(const Options& state) { mesh->communicate(vEdotGradPi, DelpPhi_2B2); ddt(Vort) -= FV::Div_a_Grad_perp(0.5 * average_atomic_mass / Bsq, vEdotGradPi); - ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi_hat, bndry_flux, - poloidal_flows); + ddt(Vort) -= + Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi_hat, bndry_flux, poloidal_flows); } } @@ -871,7 +877,7 @@ void Vorticity::finally(const Options& state) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { auto i = indexAt(phi_fa, r.ind, mesh->ystart, jz); - BoutReal phisheath = 0.5*(phi_fa[i] + phi_fa[i.ym()]); + BoutReal phisheath = 0.5 * (phi_fa[i] + phi_fa[i.ym()]); dissipation[i] = -floor(-phisheath, 0.0); } } @@ -879,7 +885,7 @@ void Vorticity::finally(const Options& state) { for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { auto i = indexAt(phi_fa, r.ind, mesh->yend, jz); - BoutReal phisheath = 0.5*(phi_fa[i] + phi_fa[i.yp()]); + BoutReal phisheath = 0.5 * (phi_fa[i] + phi_fa[i.yp()]); dissipation[i] = -floor(-phisheath, 0.0); } } diff --git a/tests/mms_operator/main.cxx b/tests/mms_operator/main.cxx index eab40edba..480cf1ef1 100644 --- a/tests/mms_operator/main.cxx +++ b/tests/mms_operator/main.cxx @@ -1,3 +1,4 @@ +#include "bout/bout.hxx" #include "bout/field_factory.hxx" // #include "hermes-2.hxx" diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 5d0d12de0..6dbbd6200 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -43,6 +44,9 @@ public: GlobalNx = nx; GlobalNy = ny; GlobalNz = nz; + GlobalNxNoBoundaries = nx - 2; + GlobalNyNoBoundaries = ny - 2; + GlobalNzNoBoundaries = nz; LocalNx = nx; LocalNy = ny; LocalNz = nz; @@ -66,15 +70,13 @@ public: xend = nx - 2; ystart = 1; yend = ny - 2; - zstart = 0; + zstart = 0; // no guards zend = nz - 1; StaggerGrids = false; // Unused variables periodicX = false; - NXPE = 1; - PE_XIND = 0; IncIntShear = false; maxregionblocksize = MAXREGIONBLOCKSIZE; @@ -109,10 +111,16 @@ public: return nullptr; } int wait(comm_handle UNUSED(handle)) override { return 0; } - int getNXPE() override { return 1; } - int getNYPE() override { return 1; } - int getXProcIndex() override { return 1; } - int getYProcIndex() override { return 1; } + int getNXPE() const override { return 1; } + int getNYPE() const override { return 1; } + int getNZPE() const override { return 1; } + int getXProcIndex() const override { return 0; } + int getYProcIndex() const override { return 0; } + int getZProcIndex() const override { return 0; } + int getProcIndex([[maybe_unused]] int X, [[maybe_unused]] int Y, + [[maybe_unused]] int Z) const override { + return 0; + } bool firstX() const override { return true; } bool lastX() const override { return true; } int sendXOut(BoutReal* UNUSED(buffer), int UNUSED(size), int UNUSED(tag)) override { @@ -131,8 +139,15 @@ public: } MPI_Comm getXcomm(int UNUSED(jy)) const override { return BoutComm::get(); } MPI_Comm getYcomm(int UNUSED(jx)) const override { return BoutComm::get(); } - bool periodicY(int UNUSED(jx)) const override { return true; } - bool periodicY(int UNUSED(jx), BoutReal& UNUSED(ts)) const override { return true; } + MPI_Comm getXZcomm() const override { return BoutComm::get(); } + + // Periodic Y + int ix_separatrix{1000000}; // separatrix index + + bool periodicY(int jx) const override { return jx < ix_separatrix; } + bool periodicY(int jx, BoutReal& UNUSED(ts)) const override { + return jx < ix_separatrix; + } int numberOfYBoundaries() const override { return 1; } std::pair hasBranchCutLower(int UNUSED(jx)) const override { return std::make_pair(false, 0.); @@ -164,20 +179,22 @@ public: } BoutReal GlobalX(int jx) const override { return jx; } BoutReal GlobalY(int jy) const override { return jy; } + BoutReal GlobalZ(int jz) const override { return jz; } BoutReal GlobalX(BoutReal jx) const override { return jx; } BoutReal GlobalY(BoutReal jy) const override { return jy; } + BoutReal GlobalZ(BoutReal jz) const override { return jz; } int getGlobalXIndex(int) const override { return 0; } int getGlobalXIndexNoBoundaries(int) const override { return 0; } int getGlobalYIndex(int y) const override { return y; } int getGlobalYIndexNoBoundaries(int y) const override { return y; } - int getGlobalZIndex(int) const override { return 0; } - int getGlobalZIndexNoBoundaries(int) const override { return 0; } + int getGlobalZIndex(int z) const override { return z; } + int getGlobalZIndexNoBoundaries(int z) const override { return z; } int getLocalXIndex(int) const override { return 0; } int getLocalXIndexNoBoundaries(int) const override { return 0; } int getLocalYIndex(int y) const override { return y; } int getLocalYIndexNoBoundaries(int y) const override { return y; } - int getLocalZIndex(int) const override { return 0; } - int getLocalZIndexNoBoundaries(int) const override { return 0; } + int getLocalZIndex(int z) const override { return z; } + int getLocalZIndexNoBoundaries(int z) const override { return z; } void initDerivs(Options* opt) { StaggerGrids = true; @@ -225,19 +242,18 @@ public: "RGN_OUTER_X"}; // Sum up and get unique points in the boundaries defined above - addRegion2D("RGN_BNDRY", - std::accumulate(begin(boundary_names), end(boundary_names), - Region{}, - [this](Region& a, const std::string& b) { - return a + getRegion2D(b); - }) - .unique()); + addRegion2D("RGN_BNDRY", std::accumulate(begin(boundary_names), end(boundary_names), + Region{}, + [&](Region a, const std::string& b) { + return std::move(a) + getRegion2D(b); + }) + .unique()); addRegion3D("RGN_BNDRY", std::accumulate(begin(boundary_names), end(boundary_names), Region{}, - [this](Region& a, const std::string& b) { - return a + getRegion3D(b); + [this](Region a, const std::string& b) { + return std::move(a) + getRegion3D(b); }) .unique()); addRegionPerp("RGN_BNDRY", diff --git a/tests/unit/reactions/C5pIzn.nc b/tests/unit/reactions/C5pIzn.nc index 7126e46c8..e74048dd0 100644 Binary files a/tests/unit/reactions/C5pIzn.nc and b/tests/unit/reactions/C5pIzn.nc differ diff --git a/tests/unit/reactions/C5pRec.nc b/tests/unit/reactions/C5pRec.nc index d07e707b7..906ab74b6 100644 Binary files a/tests/unit/reactions/C5pRec.nc and b/tests/unit/reactions/C5pRec.nc differ diff --git a/tests/unit/reactions/CIzn.nc b/tests/unit/reactions/CIzn.nc index ed539112f..67439c591 100644 Binary files a/tests/unit/reactions/CIzn.nc and b/tests/unit/reactions/CIzn.nc differ diff --git a/tests/unit/reactions/CRec.nc b/tests/unit/reactions/CRec.nc index ffa5a6543..db8b9f439 100644 Binary files a/tests/unit/reactions/CRec.nc and b/tests/unit/reactions/CRec.nc differ diff --git a/tests/unit/reactions/Cp3DCX.nc b/tests/unit/reactions/Cp3DCX.nc index dbcecb42e..d0856f488 100644 Binary files a/tests/unit/reactions/Cp3DCX.nc and b/tests/unit/reactions/Cp3DCX.nc differ diff --git a/tests/unit/reactions/Cp5TCX.nc b/tests/unit/reactions/Cp5TCX.nc index 2a1174635..6cfbb102e 100644 Binary files a/tests/unit/reactions/Cp5TCX.nc and b/tests/unit/reactions/Cp5TCX.nc differ diff --git a/tests/unit/reactions/CpHCXCX.nc b/tests/unit/reactions/CpHCXCX.nc index f9ffbd4ed..ffc1df6aa 100644 Binary files a/tests/unit/reactions/CpHCXCX.nc and b/tests/unit/reactions/CpHCXCX.nc differ diff --git a/tests/unit/reactions/CpIzn.nc b/tests/unit/reactions/CpIzn.nc index 3fa0cb037..5bd974f43 100644 Binary files a/tests/unit/reactions/CpIzn.nc and b/tests/unit/reactions/CpIzn.nc differ diff --git a/tests/unit/reactions/CpRec.nc b/tests/unit/reactions/CpRec.nc index e27d3565f..967ed72e0 100644 Binary files a/tests/unit/reactions/CpRec.nc and b/tests/unit/reactions/CpRec.nc differ diff --git a/tests/unit/reactions/DDpCX.nc b/tests/unit/reactions/DDpCX.nc index 55d17b292..dc2449842 100644 Binary files a/tests/unit/reactions/DDpCX.nc and b/tests/unit/reactions/DDpCX.nc differ diff --git a/tests/unit/reactions/DIzn.nc b/tests/unit/reactions/DIzn.nc index 7a8aa02de..e95545d21 100644 Binary files a/tests/unit/reactions/DIzn.nc and b/tests/unit/reactions/DIzn.nc differ diff --git a/tests/unit/reactions/DRec.nc b/tests/unit/reactions/DRec.nc index b7b7eaf55..0c01b4fb6 100644 Binary files a/tests/unit/reactions/DRec.nc and b/tests/unit/reactions/DRec.nc differ diff --git a/tests/unit/reactions/DTpCX.nc b/tests/unit/reactions/DTpCX.nc index d93e8472d..23acf30cb 100644 Binary files a/tests/unit/reactions/DTpCX.nc and b/tests/unit/reactions/DTpCX.nc differ diff --git a/tests/unit/reactions/HDpCX.nc b/tests/unit/reactions/HDpCX.nc index 749b73535..f8e299b3a 100644 Binary files a/tests/unit/reactions/HDpCX.nc and b/tests/unit/reactions/HDpCX.nc differ diff --git a/tests/unit/reactions/HHpCX.nc b/tests/unit/reactions/HHpCX.nc index 862ef5c3b..2eac9143e 100644 Binary files a/tests/unit/reactions/HHpCX.nc and b/tests/unit/reactions/HHpCX.nc differ diff --git a/tests/unit/reactions/HHpCX_noNeutralMomGain.nc b/tests/unit/reactions/HHpCX_noNeutralMomGain.nc index 8eed684db..063f4414e 100644 Binary files a/tests/unit/reactions/HHpCX_noNeutralMomGain.nc and b/tests/unit/reactions/HHpCX_noNeutralMomGain.nc differ diff --git a/tests/unit/reactions/HIzn.nc b/tests/unit/reactions/HIzn.nc index f121701f9..6f901eb46 100644 Binary files a/tests/unit/reactions/HIzn.nc and b/tests/unit/reactions/HIzn.nc differ diff --git a/tests/unit/reactions/HRec.nc b/tests/unit/reactions/HRec.nc index 80d0d19b9..0e3e997d4 100644 Binary files a/tests/unit/reactions/HRec.nc and b/tests/unit/reactions/HRec.nc differ diff --git a/tests/unit/reactions/HeIzn01.nc b/tests/unit/reactions/HeIzn01.nc index ed511559c..e287e4bdb 100644 Binary files a/tests/unit/reactions/HeIzn01.nc and b/tests/unit/reactions/HeIzn01.nc differ diff --git a/tests/unit/reactions/HeRec10.nc b/tests/unit/reactions/HeRec10.nc index e91d2f92b..8cfd7f0e3 100644 Binary files a/tests/unit/reactions/HeRec10.nc and b/tests/unit/reactions/HeRec10.nc differ diff --git a/tests/unit/reactions/Li2pIzn.nc b/tests/unit/reactions/Li2pIzn.nc index 14c3a472f..cf74578cb 100644 Binary files a/tests/unit/reactions/Li2pIzn.nc and b/tests/unit/reactions/Li2pIzn.nc differ diff --git a/tests/unit/reactions/Li2pRec.nc b/tests/unit/reactions/Li2pRec.nc index 4ee7bc4d3..4978d2edb 100644 Binary files a/tests/unit/reactions/Li2pRec.nc and b/tests/unit/reactions/Li2pRec.nc differ diff --git a/tests/unit/reactions/LiIzn.nc b/tests/unit/reactions/LiIzn.nc index 9876d7511..8c64453a8 100644 Binary files a/tests/unit/reactions/LiIzn.nc and b/tests/unit/reactions/LiIzn.nc differ diff --git a/tests/unit/reactions/LiRec.nc b/tests/unit/reactions/LiRec.nc index e7cc7ddad..5e590f92c 100644 Binary files a/tests/unit/reactions/LiRec.nc and b/tests/unit/reactions/LiRec.nc differ diff --git a/tests/unit/reactions/Lip2DCX.nc b/tests/unit/reactions/Lip2DCX.nc index b886daaff..054abdf13 100644 Binary files a/tests/unit/reactions/Lip2DCX.nc and b/tests/unit/reactions/Lip2DCX.nc differ diff --git a/tests/unit/reactions/Lip3TCX.nc b/tests/unit/reactions/Lip3TCX.nc index 27811144b..322ab55eb 100644 Binary files a/tests/unit/reactions/Lip3TCX.nc and b/tests/unit/reactions/Lip3TCX.nc differ diff --git a/tests/unit/reactions/LipHCX.nc b/tests/unit/reactions/LipHCX.nc index 2221e82d4..7f4e95e29 100644 Binary files a/tests/unit/reactions/LipHCX.nc and b/tests/unit/reactions/LipHCX.nc differ diff --git a/tests/unit/reactions/LipIzn.nc b/tests/unit/reactions/LipIzn.nc index 0e4f5f4ed..8fe873cf6 100644 Binary files a/tests/unit/reactions/LipIzn.nc and b/tests/unit/reactions/LipIzn.nc differ diff --git a/tests/unit/reactions/LipRec.nc b/tests/unit/reactions/LipRec.nc index ad9d4ac37..6ef030b98 100644 Binary files a/tests/unit/reactions/LipRec.nc and b/tests/unit/reactions/LipRec.nc differ diff --git a/tests/unit/reactions/Ne9pIzn.nc b/tests/unit/reactions/Ne9pIzn.nc index 9d39c5e27..69ca7d261 100644 Binary files a/tests/unit/reactions/Ne9pIzn.nc and b/tests/unit/reactions/Ne9pIzn.nc differ diff --git a/tests/unit/reactions/Ne9pRec.nc b/tests/unit/reactions/Ne9pRec.nc index be62b8c5f..e0d5bd936 100644 Binary files a/tests/unit/reactions/Ne9pRec.nc and b/tests/unit/reactions/Ne9pRec.nc differ diff --git a/tests/unit/reactions/NeIzn.nc b/tests/unit/reactions/NeIzn.nc index 24e1eac1e..848a466e3 100644 Binary files a/tests/unit/reactions/NeIzn.nc and b/tests/unit/reactions/NeIzn.nc differ diff --git a/tests/unit/reactions/NeRec.nc b/tests/unit/reactions/NeRec.nc index 5d72d2be3..30e30c455 100644 Binary files a/tests/unit/reactions/NeRec.nc and b/tests/unit/reactions/NeRec.nc differ diff --git a/tests/unit/reactions/Nep5DCX.nc b/tests/unit/reactions/Nep5DCX.nc index dfa480a00..8974797f2 100644 Binary files a/tests/unit/reactions/Nep5DCX.nc and b/tests/unit/reactions/Nep5DCX.nc differ diff --git a/tests/unit/reactions/Nep9TCX.nc b/tests/unit/reactions/Nep9TCX.nc index c47ea7e00..9d0a194b0 100644 Binary files a/tests/unit/reactions/Nep9TCX.nc and b/tests/unit/reactions/Nep9TCX.nc differ diff --git a/tests/unit/reactions/NepHCX.nc b/tests/unit/reactions/NepHCX.nc index f4238141e..6b8ed67f4 100644 Binary files a/tests/unit/reactions/NepHCX.nc and b/tests/unit/reactions/NepHCX.nc differ diff --git a/tests/unit/reactions/NepIzn.nc b/tests/unit/reactions/NepIzn.nc index cd55e89f5..1836bcbe6 100644 Binary files a/tests/unit/reactions/NepIzn.nc and b/tests/unit/reactions/NepIzn.nc differ diff --git a/tests/unit/reactions/NepRec.nc b/tests/unit/reactions/NepRec.nc index 75b8e616f..16b2b096c 100644 Binary files a/tests/unit/reactions/NepRec.nc and b/tests/unit/reactions/NepRec.nc differ diff --git a/tests/unit/reactions/THpCX.nc b/tests/unit/reactions/THpCX.nc index 2aff8ab4c..1dafa7a82 100644 Binary files a/tests/unit/reactions/THpCX.nc and b/tests/unit/reactions/THpCX.nc differ diff --git a/tests/unit/reactions/TIzn.nc b/tests/unit/reactions/TIzn.nc index 3aefd0519..f80b64f86 100644 Binary files a/tests/unit/reactions/TIzn.nc and b/tests/unit/reactions/TIzn.nc differ diff --git a/tests/unit/reactions/TRec.nc b/tests/unit/reactions/TRec.nc index 3b1cf56a7..64a668722 100644 Binary files a/tests/unit/reactions/TRec.nc and b/tests/unit/reactions/TRec.nc differ diff --git a/tests/unit/reactions/TTpCX.nc b/tests/unit/reactions/TTpCX.nc index 566926f82..559c6d7e5 100644 Binary files a/tests/unit/reactions/TTpCX.nc and b/tests/unit/reactions/TTpCX.nc differ diff --git a/tests/unit/test_reactions.hxx b/tests/unit/test_reactions.hxx index 3fb65361a..59af61b59 100644 --- a/tests/unit/test_reactions.hxx +++ b/tests/unit/test_reactions.hxx @@ -47,7 +47,7 @@ class ReactionTest : public FakeMeshFixture_tmpl<8, 8, 8> { protected: ReactionTest(std::string lbl, std::string reaction_str) - : lbl(lbl), parser(reaction_str){}; + : lbl(lbl), parser(reaction_str) {}; std::string lbl; ReactionParser parser; @@ -91,8 +91,8 @@ protected: break; case linfunc_axis::z: axis_str = "z"; - axis_min = TWOPI * (mesh->zstart) / static_cast(mesh->LocalNz); - axis_max = TWOPI * (mesh->zend) / static_cast(mesh->LocalNz); + axis_min = TWOPI * (mesh->zstart); + axis_max = TWOPI * (mesh->zend); break; default: axis_min = axis_max = 0; @@ -168,6 +168,9 @@ protected: void sources_regression_test(bool compare_all_values = true, const int ignore_last_n_sigfigs = 6) { + // Uncomment below to update the reference data. + //generate_data(); + // Read reference state Options ref_state = bout::OptionsIO::create(ref_data_path())->read();