diff --git a/CMakeLists.txt b/CMakeLists.txt index 9eae66871..2ed803704 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -209,6 +209,32 @@ configure_file( "${PROJECT_BINARY_DIR}/include/revision.hxx") + + +# Once built, copy the data and test directories +option(HERMES_COPY_EXAMPLES_TO_BUILD "Copy the examples to the build directory" + ON) + +option(HERMES_COPY_JSON_DATABASE_TO_BUILD + "Copy the json database to the build directory" ON) + +if(HERMES_COPY_EXAMPLES_TO_BUILD OR HERMES_COPY_JSON_DATABASE_TO_BUILD) + add_custom_target(setup-examples ALL) +else() + add_custom_target(setup-examples) +endif() + + +add_custom_command(TARGET setup-examples PRE_BUILD + COMMAND ${CMAKE_COMMAND} -E copy_directory + ${CMAKE_SOURCE_DIR}/examples $/examples) + +add_custom_command(TARGET setup-examples PRE_BUILD + COMMAND ${CMAKE_COMMAND} -E copy_directory + ${CMAKE_SOURCE_DIR}/json_database $/json_database) + + + # Tests option(HERMES_COPY_TESTS_TO_BUILD "Copy the tests to the build directory" ON) if(HERMES_COPY_TESTS_TO_BUILD) @@ -225,23 +251,17 @@ else() add_custom_target(setup-test) endif() -add_custom_target(setup-examples) - +# add_custom_target(setup-examples) -# copy the data and test directories -add_custom_command(TARGET setup-examples PRE_BUILD - COMMAND ${CMAKE_COMMAND} -E copy_directory - ${CMAKE_SOURCE_DIR}/examples $/examples) - -add_custom_command(TARGET setup-examples PRE_BUILD - COMMAND ${CMAKE_COMMAND} -E copy_directory - ${CMAKE_SOURCE_DIR}/json_database $/json_database) add_custom_command(TARGET setup-test PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/tests $/tests) + + + if(HERMES_TESTS) enable_testing() @@ -298,6 +318,8 @@ if(HERMES_TESTS) DEPENDS setup-test) endif() endfunction() + hermes_add_integrated_test(1D-Cimpurities-fci + REQUIRES BOUT_ENABLE_METRIC_3D) hermes_add_integrated_test(2D-vorticity-inversion-fci REQUIRES BOUT_ENABLE_METRIC_3D) hermes_add_integrated_test(1D-sheathtest-recycling-fci diff --git a/include/amjuel_reaction.hxx b/include/amjuel_reaction.hxx index c04cd37eb..9632581b0 100644 --- a/include/amjuel_reaction.hxx +++ b/include/amjuel_reaction.hxx @@ -102,16 +102,9 @@ protected: BoutReal avgNe = 0.0; BoutReal avgN1 = 0.0; BoutReal avgTe = 0.0; - if (Ne.isFci()) { - avgNe = (4.0 * Ne[i] + Ne.ydown()[iym] + Ne.yup()[iyp]) / 6.0; - avgN1 = (4.0 * N1[i] + N1.ydown()[iym] + N1.yup()[iyp]) / 6.0; - avgTe = (4.0 * Te[i] + Te.ydown()[iym] + Te.yup()[iyp]) / 6.0; - } else { - avgNe = Ne[i]; - avgN1 = N1[i]; - avgTe = Te[i]; - } - + avgNe = Ne[i]; + avgN1 = N1[i]; + avgTe = Te[i]; reaction_rate[i] = avgNe * avgN1 * evaluate(rate_coefs, avgTe * Tnorm, avgNe * Nnorm) * Nnorm / FreqNorm * rate_multiplier; } @@ -179,15 +172,11 @@ protected: BoutReal avgNe = 0.0; BoutReal avgN1 = 0.0; BoutReal avgTe = 0.0; - if (Ne.isFci()) { - avgNe = (4.0 * Ne[i] + Ne.ydown()[iym] + Ne.yup()[iyp]) / 6.0; - avgN1 = (4.0 * N1[i] + N1.ydown()[iym] + N1.yup()[iyp]) / 6.0; - avgTe = (4.0 * Te[i] + Te.ydown()[iym] + Te.yup()[iyp]) / 6.0; - } else { - avgNe = Ne[i]; - avgN1 = N1[i]; - avgTe = Te[i]; - } + + avgNe = Ne[i]; + avgN1 = N1[i]; + avgTe = Te[i]; + energy_loss[i] = avgNe * avgN1 * evaluate(radiation_coefs, avgTe * Tnorm, avgNe * Nnorm) * Nnorm / (Tnorm * FreqNorm) * radiation_multiplier; } diff --git a/include/electromagnetic.hxx b/include/electromagnetic.hxx index de0870a57..f8c2618bf 100644 --- a/include/electromagnetic.hxx +++ b/include/electromagnetic.hxx @@ -66,14 +66,14 @@ private: BoutReal beta_em; // Normalisation coefficient mu_0 e T n / B^2 std::unique_ptr aparSolver; // Laplacian solver in X-Z - + bool use_normdensity; bool const_gradient; // Set neumann boundaries by extrapolation BoutReal apar_boundary_timescale; // Relaxation timescale BoutReal last_time; // The last time the boundaries were updated bool magnetic_flutter; ///< Set the magnetic flutter term? Field3D Apar_flutter; - + Field3D zeroes; bool diagnose; ///< Output additional diagnostics? }; diff --git a/include/evolve_density.hxx b/include/evolve_density.hxx index 441a7bb67..e40fd7b59 100644 --- a/include/evolve_density.hxx +++ b/include/evolve_density.hxx @@ -64,7 +64,7 @@ private: BoutReal AA; ///< Atomic mass e.g. proton = 1 Field3D N; ///< Species density (normalised, evolving) - + bool magnetic_flutter; bool bndry_flux; ///< Allow flows through boundaries? bool exb_advection; ///< Include ExB advection? bool poloidal_flows; ///< Include ExB flow in Y direction? diff --git a/include/evolve_momentum.hxx b/include/evolve_momentum.hxx index 8f3df5278..20933b98a 100644 --- a/include/evolve_momentum.hxx +++ b/include/evolve_momentum.hxx @@ -50,7 +50,7 @@ private: BoutReal pressure_floor; bool low_p_diffuse_perp; ///< Cross-field diffusion at low pressure? BoutReal scale_ExB; - + bool magnetic_flutter; BoutReal hyper_z; ///< Hyper-diffusion BoutReal hyper_nv; std::string Vname; diff --git a/include/evolve_pressure.hxx b/include/evolve_pressure.hxx index 1cc49f078..82ff5b8b9 100644 --- a/include/evolve_pressure.hxx +++ b/include/evolve_pressure.hxx @@ -84,6 +84,8 @@ private: bool thermal_conduction; ///< Include thermal conduction? BoutReal kappa_coefficient; ///< Leading numerical coefficient in parallel heat flux calculation BoutReal kappa_limit_alpha; ///< Flux limit if >0 + bool kappa_limit_grillix; + BoutReal kappa_limit_Lpar; bool disable_ddt; bool p_div_v; ///< Use p*Div(v) form? False -> v * Grad(p) BoutReal T_lowsource; @@ -103,7 +105,7 @@ private: Field3D source, final_source; ///< External pressure source Field3D Sp; ///< Total pressure source FieldGeneratorPtr source_prefactor_function; - + bool magnetic_flutter; BoutReal adapt_source; BoutReal hyper_z; ///< Hyper-diffusion diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index 5eabc3793..48688bdcc 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -90,6 +90,9 @@ private: Field3D viscosity_core; Field3D viscosity_par; Field3D ones; + Field3D zeroes; + Field3D is_SOL; + bool core_dissipation; bool phi_dissipation; /// Parallel dissipation of potential BoutReal vort_timedissipation; diff --git a/include/vorticity.hxx b/include/vorticity.hxx index a70328d3b..dd7a432ad 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -120,12 +120,17 @@ private: BoutReal average_atomic_mass; // Weighted average atomic mass, for polarisaion current (Boussinesq approximation) bool poloidal_flows; ///< Include poloidal ExB flow? bool bndry_flux; ///< Allow flows through radial boundaries? - + bool boussinesq; bool collisional_friction; ///< Damping of vorticity due to collisional friction + std::unique_ptr phiSolver_zonalneumann; + + bool zonal_neumann; + bool sheath_boundary; ///< Set outer boundary to j=0? bool vort_dissipation; ///< Parallel dissipation of vorticity bool phi_dissipation; ///< Parallel dissipation of potential + BoutReal phi_diss_factor; bool phi_sheath_dissipation; ///< Dissipation at the sheath if phi < 0 bool damp_core_vorticity; ///< Damp axisymmetric component of vorticity @@ -133,7 +138,7 @@ private: BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised] BoutReal phi_boundary_last_update; ///< Time when last updated bool phi_core_averagey; ///< Average phi core boundary in Y? - + bool magnetic_flutter; bool split_n0; // Split phi into n=0 and n!=0 components LaplaceXY* laplacexy; // Laplacian solver in X-Y (n=0) Field3D logB; diff --git a/src/anomalous_diffusion_3d.cxx b/src/anomalous_diffusion_3d.cxx index 67306ab07..0d7e9c43a 100644 --- a/src/anomalous_diffusion_3d.cxx +++ b/src/anomalous_diffusion_3d.cxx @@ -97,7 +97,7 @@ void AnomalousDiffusion3D::transform(Options& state) { ? GET_NOBOUNDARY(Field3D, species["temperature"]) : 0.0; Field3D V = - species.isSet("velocity") ? GET_NOBOUNDARY(Field3D, species["velocity"]) : 0.0; + (species.isSet("velocity") && include_nu)? GET_NOBOUNDARY(Field3D, species["velocity"]) : 0.0; Field3D flow_xlow, flow_zlow; // Flows through cell faces diff --git a/src/electromagnetic.cxx b/src/electromagnetic.cxx index eeb1c87cb..30daca12f 100644 --- a/src/electromagnetic.cxx +++ b/src/electromagnetic.cxx @@ -38,6 +38,10 @@ Electromagnetic::Electromagnetic(std::string name, Options &alloptions, Solver* .doc("Extrapolate gradient of Apar into all radial boundaries?") .withDefault(false); + use_normdensity = options["use_normdensity"] + .doc("Use the normalized density instead of the full density in the Apar equation?") + .withDefault(false); + // Give Apar an initial value because we solve Apar by iteration // starting from the previous solution // Note: On restart the value is restored (if available) in restartVars @@ -73,6 +77,11 @@ Electromagnetic::Electromagnetic(std::string name, Options &alloptions, Solver* .doc("Output additional diagnostics?") .withDefault(false); + zeroes = 0.0; + zeroes.applyBoundary("neumann"); + bout::globals::mesh->communicate(zeroes); + zeroes.applyParallelBoundary("parallel_neumann_o1"); + magnetic_flutter = options["magnetic_flutter"] .doc("Set magnetic flutter terms (Apar_flutter)?") .withDefault(false); @@ -127,15 +136,21 @@ void Electromagnetic::transform(Options &state) { const BoutReal A = get(species["AA"]); // Coefficient in front of A_|| - alpha_em += floor(N, 1e-5) * (SQ(Z) / A); + if (use_normdensity) { + alpha_em += (SQ(Z) / A); + } else { + alpha_em += floor(N, 1e-5) * (SQ(Z) / A); + } + // Right hand side Ajpar += mom * (Z / A); } // Invert Helmholtz equation for Apar aparSolver->setCoefA((-beta_em) * alpha_em); - + aparSolver->setCoefC(1.0); + //aparSolver->setCoefA(0.0); if (const_gradient) { // Set gradient boundary condition from gradient inside boundary Field3D rhs = (-beta_em) * Ajpar; @@ -173,7 +188,7 @@ void Electromagnetic::transform(Options &state) { // Use previous value of Apar as initial guess Apar = aparSolver->solve(rhs, Apar); } else { - Apar = aparSolver->solve((-beta_em) * Ajpar, Apar); + Apar = aparSolver->solve((-beta_em) * Ajpar, zeroes); } // Save in the state @@ -194,17 +209,24 @@ void Electromagnetic::transform(Options &state) { const Field3D N = GET_NOBOUNDARY(Field3D, species["density"]); Field3D nv = getNonFinal(species["momentum"]); - nv -= Z * N * Apar; + if (use_normdensity) { + nv -= Z * Apar; + } else { + nv -= Z * N * Apar; + } // Note: velocity is momentum / (A * N) Field3D v = getNonFinal(species["velocity"]); - v -= (Z / A) * N * Apar / floor(N, 1e-5); - // Need to update the guard cells + if (use_normdensity) { + v -= (Z / A) * Apar / floor(N, 1e-5); + } else { + v -= (Z / A) * N * Apar / floor(N, 1e-5); + } + nv.applyBoundary("neumann"); v.applyBoundary("neumann"); bout::globals::mesh->communicate(nv, v); v.applyParallelBoundary("parallel_neumann_o1"); nv.applyParallelBoundary("parallel_neumann_o1"); - set(species["momentum"], nv); set(species["velocity"], v); @@ -212,32 +234,35 @@ void Electromagnetic::transform(Options &state) { if (magnetic_flutter) { // Magnetic flutter terms - Apar_flutter = Apar - DC(Apar); - + if (bout::globals::mesh->isFci()){ + Apar_flutter = Apar; + set(state["fields"]["Apar_flutter"], Apar); + } else { + Apar_flutter = Apar - DC(Apar); // Ensure that guard cells are communicated - Apar.getMesh()->communicate(Apar_flutter); + Apar.getMesh()->communicate(Apar_flutter); + set(state["fields"]["Apar_flutter"], Apar_flutter); - set(state["fields"]["Apar_flutter"], Apar_flutter); - -#if 0 - // Create a vector A from covariant components - // (A^x, A^y, A^z) - // Note: b = e_y / (JB) + #if 0 + // Create a vector A from covariant components + // (A^x, A^y, A^z) + // Note: b = e_y / (JB) const auto* coords = Apar.getCoordinates(); Vector3D A; A.covariant = true; A.x = A.z = 0.0; A.y = Apar_flutter * (coords->J * coords->Bxy); - // Perturbed magnetic field vector - // Note: Contravariant components (dB_x, dB_y, dB_z) + // Perturbed magnetic field vector + // Note: Contravariant components (dB_x, dB_y, dB_z) Vector3D delta_B = Curl(A); - // Set components of the perturbed unit vector - // Note: Options can't (yet) contain vectors + // Set components of the perturbed unit vector + // Note: Options can't (yet) contain vectors set(state["fields"]["deltab_flutter_x"], delta_B.x / coords->Bxy); set(state["fields"]["deltab_flutter_z"], delta_B.z / coords->Bxy); #endif + } } } diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index 40f413963..f6d0ecc0b 100644 --- a/src/evolve_density.cxx +++ b/src/evolve_density.cxx @@ -51,6 +51,7 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv .doc("Perpendicular diffusion at low density") .withDefault(false); + pressure_floor = density_floor * (1./get(alloptions["units"]["eV"])); @@ -131,6 +132,11 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv .withDefault(source) / source_normalisation; + magnetic_flutter= + n_options["magnetic_flutter"] + .doc("Use flutter terms??") + .withDefault(true); + disable_ddt = n_options["disable_ddt"] .withDefault(false); @@ -299,12 +305,13 @@ void EvolveDensity::finally(const Options& state) { flow_ylow = 0.0; ddt(N) -= FV::Div_par_mod(N, V, fastest_wave, flow_ylow, false, dissipative); - if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter") and magnetic_flutter) { // Magnetic flutter term const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); // Note: Using -Apar_flutter rather than reversing sign in front, // so that upwinding is handled correctly - ddt(N) -= Div_n_g_bxGrad_f_B_XZ(N, V, -Apar_flutter); + auto* coord = N.getCoordinates(); + ddt(N) -= coord->Bxy * bracket(N*V/coord->Bxy, Apar_flutter, BRACKET_ARAKAWA) * bracket_factor; } } diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index 63325998f..ffec765ee 100644 --- a/src/evolve_momentum.cxx +++ b/src/evolve_momentum.cxx @@ -105,6 +105,11 @@ EvolveMomentum::EvolveMomentum(std::string name, Options &alloptions, Solver *so disable_ddt = nv_options["disable_ddt"] .withDefault(false); + magnetic_flutter= + nv_options["magnetic_flutter"] + .doc("Use flutter terms??") + .withDefault(true); + } void EvolveMomentum::transform(Options &state) { @@ -222,7 +227,7 @@ void EvolveMomentum::finally(const Options &state) { } ddt(NV) += Z * Apar * dndt; } - if (state["fields"].isSet("Apar_flutter")) { + if (state["fields"].isSet("Apar_flutter") and magnetic_flutter) { // Magnetic flutter term const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); @@ -254,14 +259,14 @@ void EvolveMomentum::finally(const Options &state) { ddt(NV) -= Grad_par(P); } - if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter") and magnetic_flutter) { // Magnetic flutter term 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")) { Field3D P = get(species["pressure"]); - ddt(NV) -= bracket(P, Apar_flutter, BRACKET_ARAKAWA); + ddt(NV) -= bracket(P, Apar_flutter, BRACKET_ARAKAWA) * bracket_factor; } } diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index 45b7df0d0..8ab6d238d 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -168,6 +168,11 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so } } + magnetic_flutter= + p_options["magnetic_flutter"] + .doc("Use flutter terms??") + .withDefault(true); + neumann_boundary_average_z = p_options["neumann_boundary_average_z"] .doc("Apply neumann boundary with Z average?") .withDefault(false); @@ -213,6 +218,15 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so .doc("Flux limiter factor. < 0 means no limit. Typical is 0.2 for electrons, 1 for ions.") .withDefault(-1.0); + kappa_limit_grillix = options["kappa_limit_grillix"] + .doc("Use grillix-style flux limiter?") + .withDefault(false); + + kappa_limit_Lpar = options["kappa_limit_Lpar"] + .doc("Parallel decay length, given by T / grad_par Te") + .withDefault(1.0) / Lnorm; + + if (mesh->isFci()) { const auto coord = mesh->getCoordinates(); // Note: This is 1 for a Clebsch coordinate system @@ -349,11 +363,11 @@ void EvolvePressure::finally(const Options& state) { flow_ylow *= 5. / 2; // Energy flow } - if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter") and magnetic_flutter) { // Magnetic flutter term const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); ddt(P) -= (5. / 3) * Div_n_g_bxGrad_f_B_XZ(P, V, -Apar_flutter); - ddt(P) += (2. / 3) * V * bracket(P, Apar_flutter, BRACKET_ARAKAWA); + ddt(P) += (2. / 3) * V * bracket(P, Apar_flutter, BRACKET_ARAKAWA) * bracket_factor; } if (numerical_viscous_heating || diagnose) { @@ -414,7 +428,15 @@ void EvolvePressure::finally(const Options& state) { // Note: Coefficient is slightly different for electrons (3.16) and ions (3.9) kappa_par = kappa_coefficient * Pfloor * tau / AA; - if (kappa_limit_alpha > 0.0) { + if (kappa_limit_alpha > 0.0 && kappa_limit_grillix) { + + // Based on https://iopscience.iop.org/article/10.1088/1741-4326/ac1e61 + // Equation 11 + // kappa_limit_Lpar = q_95 * R0 + Field3D denom = 1.0 + kappa_par / (kappa_limit_alpha * sqrt(T / AA) * N * kappa_limit_Lpar); + kappa_par /= denom; + + } else if (kappa_limit_alpha > 0.0) { /* * Flux limiter, as used in SOLPS. * @@ -464,7 +486,7 @@ void EvolvePressure::finally(const Options& state) { } } - if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter") and magnetic_flutter) { // Magnetic flutter term. The operator splits into 4 pieces: // Div(k b b.Grad(T)) = Div(k b0 b0.Grad(T)) + Div(k d0 db.Grad(T)) // + Div(k db b0.Grad(T)) + Div(k db db.Grad(T)) @@ -473,9 +495,9 @@ void EvolvePressure::finally(const Options& state) { const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); Field3D db_dot_T = bracket(T, Apar_flutter, BRACKET_ARAKAWA); Field3D b0_dot_T = Grad_par(T); - mesh->communicate(db_dot_T, b0_dot_T); db_dot_T.applyBoundary("neumann"); b0_dot_T.applyBoundary("neumann"); + mesh->communicate(db_dot_T, b0_dot_T); ddt(P) += (2. / 3) * (Div_par(kappa_par * db_dot_T) - Div_n_g_bxGrad_f_B_XZ(kappa_par, db_dot_T + b0_dot_T, Apar_flutter)); } diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index c5ef62ef0..24482d07c 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -161,6 +161,20 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so ones.applyBoundary("neumann"); mesh->communicate(ones); ones.applyParallelBoundary("parallel_neumann_o1"); + + zeroes = 0.0; + zeroes.applyBoundary("neumann"); + mesh->communicate(zeroes); + zeroes.applyParallelBoundary("parallel_neumann_o1"); + + core_dissipation = options["core_dissipation"] + .doc("Use core dissipation of vorticity? Grid field required!") + .withDefault(false); + + if (core_dissipation || (vort_timedissipation > 0.0)) { + mesh->get(is_SOL, "is_SOL", 0.0); + } + } @@ -366,8 +380,15 @@ void RelaxPotential::finally(const Options& state) { // Solve diffusion equation for potential if (vort_timedissipation > 0.0) { - ddt(Vort) -= vort_timedissipation * Vort; + ddt(Vort) -= (1.0 - is_SOL) * Vort / vort_timedissipation; } + + if (core_dissipation) { + Field3D sound_speed = get(state["sound_speed"]); + Field3D dummy; + ddt(Vort) -= (1.0 - is_SOL) * FV::Div_par_mod(Vort, zeroes, sound_speed, dummy); + } + if (boussinesq) { @@ -397,8 +418,8 @@ void RelaxPotential::finally(const Options& state) { } } else { // Non-Boussinesq. Calculate mass density by summing over species - throw BoutException("Non_boussinesq not implemented"); // Calculate vorticity from potential phi + Field3D flow_xlow_phi,flow_zlow_phi; Field3D phi_vort = 0.0; for (auto& kv : allspecies.getChildren()) { const Options& species = kv.second; @@ -413,12 +434,12 @@ void RelaxPotential::finally(const Options& state) { const BoutReal Ai = get(species["AA"]); const Field3D Ni = get(species["density"]); - phi_vort += Div_a_Grad_perp((Ai / Bsq) * Ni, phi); + phi_vort += (*dagp)((Ai / Bsq) * Ni, phi, flow_xlow_phi, flow_zlow_phi, false); if (diamagnetic_polarisation and species.isSet("pressure")) { // Calculate the diamagnetic flow contribution const Field3D Pi = get(species["pressure"]); - phi_vort += Div_a_Grad_perp(Ai / Bsq, Pi); + phi_vort += (*dagp)(Ai / Bsq / Zi, Pi, flow_xlow_phi, flow_zlow_phi, false); } } diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 10f02f332..069e0980d 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -1,6 +1,9 @@ #include "../include/vorticity.hxx" #include "../include/div_ops.hxx" +#include "../include/hermes_utils.hxx" +#include "../include/hermes_build_config.hxx" + #include #include @@ -14,11 +17,6 @@ using bout::globals::mesh; namespace { -BoutReal floor(BoutReal value, BoutReal min) { - if (value < min) - return min; - return value; -} Ind3D indexAt(const Field3D& f, int x, int y, int z) { int ny = f.getNy(); @@ -51,6 +49,10 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { output_ddt = options["output_ddt"] .doc("Include ExB advection?") .withDefault(false); + + boussinesq = options["boussinesq"] + .doc("Use Boussinesq approximation?") + .withDefault(true); diamagnetic = options["diamagnetic"].doc("Include diamagnetic current?").withDefault(true); @@ -64,6 +66,11 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { .doc("Include diamagnetic drift in polarisation current?") .withDefault(true); + magnetic_flutter= + options["magnetic_flutter"] + .doc("Use flutter terms??") + .withDefault(true); + diamagnetic_bracketform = options["diamagnetic_bracketform"] .doc("Include diamagnetic form that uses arakawa brackets? FCI version") .withDefault(mesh->isFci()); @@ -146,6 +153,10 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { .doc("Relax x boundaries of phi towards Neumann?") .withDefault(false); + phi_diss_factor = options["phi_diss_factor"] + .doc("Prefactor to increase phi dissipation") + .withDefault(1.0); + phi_sheath_dissipation = options["phi_sheath_dissipation"] .doc("Add dissipation when phi < 0.0 at the sheath") .withDefault(false); @@ -176,6 +187,15 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { // Set coefficients for Boussinesq solve phiSolver->setCoefC(average_atomic_mass / SQ(coord->Bxy)); + zonal_neumann = options["zonal_neumann"] + .doc("Do a second laplace solve for the zonal neumann?") + .withDefault(false); + + if (zonal_neumann) { + phiSolver_zonalneumann = Laplacian::create(&options["laplacian_zonalneumann"]); + } + + if (phi_boundary_relax) { // Set the last update time to -1, so it will reset // the first time RHS function is called @@ -276,131 +296,335 @@ void Vorticity::transform(Options& state) { Vort.applyParallelBoundary(); + auto coord = mesh->getCoordinates(); + + + Field3D phi_plus_pi = 0.0; + + if (boussinesq) { - // Set the boundary of phi. Both 2D and 3D fields are kept, though the 3D field - // is constant in Z. This is for efficiency, to reduce the number of conversions. - // Note: For now the boundary values are all at the midpoint, - // and only phi is considered, not phi + Pi which is handled in Boussinesq solves - Pi_hat = 0.0; // Contribution from ion pressure, weighted by atomic mass / charge - if (diamagnetic_polarisation) { - // Diamagnetic term in vorticity. Note this is weighted by the mass - // This includes all species, including electrons - Options& allspecies = state["species"]; - for (auto& kv : allspecies.getChildren()) { - Options& species = allspecies[kv.first]; // Note: need non-const + phiSolver->setCoefC(average_atomic_mass / SQ(coord->Bxy)); + phiSolver->setCoefA(0.0); + // Set the boundary of phi. Both 2D and 3D fields are kept, though the 3D field + // is constant in Z. This is for efficiency, to reduce the number of conversions. + // Note: For now the boundary values are all at the midpoint, + // and only phi is considered, not phi + Pi which is handled in Boussinesq solves + Pi_hat = 0.0; // Contribution from ion pressure, weighted by atomic mass / charge + if (diamagnetic_polarisation) { + // Diamagnetic term in vorticity. Note this is weighted by the mass + // This includes all species, including electrons + Options& allspecies = state["species"]; + for (auto& kv : allspecies.getChildren()) { + Options& species = allspecies[kv.first]; // Note: need non-const - if (!(IS_SET_NOBOUNDARY(species["pressure"]) and species.isSet("charge") - and species.isSet("AA"))) { - continue; // No pressure, charge or mass -> no polarisation current + if (!(IS_SET_NOBOUNDARY(species["pressure"]) and species.isSet("charge") + and species.isSet("AA"))) { + continue; // No pressure, charge or mass -> no polarisation current + } + + const auto charge = get(species["charge"]); + if (fabs(charge) < 1e-5) { + // No charge + continue; + } + + // Don't need sheath boundary + const auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); + const auto AA = get(species["AA"]); + + Pi_hat += P * (AA / average_atomic_mass / charge); } + } - const auto charge = get(species["charge"]); - if (fabs(charge) < 1e-5) { - // No charge - continue; + Pi_hat.applyBoundary("neumann"); + + if (phi_boundary_relax) { + // Update the boundary regions by relaxing towards zero gradient + // on a given timescale. + + BoutReal time = get(state["time"]); + + if (phi_boundary_last_update < 0.0) { + // First time this has been called. + phi_boundary_last_update = time; + + } else if (time > phi_boundary_last_update) { + // Only update if time has advanced + // Uses an exponential decay of the weighting of the value in the boundary + // so that the solution is well behaved for arbitrary steps + BoutReal weight = exp(-(time - phi_boundary_last_update) / phi_boundary_timescale); + phi_boundary_last_update = time; + + if (mesh->firstX()) { + BoutReal phivalue = 0.0; + if (phi_core_averagey) { + BoutReal philocal = 0.0; + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + philocal += phi(mesh->xstart, j, k); + } + } + 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); + phivalue /= (np * mesh->LocalNz * mesh->LocalNy); + } + + for (int j = mesh->ystart; j <= mesh->yend; j++) { + if (!phi_core_averagey) { + phivalue = 0.0; // Calculate phi boundary for each Y index separately + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xstart, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + } + + // Old value of phi at boundary + BoutReal oldvalue = + 0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xstart - 1, j, k) = 2. * newvalue - phi(mesh->xstart, j, 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) + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal phivalue = 0.0; + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xend, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + + // Old value of phi at boundary + BoutReal oldvalue = 0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xend + 1, j, k) = 2. * newvalue - phi(mesh->xend, j, 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) + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } + } } + } else { + // phi_boundary_relax = false + // + // Set boundary from temperature, to be consistent with j=0 at sheath - // Don't need sheath boundary - const auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); - const auto AA = get(species["AA"]); + // Sheath multiplier Te -> phi (2.84522 for Deuterium) + BoutReal sheathmult = 0.0; + if (sheath_boundary) { + BoutReal Me_Mp = get(state["species"]["e"]["AA"]); + sheathmult = log(0.5 * sqrt(1. / (Me_Mp * PI))); + } - Pi_hat += P * (AA / average_atomic_mass / charge); + Field3D Te; // Electron temperature, use for outer boundary conditions + if (state["species"]["e"].isSet("temperature")) { + // Electron temperature set + Te = GET_NOBOUNDARY(Field3D, state["species"]["e"]["temperature"]); + } else { + Te = 0.0; + } + + // Sheath multiplier Te -> phi (2.84522 for Deuterium if Ti = 0) + if ( mesh->firstX() && !mesh->isFci()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal teavg = 0.0; // Average Te in Z + + for (int k = 0; k < mesh->LocalNz; k++) { + teavg += Te(mesh->xstart, j, k); + } + teavg /= mesh->LocalNz; + BoutReal phivalue = sheathmult * teavg; + + // Set midpoint (boundary) value + for (int k = 0; k < mesh->LocalNz; k++) { + BoutReal phivalue = sheathmult * teavg; + phi(mesh->xstart - 1, j, k) = 2. * phivalue - phi(mesh->xstart, j, 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) + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } + } + } + + if ( mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + // Set midpoint (boundary) value + for (int k = 0; k < mesh->LocalNz; k++) { + BoutReal phivalue = sheathmult * 0.5 * (Te(mesh->xend, j, k) + Te(mesh->xend+1, j, k)); + phi(mesh->xend + 1, j, k) = 2. * phivalue - phi(mesh->xend, j, 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) + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } + } } - } + phi.name = "phi"; - Pi_hat.applyBoundary("neumann"); + // Update boundary conditions. Two issues: + // 1) Solving here for phi + Pi, and then subtracting Pi from the result + // The boundary values should therefore include Pi + // 2) The INVERT_SET flag takes the value in the guard (boundary) cell + // and sets the boundary between cells to this value. + // This shift by 1/2 grid cell is important. - if (phi_boundary_relax) { - // Update the boundary regions by relaxing towards zero gradient - // on a given timescale. - - BoutReal time = get(state["time"]); - - if (phi_boundary_last_update < 0.0) { - // First time this has been called. - phi_boundary_last_update = time; - - } else if (time > phi_boundary_last_update) { - // Only update if time has advanced - // Uses an exponential decay of the weighting of the value in the boundary - // so that the solution is well behaved for arbitrary steps - BoutReal weight = exp(-(time - phi_boundary_last_update) / phi_boundary_timescale); - phi_boundary_last_update = time; - - if (mesh->firstX()) { - BoutReal phivalue = 0.0; - if (phi_core_averagey) { - BoutReal philocal = 0.0; - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - philocal += phi(mesh->xstart, j, k); - } - } - 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); - phivalue /= (np * mesh->LocalNz * mesh->LocalNy); - } + phi_plus_pi = phi + Pi_hat; - for (int j = mesh->ystart; j <= mesh->yend; j++) { - if (!phi_core_averagey) { - phivalue = 0.0; // Calculate phi boundary for each Y index separately - for (int k = 0; k < mesh->LocalNz; k++) { - phivalue += phi(mesh->xstart, j, k); - } - phivalue /= mesh->LocalNz; // Average in Z of point next to boundary - } + if ( mesh->firstX() ) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Average phi + Pi at the boundary, and set the boundary cell + // to this value. The phi solver will then put the value back + // onto the cell mid-point + phi_plus_pi(mesh->xstart - 1, j, k) = + 0.5 * (phi_plus_pi(mesh->xstart - 1, j, k) + phi_plus_pi(mesh->xstart, j, k)); + phi_plus_pi(mesh->xstart - 2, j, k) = phi_plus_pi(mesh->xstart - 1, j, k); + } + } + } - // Old value of phi at boundary - BoutReal oldvalue = - 0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0)); + if ( mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + phi_plus_pi(mesh->xend + 1, j, k) = + 0.5 * (phi_plus_pi(mesh->xend + 1, j, k) + phi_plus_pi(mesh->xend, j, k)); + phi_plus_pi(mesh->xend + 2, j, k) = phi_plus_pi(mesh->xend + 1, j, k); + } + } + } - // New value of phi at boundary, relaxing towards phivalue - BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + // Calculate potential + if (split_n0) { + //////////////////////////////////////////// + // Split into axisymmetric and non-axisymmetric components + Field2D Vort2D = DC(Vort); // n=0 component + Field2D phi_plus_pi_2d = DC(phi_plus_pi); + phi_plus_pi -= phi_plus_pi_2d; - // Set phi at the boundary to this value - for (int k = 0; k < mesh->LocalNz; k++) { - phi(mesh->xstart - 1, j, k) = 2. * newvalue - phi(mesh->xstart, j, k); + phi_plus_pi_2d = laplacexy->solve(Vort2D, phi_plus_pi_2d); - // 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) - phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); - } - } + // Solve non-axisymmetric part using X-Z solver + phi = phi_plus_pi_2d + + phiSolver->solve((Vort - Vort2D) * (Bsq / average_atomic_mass), phi_plus_pi) + - Pi_hat; + + } else { + const auto tosolve = Vort * (Bsq / average_atomic_mass); + checkData(tosolve); + checkData(phi_plus_pi); + try { + phi = phiSolver->solve(tosolve, phi_plus_pi) - Pi_hat; + } catch (const BoutException& e) { + Options debug; + debug["tosolve"] = tosolve; + debug["guess"] = phi_plus_pi; + debug["Vort"] = Vort; + debug["Bsq"] = Bsq; + debug["Pi_hat"] = Pi_hat; + mesh->outputVars(debug); + const std::string outname = + fmt::format("{}/BOUT.debug_vorticity.{}.nc", + Options::root()["datadir"].withDefault("data"), + BoutComm::rank()); + + bout::OptionsIO::create(outname)->write(debug); + MPI_Barrier(BoutComm::get()); + throw e; } + } - if (mesh->lastX()) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - BoutReal phivalue = 0.0; - for (int k = 0; k < mesh->LocalNz; k++) { - phivalue += phi(mesh->xend, j, k); - } - phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + if (zonal_neumann) { + Field2D avg_phi = DC(phi); + if ( mesh->firstX() ) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + phi_plus_pi(mesh->xstart - 1, j, k) = 0.5 * ( avg_phi(mesh->xstart - 1 ,j) + avg_phi(mesh->xstart ,j) ) + + 0.5 * (Pi_hat(mesh->xstart - 1, j, k) + Pi_hat(mesh->xstart, j, k)); + phi_plus_pi(mesh->xstart - 2, j, k) = phi_plus_pi(mesh->xstart - 1, j, k); + } + } + } + phiSolver_zonalneumann->setCoefC(average_atomic_mass / SQ(coord->Bxy)); + const auto tosolve = Vort * (Bsq / average_atomic_mass); + phi = phiSolver_zonalneumann->solve(tosolve, phi_plus_pi) - Pi_hat; + } - // Old value of phi at boundary - BoutReal oldvalue = 0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0)); - // New value of phi at boundary, relaxing towards phivalue - BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Set phi at the boundary to this value - for (int k = 0; k < mesh->LocalNz; k++) { - phi(mesh->xend + 1, j, k) = 2. * newvalue - phi(mesh->xend, j, k); + } else { // Non-boussinesq part - // 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) - phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); - } - } + ////////////////////////////////////////////////////////////////// + + Field3D AN_Bsq = 0.0; + Field3D RHS = 0.0; + Options& allspecies = state["species"]; + for (auto& kv : allspecies.getChildren()) { + Options& species = allspecies[kv.first]; // Note: need non-const + + if (!(IS_SET_NOBOUNDARY(species["pressure"]) and species.isSet("charge") + and species.isSet("AA") and IS_SET_NOBOUNDARY(species["density"]) )) { + continue; // No pressure, charge or mass -> no polarisation current + } + + const auto charge = get(species["charge"]); + if (fabs(charge) < 1e-5) { + // No charge + continue; + } + + if (charge < 0.0) { + // No electrons + continue; } + + // Don't need sheath boundary + const auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); + const auto N = GET_NOBOUNDARY(Field3D, species["density"]); + const auto AA = get(species["AA"]); + + AN_Bsq += AA*N / Bsq; + + Field3D dummy1; + Field3D dummy2; + RHS += (*dagp)(AA / Bsq / charge, P, dummy1, dummy2, false); + } - } else { + RHS.applyBoundary("free_o2"); + mesh->communicate(RHS); + + ////////////////////////////////////////////////////////////////// + + // phi_boundary_relax = false // // Set boundary from temperature, to be consistent with j=0 at sheath @@ -421,7 +645,7 @@ void Vorticity::transform(Options& state) { } // Sheath multiplier Te -> phi (2.84522 for Deuterium if Ti = 0) - if ( mesh->firstX()) { + if ( mesh->firstX() && !mesh->isFci()) { for (int j = mesh->ystart; j <= mesh->yend; j++) { BoutReal teavg = 0.0; // Average Te in Z @@ -445,99 +669,79 @@ void Vorticity::transform(Options& state) { if ( mesh->lastX()) { for (int j = mesh->ystart; j <= mesh->yend; j++) { - BoutReal teavg = 0.0; // Average Te in Z - - for (int k = 0; k < mesh->LocalNz; k++) { - teavg += Te(mesh->xend, j, k); - } - teavg /= mesh->LocalNz; - BoutReal phivalue = sheathmult * teavg; // Set midpoint (boundary) value for (int k = 0; k < mesh->LocalNz; k++) { + BoutReal phivalue = sheathmult * 0.5 * (Te(mesh->xend, j, k) + Te(mesh->xend+1, j, k)); phi(mesh->xend + 1, j, k) = 2. * phivalue - phi(mesh->xend, j, 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) phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } } } - } - phi.name = "phi"; - // Update boundary conditions. Two issues: - // 1) Solving here for phi + Pi, and then subtracting Pi from the result - // The boundary values should therefore include Pi - // 2) The INVERT_SET flag takes the value in the guard (boundary) cell - // and sets the boundary between cells to this value. - // This shift by 1/2 grid cell is important. - - Field3D phi_plus_pi = phi + Pi_hat; - - if ( mesh->firstX() ) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - // Average phi + Pi at the boundary, and set the boundary cell - // to this value. The phi solver will then put the value back - // onto the cell mid-point - phi_plus_pi(mesh->xstart - 1, j, k) = - 0.5 * (phi_plus_pi(mesh->xstart - 1, j, k) + phi_plus_pi(mesh->xstart, j, k)); - phi_plus_pi(mesh->xstart - 2, j, k) = phi_plus_pi(mesh->xstart - 1, j, k); + phi.name = "phi"; + + phiSolver->setCoefC(AN_Bsq); + + Field3D tosolve = (Vort-RHS)/AN_Bsq; + + phi = phiSolver->solve(tosolve, phi) ; + + + if (zonal_neumann) { + phiSolver_zonalneumann->setCoefC(AN_Bsq); + Field2D avg_phi = DC(phi); + if ( mesh->firstX() ) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xstart - 1, j, k) = 0.5 * ( avg_phi(mesh->xstart - 1 ,j) + avg_phi(mesh->xstart ,j) ) ; + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } + } } - } - } - if ( mesh->lastX()) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - phi_plus_pi(mesh->xend + 1, j, k) = - 0.5 * (phi_plus_pi(mesh->xend + 1, j, k) + phi_plus_pi(mesh->xend, j, k)); - phi_plus_pi(mesh->xend + 2, j, k) = phi_plus_pi(mesh->xend + 1, j, k); + if ( mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xend + 1, j, k) = 0.5 * (phi(mesh->xend + 1, j, k) + phi(mesh->xend, j, k)); + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } } + + + + phi = phiSolver_zonalneumann->solve(tosolve, phi); } - } - // Calculate potential - if (split_n0) { - //////////////////////////////////////////// - // Split into axisymmetric and non-axisymmetric components - Field2D Vort2D = DC(Vort); // n=0 component - Field2D phi_plus_pi_2d = DC(phi_plus_pi); - phi_plus_pi -= phi_plus_pi_2d; - phi_plus_pi_2d = laplacexy->solve(Vort2D, phi_plus_pi_2d); + + } // END non-boussinesq - // Solve non-axisymmetric part using X-Z solver - phi = phi_plus_pi_2d - + phiSolver->solve((Vort - Vort2D) * (Bsq / average_atomic_mass), phi_plus_pi) - - Pi_hat; + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + - } else { - const auto tosolve = Vort * (Bsq / average_atomic_mass); - checkData(tosolve); - checkData(phi_plus_pi); - try { - phi = phiSolver->solve(tosolve, phi_plus_pi) - Pi_hat; - } catch (const BoutException& e) { - Options debug; - debug["tosolve"] = tosolve; - debug["guess"] = phi_plus_pi; - debug["Vort"] = Vort; - debug["Bsq"] = Bsq; - debug["Pi_hat"] = Pi_hat; - mesh->outputVars(debug); - const std::string outname = - fmt::format("{}/BOUT.debug_vorticity.{}.nc", - Options::root()["datadir"].withDefault("data"), - BoutComm::rank()); - - bout::OptionsIO::create(outname)->write(debug); - MPI_Barrier(BoutComm::get()); - throw e; + if (zonal_neumann && boussinesq) { + Field2D avg_phi = DC(phi); + if ( mesh->firstX() ) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + phi_plus_pi(mesh->xstart - 1, j, k) = 0.5 * ( avg_phi(mesh->xstart - 1 ,j) + avg_phi(mesh->xstart ,j) ) + + 0.5 * (Pi_hat(mesh->xstart - 1, j, k) + Pi_hat(mesh->xstart, j, k)); + phi_plus_pi(mesh->xstart - 2, j, k) = phi_plus_pi(mesh->xstart - 1, j, k); + } + } } + const auto tosolve = Vort * (Bsq / average_atomic_mass); + phi = phiSolver_zonalneumann->solve(tosolve, phi_plus_pi) - Pi_hat; } - + // Ensure that potential is set in the communication guard cells mesh->communicate(phi); @@ -774,26 +978,43 @@ void Vorticity::finally(const Options& state) { const Field3D N = get(species["density"]); const Field3D NV = get(species["momentum"]); + const Field3D V = get(species["velocity"]); const BoutReal A = get(species["AA"]); - // Note: Using NV rather than N*V so that the cell boundary flux is correct - const Field3D jpar = (Z / A) * NV; - ddt(Vort) += Div_par(jpar); + Field3D fastest_wave; + if (state.isSet("fastest_wave")) { + fastest_wave = get(state["fastest_wave"]); + } else { + Field3D T = get(species["temperature"]); + BoutReal AA = get(species["AA"]); + fastest_wave = sqrt(T / AA); + } - if (state["fields"].isSet("Apar_flutter")) { + Field3D flow_ylow = 0.0; + ddt(Vort) += Z * FV::Div_par_mod(N, V, fastest_wave, flow_ylow, false, + false, true); + + if (state["fields"].isSet("Apar_flutter") && magnetic_flutter) { // Magnetic flutter term const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); // Div_par(jpar) = B * Grad_par(jpar / B) // Using the approximation for small delta-B/B // b dot Grad(jpar) = Grad_par(jpar) + [jpar, Apar] - ddt(Vort) += coord->Bxy * bracket(jpar / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA); + ddt(Vort) += coord->Bxy * bracket((Z/A)*NV / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA)*bracket_factor; } } // Viscosity if (has_viscosity) { - ddt(Vort) += Div_a_Grad_perp(viscosity, Vort); + if (Vort.isFci()) { + Field3D dummy1; + Field3D dummy2; + + ddt(Vort) += (*dagp)(viscosity, Vort, dummy1, dummy2, false); + } else { + ddt(Vort) += Div_a_Grad_perp(viscosity, Vort); + } } Field3D dummy; @@ -812,7 +1033,8 @@ void Vorticity::finally(const Options& state) { // Adds dissipation term like in other equations, but depending on gradient of // potential Field3D sound_speed = get(state["sound_speed"]); - ddt(Vort) -= FV::Div_par(-phi, 0.0, sound_speed); + Field3D dummy; + ddt(Vort) -= phi_diss_factor * FV::Div_par_mod(-phi, 0.0, sound_speed, dummy, true, false); } if (hyper > 0) { diff --git a/tests/integrated/1D-Cimpurities-fci/data/BOUT.inp b/tests/integrated/1D-Cimpurities-fci/data/BOUT.inp new file mode 100644 index 000000000..978d9b8cd --- /dev/null +++ b/tests/integrated/1D-Cimpurities-fci/data/BOUT.inp @@ -0,0 +1,320 @@ + +nout = 100 +timestep = 1e6 / nout + +[mesh] + +file = "../../../../grids/straight_slab_120m/MMS_straight_slab_8_128_4_True.fci.grid.nc" + +extrapolate_y = false + +[mesh:paralleltransform] +type = fci + +[input] + +error_on_unused_options = false + +[solver] +mxstep = 10000 +rtol = 1e-6 +atol = 1e-12 +type = cvode + + +[hermes] +components = h+, e, c, c+, c+2, c+3, c+4,c+5, c+6, collisions, sound_speed, sheath_boundary_parallel, reactions, ion_viscosity, recycling_fci + + +Nnorm = 1e18 +Bnorm = 1 +Tnorm = 5 + +[collisions] +electron_ion = true +ion_ion = true +neutral_neutral = false +ion_neutral = false +electron_neutral = true + +[recycling_fci] +species = h+ + + +[h+] # Ions +type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion_3d + +charge = 1.0 +AA = 1.0 + +diagnose = true + +thermal_conduction = true # in evolve_pressure +p_div_v = true + +dissipative = false + +recycle_as = c +target_recycle_multiplier = 0.01 +target_recycle = true + +anomalous_D_par = 1e4 + +[Nh+] +function = (1e19 / `hermes`:Nnorm) + +flux = 1e23 + +adapt_source = 2e19 + + +source = flux * gauss(y - pi, 0.25) + +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + +[NVh+] +function = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + +[Ph+] +function = (10.0 / `hermes`:Tnorm) * `Nh+`:function + +adapt_source = 50.0 + +heating = 1e7 +T_source = heating / (`Nh+`:flux * 1.602e-19 * 1.5) +source = `Nh+`:source * T_source * 1.602176634e-19 + + + +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + + +[e] # electrons +type = quasineutral, zero_current, evolve_pressure + +charge = -1.0 +AA = 1.0 / 1836.0 + +thermal_conduction = true # in evolve_pressure +p_div_v = true + +dissipative = false + +[Pe] + +adapt_source = 50.0 + +function = `Ph+`:function +source = `Ph+`:source +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + + + +################################################################################################################## + +[reactions] +type = ( + c + e -> c+ + 2e, + c+ + e -> c+2 + 2e, + c+2 + e -> c+3 + 2e, + c+3 + e -> c+4 + 2e, + c+4 + e -> c+5 + 2e, + c+5 + e -> c+6 + 2e + ) + + + + + +[c] +type = neutral_mixed, anomalous_diffusion_3d +AA = 6.0 +anomalous_D_par = 1e4 +evolve_momentum = true +neutral_conduction = true +neutral_viscosity = true +neutral_lmax = 0.05 +precondition = true +dissipative = false +density_floor = 1e-6 +parallel_dirichlet = false +diagnose = true +n_lowsource = 1e15 +[Nc] +function = 1e-6 +source = 0.0 +bndry_all = neumann # All other boundaries low density +bndry_par_all = parallel_neumann_o1 +bndry_xout = neumann +bndry_xin = neumann +[NVc] +bndry_all = neumann +bndry_xin = neumann +bndry_par_all = parallel_neumann_o1 +[Pc] +Teprofile = 10 / `hermes:Tnorm` +function = `Nc`:function * Teprofile +source = 0.0 +bndry_par_all = parallel_neumann_o1 +bndry_all = neumann + + +[c+] # Ions +type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion_3d +charge = 1.0 +AA = 6.0 +thermal_conduction = true # in evolve_pressure +p_div_v = true +dissipative = false +anomalous_D_par = 1e4 +n_lowsource = 1e15 +[Nc+] +function = (1e16 / `hermes`:Nnorm) +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[NVc+] +function = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[Pc+] +function = (10.0 / `hermes`:Tnorm) * `Nc`:function +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + +[c+2] # Ions +type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion_3d +charge = 2.0 +AA = 6.0 +thermal_conduction = true # in evolve_pressure +p_div_v = true +dissipative = false +anomalous_D_par = 1e4 +n_lowsource = 1e15 + +[Nc+2] +function = (1e16 / `hermes`:Nnorm) +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[NVc+2] +function = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[Pc+2] +function = (10.0 / `hermes`:Tnorm) * `Nc+2`:function +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + +[c+3] # Ions +type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion_3d +charge = 3.0 +AA = 6.0 +thermal_conduction = true # in evolve_pressure +p_div_v = true +dissipative = false +anomalous_D_par = 1e4 +n_lowsource = 1e15 + +[Nc+3] +function = (1e16 / `hermes`:Nnorm) +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[NVc+3] +function = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[Pc+3] +function = (10.0 / `hermes`:Tnorm) * `Nc+3`:function +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + +[c+4] # Ions +type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion_3d +charge = 4.0 +AA = 6.0 +thermal_conduction = true # in evolve_pressure +p_div_v = true +dissipative = false +anomalous_D_par = 1e4 +n_lowsource = 1e15 + +[Nc+4] +function = (1e16 / `hermes`:Nnorm) +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[NVc+4] +function = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[Pc+4] +function = (10.0 / `hermes`:Tnorm) * `Nc+4`:function +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + + + + +[c+5] # Ions +type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion_3d +charge = 5.0 +AA = 6.0 +thermal_conduction = true # in evolve_pressure +p_div_v = true +dissipative = false +anomalous_D_par = 1e4 +n_lowsource = 1e15 + +[Nc+5] +function = (1e16 / `hermes`:Nnorm) +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[NVc+5] +function = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[Pc+5] +function = (10.0 / `hermes`:Tnorm) * `Nc+5`:function +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 + + + + +[c+6] # Ions +type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion_3d +charge = 6.0 +AA = 6.0 +thermal_conduction = true # in evolve_pressure +p_div_v = true +dissipative = false +anomalous_D_par = 1e4 +n_lowsource = 1e15 + +[Nc+6] +function = (1e16 / `hermes`:Nnorm) +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[NVc+6] +function = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 +[Pc+6] +function = (10.0 / `hermes`:Tnorm) * `Nc+6`:function +source = 0.0 +bndry_all = neumann +bndry_par_all = parallel_neumann_o1 diff --git a/tests/integrated/1D-Cimpurities-fci/runtest b/tests/integrated/1D-Cimpurities-fci/runtest new file mode 100755 index 000000000..c0942bc79 --- /dev/null +++ b/tests/integrated/1D-Cimpurities-fci/runtest @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 + +# Python script to run and analyse MMS test + +from __future__ import division +from __future__ import print_function + +try: + from builtins import str +except: + pass + +from boututils.run_wrapper import shell, launch_safe, getmpirun +from boutdata.collect import collect + +from numpy import sqrt, max, abs, mean, array, log, concatenate +import numpy as np +import tarfile + +gridpath = "../../../../grids/straight_slab_120m/" +archive_path = "Straight_slab_120m.xz" + +try: + with tarfile.open(gridpath + archive_path, "r:xz") as tar: + tar.extractall(path=gridpath) +except: + raise FileNotFoundError("Could not fine mms slab y file to extract") + + +def gridname(ny): + return gridpath + f"MMS_straight_slab_8_{ny}_4_True.fci.grid.nc" + + +# Copy the json files into the data folder +import shutil + +src = "../../../json_database" +dst = "json_database" + +shutil.copytree(src, dst, dirs_exist_ok=True) + +ny = 128 +nout = 5 +timestep = 1e3 / nout +nproc = 1 + +args = ( + "mesh:file=" + + gridname(ny) + + " nout=" + + str(nout) + + " timestep=" + + str(timestep) +) + +print("Running with " + args) + +# Delete old data +shell("rm data/BOUT.dmp.*.nc") + +# Command to run +cmd = "../../../hermes-3 " + args +# Launch using MPI +s, out = launch_safe(cmd, nproc=nproc, pipe=True) + +# Save output to log file +f = open("run.log." + str(ny), "w") +f.write(out) +f.close() + +success = True + +Pi = np.mean(collect( "Ph+", tind=[nout, nout], path="data", info=False)[-1],axis=(0,2)) +Ni = np.mean(collect( "Nh+", tind=[nout, nout], path="data", info=False)[-1],axis=(0,2)) + +if success: + print(" => Test passed") + exit(0) +else: + print(" => Test failed") + exit(1) diff --git a/tests/integrated/1D-neutral-fci/data/BOUT.inp b/tests/integrated/1D-neutral-fci/data/BOUT.inp index b9fdd1daa..2db7fffe9 100644 --- a/tests/integrated/1D-neutral-fci/data/BOUT.inp +++ b/tests/integrated/1D-neutral-fci/data/BOUT.inp @@ -12,8 +12,8 @@ type = fci [solver] mxstep = 10000 -rtol = 1e-5 -atol = 1e-9 +rtol = 1e-6 +atol = 1e-10 mms = true # Run with MMS sources and output diagnostics mms_initialise = true diff --git a/tests/unit/test_adas_reactions.hxx b/tests/unit/test_adas_reactions.hxx new file mode 100644 index 000000000..d3848d5a7 --- /dev/null +++ b/tests/unit/test_adas_reactions.hxx @@ -0,0 +1,244 @@ +#ifndef TEST_ADAS_REACTIONS_H__ +#define TEST_ADAS_REACTIONS_H__ + +#include "test_reactions.hxx" + +#include "../../include/adas_carbon.hxx" +#include "../../include/adas_lithium.hxx" +#include "../../include/adas_neon.hxx" +#include + +/** + * @brief Base fixture for tests of OpenADAS subclasses. + * + * @tparam RTYPE The reaction class; must derive from OpenADAS. + */ +template +class ADASIznRecReactionTest : public IznRecReactionTest { + + static_assert(std::is_base_of(), + "Template arg to ADASReactionTest must derive from OpenADAS"); + +public: + ADASIznRecReactionTest(std::string lbl, std::string reaction_str) + : IznRecReactionTest(lbl, reaction_str) {} +}; + +// C izn +class ADASCIznTest : public ADASIznRecReactionTest> { +public: + ADASCIznTest() + : ADASIznRecReactionTest>("CIzn", "c + e -> c+ + 2e") {} +}; + +class ADASCpIznTest : public ADASIznRecReactionTest> { +public: + ADASCpIznTest() : ADASIznRecReactionTest("CpIzn", "c+ + e -> c+2 + 2e") {} +}; + +class ADASC5pIznTest : public ADASIznRecReactionTest> { +public: + ADASC5pIznTest() : ADASIznRecReactionTest("C5pIzn", "c+5 + e -> c+6 + 2e") {} +}; + +// C rec +class ADASCRecTest : public ADASIznRecReactionTest> { +public: + ADASCRecTest() + : ADASIznRecReactionTest>("CRec", "c+ + e -> c") {} +}; + +class ADASCpRecTest : public ADASIznRecReactionTest> { +public: + ADASCpRecTest() + : ADASIznRecReactionTest>("CpRec", "c+2 + e -> c+") {} +}; + +class ADASC5pRecTest : public ADASIznRecReactionTest> { +public: + ADASC5pRecTest() + : ADASIznRecReactionTest>("C5pRec", "c+6 + e -> c+5") {} +}; + +// Li izn +class ADASLiIznTest : public ADASIznRecReactionTest> { +public: + ADASLiIznTest() + : ADASIznRecReactionTest>("LiIzn", "li + e -> li+ + 2e") {} +}; + +class ADASLipIznTest : public ADASIznRecReactionTest> { +public: + ADASLipIznTest() + : ADASIznRecReactionTest>("LipIzn", + "li+ + e -> li+2 + 2e") {} +}; + +class ADASLi2pIznTest : public ADASIznRecReactionTest> { +public: + ADASLi2pIznTest() + : ADASIznRecReactionTest>("Li2pIzn", + "li+2 + e -> li+3 + 2e") {} +}; + +// Li rec +class ADASLiRecTest : public ADASIznRecReactionTest> { +public: + ADASLiRecTest() + : ADASIznRecReactionTest>("LiRec", "li+ + e -> li") {} +}; + +class ADASLipRecTest : public ADASIznRecReactionTest> { +public: + ADASLipRecTest() + : ADASIznRecReactionTest>("LipRec", "li+2 + e -> li+") { + } +}; + +class ADASLi2pRecTest : public ADASIznRecReactionTest> { +public: + ADASLi2pRecTest() + : ADASIznRecReactionTest>("Li2pRec", + "li+3 + e -> li+2") {} +}; + +// Ne izn +class ADASNeIznTest : public ADASIznRecReactionTest> { +public: + ADASNeIznTest() + : ADASIznRecReactionTest>("NeIzn", "ne + e -> ne+ + 2e") {} +}; + +class ADASNepIznTest : public ADASIznRecReactionTest> { +public: + ADASNepIznTest() + : ADASIznRecReactionTest>("NepIzn", "ne+ + e -> ne+2 + 2e") {} +}; + +class ADASNe9pIznTest : public ADASIznRecReactionTest> { +public: + ADASNe9pIznTest() + : ADASIznRecReactionTest>("Ne9pIzn", + "ne+9 + e -> ne+10 + 2e") {} +}; + +// Ne rec +class ADASNeRecTest : public ADASIznRecReactionTest> { +public: + ADASNeRecTest() + : ADASIznRecReactionTest>("NeRec", "ne+ + e -> ne") {} +}; + +class ADASNepRecTest : public ADASIznRecReactionTest> { +public: + ADASNepRecTest() + : ADASIznRecReactionTest>("NepRec", "ne+2 + e -> ne+") {} +}; + +class ADASNe9pRecTest : public ADASIznRecReactionTest> { +public: + ADASNe9pRecTest() + : ADASIznRecReactionTest>("Ne9pRec", "ne+10 + e -> ne+9") { + } +}; + +/** + * @brief Base fixture class to test ADAS CX reactions + * + * @tparam RTYPE the reaction class + */ +template +class ADASCXReactionTest : public CXReactionTest { +protected: + ADASCXReactionTest(const std::string& lbl, const std::string& reaction_str, + const std::string& neutral_sp_in, const std::string& ion_sp_in, + const std::string& neutral_sp_out, const std::string& ion_sp_out) + : CXReactionTest(lbl, reaction_str, neutral_sp_in, ion_sp_in, neutral_sp_out, + ion_sp_out) {} + // Add n_e, T_e to the input state + virtual Options generate_state() override { + Options state = CXReactionTest::generate_state(); + state["species"]["e"]["density"] = FieldFactory::get()->create3D( + this->gen_lin_field_str(this->logn_min, this->logn_max, linfunc_axis::z), &state, + mesh); + state["species"]["e"]["temperature"] = FieldFactory::get()->create3D( + this->gen_lin_field_str(this->logT_min, this->logT_max, linfunc_axis::x), &state, + mesh); + return state; + } +}; + +// C CX +template +class ADASCCXReactionTest : public ADASCXReactionTest> { +public: + ADASCCXReactionTest(std::string lbl, std::string reaction_str) + : ADASCXReactionTest>( + lbl, reaction_str, {carbon_species_name}, {Hisotope}, + {carbon_species_name}, {Hisotope, '+'}) {} +}; + +class ADASCpHCXTest : public ADASCCXReactionTest<0, 'h'> { +public: + ADASCpHCXTest() : ADASCCXReactionTest<0, 'h'>("CpHCXCX", "c+ + h -> c + h+") {} +}; + +class ADASCp3DCXTest : public ADASCCXReactionTest<2, 'd'> { +public: + ADASCp3DCXTest() : ADASCCXReactionTest<2, 'd'>("Cp3DCX", "c+3 + d -> c+2 + d+") {} +}; +class ADASCp5TCXTest : public ADASCCXReactionTest<4, 't'> { +public: + ADASCp5TCXTest() : ADASCCXReactionTest<4, 't'>("Cp5TCX", "c+5 + t -> c+4 + t+") {} +}; + +// Li CX +template +class ADASLiCXReactionTest : public ADASCXReactionTest> { +public: + ADASLiCXReactionTest(std::string lbl, std::string reaction_str) + : ADASCXReactionTest>( + lbl, reaction_str, {lithium_species_name}, {Hisotope}, + {lithium_species_name}, {Hisotope, '+'}) {} +}; + +class ADASLipHCXTest : public ADASLiCXReactionTest<0, 'h'> { +public: + ADASLipHCXTest() : ADASLiCXReactionTest<0, 'h'>("LipHCX", "li+ + h -> li + h+") {} +}; + +class ADASLip2DCXTest : public ADASLiCXReactionTest<1, 'd'> { +public: + ADASLip2DCXTest() : ADASLiCXReactionTest<1, 'd'>("Lip2DCX", "li+2 + d -> li+ + d+") {} +}; + +class ADASLip3TCXTest : public ADASLiCXReactionTest<2, 't'> { +public: + ADASLip3TCXTest() : ADASLiCXReactionTest<2, 't'>("Lip3TCX", "li+3 + t -> li+2 + t+") {} +}; + +// Ne CX +template +class ADASNeCXReactionTest : public ADASCXReactionTest> { +public: + ADASNeCXReactionTest(std::string lbl, std::string reaction_str) + : ADASCXReactionTest>( + lbl, reaction_str, {neon_species_name}, {Hisotope}, + {neon_species_name}, {Hisotope, '+'}) {} +}; + +class ADASNepHCXTest : public ADASNeCXReactionTest<0, 'h'> { +public: + ADASNepHCXTest() : ADASNeCXReactionTest<0, 'h'>("NepHCX", "ne+ + h -> ne + h+") {} +}; + +class ADASNep5DCXTest : public ADASNeCXReactionTest<4, 'd'> { +public: + ADASNep5DCXTest() : ADASNeCXReactionTest<4, 'd'>("Nep5DCX", "ne+5 + d -> ne+4 + d+") {} +}; + +class ADASNep9TCXTest : public ADASNeCXReactionTest<8, 't'> { +public: + ADASNep9TCXTest() : ADASNeCXReactionTest<8, 't'>("Nep9TCX", "ne+9 + t -> ne+8 + t+") {} +}; +#endif