diff --git a/CMakeLists.txt b/CMakeLists.txt index 907f3fef5..ee7a1ddb8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -69,6 +69,7 @@ set(HERMES_SOURCES src/neutral_full_velocity.cxx src/electromagnetic.cxx src/electron_force_balance.cxx + src/evolve_anisotropic_pressure.cxx src/evolve_density.cxx src/evolve_energy.cxx src/evolve_pressure.cxx @@ -128,6 +129,7 @@ set(HERMES_SOURCES include/div_ops.hxx include/electromagnetic.hxx include/electron_force_balance.hxx + include/evolve_anisotropic_pressure.hxx include/evolve_density.hxx include/evolve_energy.hxx include/evolve_momentum.hxx diff --git a/include/evolve_anisotropic_pressure.hxx b/include/evolve_anisotropic_pressure.hxx new file mode 100644 index 000000000..e4ba3f14c --- /dev/null +++ b/include/evolve_anisotropic_pressure.hxx @@ -0,0 +1,109 @@ +#pragma once +#ifndef EVOLVE_ANISOTROPIC_PRESSURE_H +#define EVOLVE_ANISOTROPIC_PRESSURE_H + +#include "../include/hermes_utils.hxx" +#include "component.hxx" +#include + +/// Evolves species anisotropic pressure (P_parallel & P_perp) in time. +/// This formulation evolves +/// E Energy, combining thermal and parallel kinetic +/// PA Pressure anisotropy, PA = P_par - P_perp +/// +/// # Mesh inputs +/// +/// P_src A source of pressure, in Pascals per second +/// This can be over-ridden by the `source` option setting. +/// +struct EvolveAnisotropicPressure : public Component { + /// + /// # Inputs + /// + /// - + /// - bndry_flux Allow flows through radial boundaries? Default is true + /// - density_floor Minimum density floor. Default 1e-5 normalised units. + /// - diagnose Output additional diagnostic fields? + /// - poloidal_flows Include poloidal ExB flows? Default is true + /// - precon Enable preconditioner? Note: solver may not use it even if + /// enabled. + /// - p_div_v Use p * Div(v) form? Default is v * Grad(p) form + /// - thermal_conduction Include parallel heat conduction? Default is true + /// + /// - P e.g. "Pe", "Pd+" + /// - source Source of pressure [Pa / s]. + /// NOTE: This overrides mesh input P_src + /// - source_only_in_core Zero the source outside the closed field-line + /// region? + /// - neumann_boundary_average_z Apply Neumann boundaries with Z average? + /// + EvolveAnisotropicPressure(std::string name, Options& options, Solver* solver); + + /// + /// Optional inputs + /// + /// - species + /// - + /// - velocity. Must have sound_speed or temperature + /// - energy_source + /// - collision_rate (needed if thermal_conduction on) + /// - fields + /// - phi Electrostatic potential -> ExB drift + /// + void finally(const Options& state) override; + + void outputVars(Options& state) override; + +private: + std::string name; ///< Short name of the species e.g. h+ + + BoutReal adiabatic_index; ///< Ratio of specific heats, γ = Cp / Cv + BoutReal Cv; ///< /// Heat capacity at constant volume (3/2 for ideal monatomic gas) + + Field3D E; ///< Energy (normalised): P + 1/2 m n v^2 + Field3D PA; ///< Pressure anisotropy: P_par - P_perp + + Field3D P; ///< Pressure (normalised): P = (2 P_perp + P_par) / 3 + Field3D P_perp, P_par; ///< Perpendicular and parallel pressure (normalised) + + Field3D T, N; ///< Temperature, density + + bool bndry_flux; + bool neumann_boundary_average_z; ///< Apply neumann boundary with Z average? + bool poloidal_flows; + + BoutReal density_floor; ///< Minimum density for calculating T + BoutReal temperature_floor; ///< Low temperature scale for low_T_diffuse_perp + BoutReal pressure_floor; ///< When non-zero pressure is needed + + Field3D source, final_source; ///< External pressure source + Field3D Sp; ///< Total pressure source + FieldGeneratorPtr source_prefactor_function; + + bool diagnose; ///< Output additional diagnostics? + BoutReal source_normalisation; ///< Normalisation factor [Pa/s] + BoutReal time_normalisation; ///< Normalisation factor [s] + bool source_time_dependent; ///< Is the input source time dependent? + + /// Inputs + /// - species + /// - + /// - density + /// + /// Sets + /// - species + /// - + /// - pressure <- This is (2pressure_perp + pressure_par)/3 + /// - temperature Requires density + /// - pressure_par + /// - pressure_perp + /// + void transform_impl(GuardedOptions& state) override; +}; + +namespace { +RegisterComponent + registercomponentevolveanisotropicpressure("evolve_anisotropic_pressure"); +} + +#endif // EVOLVE_ANISOTROPIC_PRESSURE_H diff --git a/src/evolve_anisotropic_pressure.cxx b/src/evolve_anisotropic_pressure.cxx new file mode 100644 index 000000000..ba284d6c8 --- /dev/null +++ b/src/evolve_anisotropic_pressure.cxx @@ -0,0 +1,351 @@ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../include/div_ops.hxx" +#include "../include/evolve_anisotropic_pressure.hxx" +#include "../include/hermes_build_config.hxx" +#include "../include/hermes_utils.hxx" + +using bout::globals::mesh; + +EvolveAnisotropicPressure::EvolveAnisotropicPressure(std::string name, + Options& alloptions, Solver* solver) + : Component({ + readOnly("species:{name}:density", Regions::Interior), + readOnly("species:{name}:velocity", Regions::Interior), + readOnly("species:{name}:AA"), + readWrite("species:{name}:pressure"), + readWrite("species:{name}:pressure_perp"), + readWrite("species:{name}:pressure_par"), + readWrite("species:{name}:temperature"), + readWrite("species:{name}:temperature_perp"), + readWrite("species:{name}:temperature_par"), + }), + name(name) { + AUTO_TRACE(); + + auto& options = alloptions[name]; + + adiabatic_index = + options["adiabatic_index"] + .doc("Ratio of specific heats γ = Cp/Cv [5/3 for monatomic ideal gas]") + .withDefault(5. / 3); + Cv = 1. / (adiabatic_index - 1.); + + density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7); + + temperature_floor = options["temperature_floor"] + .doc("Low temperature scale for low_T_diffuse_perp") + .withDefault(0.1) + / get(alloptions["units"]["eV"]); + + pressure_floor = density_floor * temperature_floor; + + // Evolve the energy + solver->add(E, std::string("E") + name); + // Evolve the pressure anisotropy Ppar - Pperp + solver->add(PA, std::string("PA") + name); + + bndry_flux = options["bndry_flux"] + .doc("Allow flows through radial boundaries") + .withDefault(true); + + poloidal_flows = + options["poloidal_flows"].doc("Include poloidal ExB flow").withDefault(true); + + diagnose = options["diagnose"] + .doc("Save additional output diagnostics") + .withDefault(false); + + const Options& units = alloptions["units"]; + const BoutReal Nnorm = units["inv_meters_cubed"]; + const BoutReal Tnorm = units["eV"]; + const BoutReal Omega_ci = 1. / units["seconds"].as(); + + auto& p_options = alloptions[std::string("P") + name]; + source_normalisation = + SI::qe * Nnorm * Tnorm * Omega_ci; // [Pa/s] or [W/m^3] if converted to energy + time_normalisation = 1. / Omega_ci; // [s] + + // Try to read the pressure source from the mesh + // Units of Pascals per second + source = 0.0; + mesh->get(source, std::string("P") + name + "_src"); + // Allow the user to override the source + source = p_options["source"] + .doc(std::string("Source term in ddt(P") + name + + std::string("). Units [Pa/s], note P = 2/3 E")) + .withDefault(source) + / (source_normalisation); + + source_time_dependent = p_options["source_time_dependent"] + .doc("Use a time-dependent source?") + .withDefault(false); + + // If time dependent, parse the function with respect to time from the input file + if (source_time_dependent) { + auto str = p_options["source_prefactor"] + .doc("Time-dependent function of multiplier on ddt(P" + name + + std::string(") source.")) + .as(); + source_prefactor_function = FieldFactory::get()->parse(str, &p_options); + } + + if (p_options["source_only_in_core"] + .doc("Zero the source outside the closed field-line region?") + .withDefault(false)) { + for (int x = mesh->xstart; x <= mesh->xend; x++) { + if (!mesh->periodicY(x)) { + // Not periodic, so not in core + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { + source(x, y, z) = 0.0; + } + } + } + } + } + + neumann_boundary_average_z = p_options["neumann_boundary_average_z"] + .doc("Apply neumann boundary with Z average?") + .withDefault(false); +} + +void EvolveAnisotropicPressure::transform_impl(GuardedOptions& state) { + AUTO_TRACE(); + + mesh->communicate(E, PA); + + auto species = state["species"][name]; + N = getNoBoundary(species["density"]); + const Field3D V = getNoBoundary(species["velocity"]); + const BoutReal AA = get(species["AA"]); + + // Calculate pressure + // E = Cv * P + (1/2) m n v^2 + // P = (2 P_perp + P_par) / 3 + // PA = P_par - P_perp + // -> P = ( 2 (P_par - PA) + P_par ) / 3 = P_par - (2/3) PA + // P_par = P + (2/3) PA + // + // -> P = ( 2 P_perp + (PA + P_perp) ) / 3 = P_perp + (1/3) PA + // P_perp = P - (1/3) PA + // + // => Limit -(3/2) P < PA < 3 P + P.allocate(); + P_par.allocate(); + P_perp.allocate(); + BOUT_FOR(i, P.getRegion("RGN_ALL")) { + P[i] = (E[i] - 0.5 * AA * N[i] * SQ(V[i])) / Cv; + if (P[i] < 0.0) { + P[i] = 0.0; + } + + // Limited anisotropy, ensuring that P_perp, P_par >= 0 + BoutReal pa = BOUTMAX(-(3. / 2) * P[i], BOUTMIN(PA[i], 3. * P[i])); + + P_par[i] = P[i] + (2. / 3) * pa; + P_perp[i] = P[i] - (1. / 3) * pa; + } + P.applyBoundary("neumann"); + P_par.applyBoundary("neumann"); + P_perp.applyBoundary("neumann"); + + if (neumann_boundary_average_z) { + // Take Z (usually toroidal) average and apply as X (radial) boundary condition + + auto applyBoundaryAverageZ = [](Field3D& f) { + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal favg = 0.0; // Average P in Z + for (int k = 0; k < mesh->LocalNz; k++) { + favg += f(mesh->xstart, j, k); + } + favg /= mesh->LocalNz; + + // Apply boundary condition + for (int k = 0; k < mesh->LocalNz; k++) { + f(mesh->xstart - 1, j, k) = 2. * favg - f(mesh->xstart, j, k); + f(mesh->xstart - 2, j, k) = f(mesh->xstart - 1, j, k); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal favg = 0.0; // Average P in Z + for (int k = 0; k < mesh->LocalNz; k++) { + favg += f(mesh->xend, j, k); + } + favg /= mesh->LocalNz; + + for (int k = 0; k < mesh->LocalNz; k++) { + f(mesh->xend + 1, j, k) = 2. * favg - f(mesh->xend, j, k); + f(mesh->xend + 2, j, k) = f(mesh->xend + 1, j, k); + } + } + } + }; + applyBoundaryAverageZ(P); + applyBoundaryAverageZ(P_par); + applyBoundaryAverageZ(P_perp); + } + + // Calculate temperature + // Not using density boundary condition + N = getNoBoundary(species["density"]); + + T = P / floor(N, density_floor); + Field3D T_par = P_par / floor(N, density_floor); + Field3D T_perp = P_perp / floor(N, density_floor); + + set(species["pressure"], P); + set(species["pressure_perp"], P_perp); + set(species["pressure_par"], P_par); + set(species["temperature"], T); + set(species["temperature_perp"], T_perp); + set(species["temperature_par"], T_par); +} + +void EvolveAnisotropicPressure::finally(const Options& state) { + AUTO_TRACE(); + + /// Get the section containing this species + const auto& species = state["species"][name]; + + // Get updated pressure and temperature with boundary conditions + // Note: Retain pressures which fall below zero + P.clearParallelSlices(); + P.setBoundaryTo(get(species["pressure"])); + Field3D Pfloor = floor(P, 0.0); // Restricted to never go below zero + + T = get(species["temperature"]); + N = get(species["density"]); + + if (species.isSet("charge") and (fabs(get(species["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"]); + + ddt(E) = -Div_n_bxGrad_f_B_XPPM(E, phi, bndry_flux, poloidal_flows, true); + + ddt(PA) = -Div_n_bxGrad_f_B_XPPM(PA, phi, bndry_flux, poloidal_flows, true); + } else { + ddt(E) = 0.0; + ddt(PA) = 0.0; + } + + if (species.isSet("velocity")) { + Field3D V = get(species["velocity"]); + + // Typical wave speed used for numerical diffusion + Field3D fastest_wave; + if (state.isSet("fastest_wave")) { + fastest_wave = get(state["fastest_wave"]); + } else { + BoutReal AA = get(species["AA"]); + fastest_wave = sqrt(T / AA); + } + + Field3D flow_ylow; + ddt(E) -= FV::Div_par_mod(E + P_par, V, fastest_wave, flow_ylow); + + ddt(PA) -= FV::Div_par_mod(PA, V, fastest_wave, flow_ylow) + + (2. * P_par + P_perp) * Grad_par(V) - P_perp * Div_par(V); + } + + ////////////////////// + // Other sources + + 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)); + final_source = source * source_prefactor; + } else { + final_source = source; + } + + Sp = final_source; + if (species.isSet("energy_source")) { + Sp += (2. / 3) * get(species["energy_source"]); // For diagnostic output + } + + ddt(E) += (3. / 2) * Sp; +} + +void EvolveAnisotropicPressure::outputVars(Options& state) { + AUTO_TRACE(); + // Normalisations + auto Nnorm = get(state["Nnorm"]); + auto Tnorm = get(state["Tnorm"]); + auto Omega_ci = get(state["Omega_ci"]); + + BoutReal Pnorm = SI::qe * Tnorm * Nnorm; // Pressure normalisation + + state[std::string("E") + name].setAttributes( + {{"time_dimension", "t"}, + {"units", "J/m^3"}, + {"conversion", Pnorm}, + {"standard_name", "energy density"}, + {"long_name", name + " energy density"}, + {"species", name}, + {"source", "evolve_anisotropic_pressure"}}); + + state[std::string("PA") + name].setAttributes( + {{"time_dimension", "t"}, + {"units", "Pa"}, + {"conversion", Pnorm}, + {"standard_name", "pressure anisotropy"}, + {"long_name", name + " pressure anisotropy"}, + {"species", name}, + {"source", "evolve_anisotropic_pressure"}}); + + if (diagnose) { + set_with_attrs(state[std::string("T") + name], T, + {{"time_dimension", "t"}, + {"units", "eV"}, + {"conversion", Tnorm}, + {"standard_name", "temperature"}, + {"long_name", name + " temperature"}, + {"species", name}, + {"source", "evolve_anisotropic_pressure"}}); + + set_with_attrs(state[std::string("ddt(E") + name + std::string(")")], ddt(E), + {{"time_dimension", "t"}, + {"units", "J m^-3 s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", std::string("Rate of change of ") + name + " energy"}, + {"species", name}, + {"source", "evolve_anisotropic_pressure"}}); + + set_with_attrs(state[std::string("SP") + name], Sp, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"standard_name", "pressure source"}, + {"long_name", name + " pressure source"}, + {"species", name}, + {"source", "evolve_anisotropic_pressure"}}); + + set_with_attrs(state[std::string("P") + name + std::string("_src")], final_source, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"standard_name", "pressure source"}, + {"long_name", name + " pressure source"}, + {"species", name}, + {"source", "evolve_anisotropic_pressure"}}); + } +} diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index 76c1a3783..00d8c3c2e 100644 --- a/src/evolve_momentum.cxx +++ b/src/evolve_momentum.cxx @@ -1,12 +1,12 @@ +#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" @@ -18,7 +18,7 @@ EvolveMomentum::EvolveMomentum(std::string name, Options& alloptions, Solver* so readWrite("species:{name}:{outputs}")}), name(name) { AUTO_TRACE(); - + // Evolve the momentum in time solver->add(NV, std::string("NV") + name); @@ -26,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") @@ -40,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; @@ -82,14 +83,14 @@ 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 set(species["momentum"], NV); } -void EvolveMomentum::finally(const Options &state) { +void EvolveMomentum::finally(const Options& state) { AUTO_TRACE(); auto& species = state["species"][name]; @@ -137,17 +138,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); @@ -173,10 +173,15 @@ 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")) { + if (species.isSet("pressure_par")) { // Parallel pressure if available + Field3D P = get(species["pressure_par"]); + ddt(NV) -= Grad_par(P); + } else if (species.isSet("pressure")) { // Isotropic pressure Field3D P = get(species["pressure"]); ddt(NV) -= Grad_par(P); } @@ -186,7 +191,10 @@ void EvolveMomentum::finally(const Options &state) { const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); ddt(NV) -= Div_n_g_bxGrad_f_B_XZ(NV, V, -Apar_flutter); - if (species.isSet("pressure")) { + if (species.isSet("pressure_par")) { // Parallel pressure if available + Field3D P = get(species["pressure_par"]); + ddt(NV) -= bracket(P, Apar_flutter, BRACKET_ARAKAWA); + } else if (species.isSet("pressure")) { Field3D P = get(species["pressure"]); ddt(NV) -= bracket(P, Apar_flutter, BRACKET_ARAKAWA); } @@ -195,12 +203,13 @@ 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) { - ddt(NV) += - Div_Perp_Lap_FV_Index(density_floor / softFloor(N, 1e-3 * density_floor), NV); + ddt(NV) += Div_Perp_Lap_FV_Index(density_floor / softFloor(N, 1e-3 * density_floor), + NV); } if (low_p_diffuse_perp) { @@ -255,7 +264,7 @@ void EvolveMomentum::finally(const Options &state) { NV = NV_solver; } -void EvolveMomentum::outputVars(Options &state) { +void EvolveMomentum::outputVars(Options& state) { AUTO_TRACE(); // Normalisations auto Nnorm = get(state["Nnorm"]); @@ -302,24 +311,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"}}); } } }