diff --git a/CMakeLists.txt b/CMakeLists.txt index 9eae66871..d5f36ee5c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,6 +58,7 @@ set(HERMES_SOURCES src/div_ops.cxx src/fci_gridcheck.cxx src/neutral_mixed.cxx + src/neutral_full_velocity_curv.cxx src/full_velocity.cxx src/electromagnetic.cxx src/electron_force_balance.cxx @@ -149,6 +150,7 @@ set(HERMES_SOURCES include/ionisation.hxx include/neutral_boundary.hxx include/neutral_mixed.hxx + include/neutral_full_velocity_curv.hxx include/neutral_parallel_diffusion.hxx include/solkit_neutral_parallel_diffusion.hxx include/noflow_boundary.hxx @@ -298,6 +300,10 @@ if(HERMES_TESTS) DEPENDS setup-test) endif() endfunction() + hermes_add_integrated_test(2D-neutral-full-density-fci + REQUIRES BOUT_ENABLE_METRIC_3D) + hermes_add_integrated_test(2D-neutral-full-momentum-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/hermes-3.cxx b/hermes-3.cxx index 52df94705..b297b9b86 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -273,7 +273,7 @@ int Hermes::init(bool restarting) { // try loading J from the grid, otherwise use the one calculated from the metric coefficients - Field3D Jtmp = 0.0; + Coordinates::FieldMetric Jtmp = 0.0; if (mesh->get(Jtmp, "J_new")==0){ mesh->communicate(Jtmp); coord->J = Jtmp; @@ -291,7 +291,7 @@ int Hermes::init(bool restarting) { coord->J_perp = sqrt(coord->g_11 * coord->g_33 - coord->g_13 * coord->g_13); // try loading g_22 at the lower and upper cell interface from the grid, otherwise caluculate from the mean of the two cellcentered ones - Field3D loadtmp = 0.0; + Coordinates::FieldMetric loadtmp = 0.0; if (mesh->get(loadtmp, "g_22_cell_ylow")==0) { coord->g_22_cell_ylow = loadtmp / SQ(rho_s0); } else { 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/neutral_full_velocity_curv.hxx b/include/neutral_full_velocity_curv.hxx new file mode 100644 index 000000000..f2cdcf7dc --- /dev/null +++ b/include/neutral_full_velocity_curv.hxx @@ -0,0 +1,331 @@ + +#pragma once +#ifndef NEUTRAL_FULL_VELOCITY_CURV_H +#define NEUTRAL_FULL_VELOCITY_CURV_H + +#include +#include + +#include + +#include "component.hxx" +#include + + + +/// Evolve density, parallel momentum and pressure +/// for a neutral gas species with cross-field diffusion +struct NeutralFullVelocityCurv : public Component { + /// + /// @param name The name of the species e.g. "h" + /// @param options Top-level options. Settings will be taken from options[name] + /// @param solver Time-integration solver to be used + NeutralFullVelocityCurv(const std::string& name, Options& options, Solver *solver); + + /// Modify the given simulation state + void transform(Options &state) override; + + /// Use the final simulation state to update internal state + /// (e.g. time derivatives) + void finally(const Options &state) override; + + /// Add extra fields for output, or set attributes e.g docstrings + void outputVars(Options &state) override; + + /// Preconditioner + void precon(const Options &state, BoutReal gamma) override; +private: + std::string name; ///< Species name + YBoundary yboundary; + std::shared_ptr dagp; + Field3D Nn, Pn, NVn; // Density, pressure and parallel momentum + Field3D NVn_x, NVn_z; + Field3D Vn, Vn_x, Vn_z; ///< Neutral parallel velocity + Field3D Tn; ///< Neutral temperature + Field3D Nnlim, Pnlim, logPnlim, Vnlim, Tnlim, Vn_xlim, Vn_zlim; // Limited in regions of low density + bool isMMS; + BoutReal AA; ///< Atomic mass (proton = 1) + BoutReal n_lowsource, T_lowsource, lowsource_scale; + + BoutReal flux_limit; + + Field3D Rnn; + Field3D Dnn; + BoutReal neutral_lmax; + + BoutReal temperature_floor; + + bool dissipative; + BoutReal density_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. + BoutReal pressure_floor; ///< Minimum Pn used when dividing Pn by Nn to get Tn. + bool disable_dndt; + bool include_D, include_nu; + Field3D anomalous_D, anomalous_nu; + bool inherited_T; + bool parallel_dirichlet; + bool neutral_viscosity; ///< include viscosity? + bool neutral_conduction; ///< Include heat conduction? + bool evolve_momentum; ///< Evolve parallel momentum? + bool evolve_momentum_xz; + bool evolve_pressure; + bool momentum_advection; + bool momentum_loss; + Field3D initial_Tn; + Field3D initial_Vn, initial_Vn_x,initial_Vn_z; + + bool cross_terms; + + bool use_finite_difference; + Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity + + bool precondition {true}; ///< Enable preconditioner? + bool lax_flux; ///< Use Lax flux for advection terms + std::unique_ptr inv; ///< Laplacian inversion used for preconditioning + + Field3D density_source, pressure_source; ///< External input source + Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source + Field3D sound_speed; ///< Sound speed for use with Lax flux + + bool output_ddt; ///< Save time derivatives? + bool diagnose; ///< Save additional diagnostics? + + + // Flow diagnostics + Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; + Field3D mf_adv_perp_xlow, mf_adv_perp_ylow, mf_adv_par_ylow; + Field3D mf_visc_perp_xlow, mf_visc_perp_ylow, mf_visc_par_ylow; + Field3D ef_adv_perp_xlow, ef_adv_perp_ylow, ef_adv_par_ylow; + Field3D ef_cond_perp_xlow, ef_cond_perp_ylow, ef_cond_par_ylow; + + + + const Field3D Grad_x(Field3D& a) { + Mesh* mesh = a.getMesh(); + Coordinates* coord = mesh->getCoordinates(); + return DDX(a) / sqrt(coord->g_11); + } + + + const Field3D Grad_z(Field3D& a) { + Mesh* mesh = a.getMesh(); + Coordinates* coord = mesh->getCoordinates(); + return DDZ(a) / sqrt(coord->g_33); + } + + + struct Stencil1D { + // Cell centre values + BoutReal c, m, p, mm, pp; + // Left and right cell face values + BoutReal L, R; + }; + + BoutReal minmod(BoutReal a, BoutReal b) { + if (a * b <= 0.0) + return 0.0; + if (fabs(a) < fabs(b)) + return a; + return b; + } + + BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { + // If any of the signs are different, return zero gradient + if ((a * b <= 0.0) || (a * c <= 0.0)) { + return 0.0; + } + // Return the minimum absolute value + return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c)); + } + + void MinMod(Stencil1D& n, const BoutReal h) { + // Choose the gradient within the cell + // as the minimum (smoothest) solution + BoutReal slope = minmod(n.p - n.c, n.c - n.m); + n.L = n.c - 0.5 * slope; // 0.25*(n.p - n.m); + n.R = n.c + 0.5 * slope; // 0.25*(n.p - n.m) + } + + // Monotonized Central limiter (Van-Leer) + void MC(Stencil1D& n, const BoutReal h) { + BoutReal slope = minmod(2. * (n.p - n.c), 0.5 * (n.p - n.m), 2. * (n.c - n.m)); + n.L = n.c - 0.5 * slope; + n.R = n.c + 0.5 * slope; + } + + + + const Field3D Div_perp(Field3D& f, Field3D& vx, Field3D& vz, Field3D& spd, const bool& dissipation=false) { + + Mesh* mesh = f.getMesh(); + Coordinates* coord = mesh->getCoordinates(); + + Field3D cellcross_x = coord->cellvolume / (coord->dx * sqrt(coord->g_11)); + Field3D cellcross_z = coord->cellvolume / (coord->dz * sqrt(coord->g_33)); + + Field3D result{zeroFrom(f)}; + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + const auto ixp = i.xp(); + const auto ixpp = i.xpp(); + const auto ixm = i.xm(); + const auto ixmm = i.xmm(); + + const auto izp = i.zp(); + const auto izpp = i.zpp(); + const auto izm = i.zm(); + const auto izmm = i.zmm(); + + BoutReal cellarea_xup = 0.5 * (cellcross_x[i] + cellcross_x[ixp]); + BoutReal cellarea_xdown = 0.5 * (cellcross_x[i] + cellcross_x[ixm]); + + BoutReal cellarea_zup = 0.5 * (cellcross_z[i] + cellcross_z[izp]); + BoutReal cellarea_zdown = 0.5 * (cellcross_z[i] + cellcross_z[izm]); + + BoutReal v_xup = 0.5 * (vx[i] + vx[ixp]); + BoutReal v_xdown = 0.5 * (vx[i] + vx[ixm]); + + BoutReal v_zup = 0.5 * (vz[i] + vz[izp]); + BoutReal v_zdown = 0.5 * (vz[i] + vz[izm]); + + ////////////////////////////////////////////////////////////////////////7 + // x- direction + + Stencil1D sx; + sx.c = f[i]; + sx.m = f[ixm]; + sx.p = f[ixp]; + MinMod(sx, coord->dx[i]); + + Stencil1D sxm; + sxm.c = f[ixm]; + sxm.m = f[ixmm]; + sxm.p = f[i]; + MinMod(sxm, coord->dx[i]); + + Stencil1D sxp; + sxp.c = f[ixp]; + sxp.m = f[i]; + sxp.p = f[ixpp]; + MinMod(sxp, coord->dx[i]); + + BoutReal amax_xup = BOUTMAX(fabs(vx[i]),fabs(vx[ixp]),spd[i],spd[ixp]); + BoutReal flux_xup = (0.5 * (sx.R + sxp.L) * v_xup + 0.5 * amax_xup * (sx.R - sxp.L)) * cellarea_xup; + + BoutReal amax_xdown = BOUTMAX(fabs(vx[i]),fabs(vx[ixm]),spd[i],spd[ixm]); + BoutReal flux_xdown = (0.5 * (sx.L + sxm.R) * v_xdown + 0.5 * amax_xdown * (sxm.R - sx.L)) * cellarea_xdown; + + + + ///////////////////////////////////////////////////////// + // Z direction + + Stencil1D sz; + sz.c = f[i]; + sz.m = f[izm]; + sz.p = f[izp]; + MinMod(sz, coord->dz[i]); + + Stencil1D szm; + szm.c = f[izm]; + szm.m = f[izmm]; + szm.p = f[i]; + MinMod(szm, coord->dz[i]); + + Stencil1D szp; + szp.c = f[izp]; + szp.m = f[i]; + szp.p = f[izpp]; + MinMod(szp, coord->dz[i]); + + + BoutReal amax_zup = BOUTMAX(fabs(vz[i]),fabs(vz[izp]),spd[i],spd[izp]); + BoutReal flux_zup = (0.5 * (sz.R + szp.L) * v_zup + 0.5 * amax_zup * (sz.R - szp.L)) * cellarea_zup; + + BoutReal amax_zdown = BOUTMAX(fabs(vz[i]),fabs(vz[izm]),spd[i],spd[izm]); + BoutReal flux_zdown = (0.5 * (sz.L + szm.R) * v_zdown + 0.5 * amax_zdown * (szm.R - sz.L)) * cellarea_zdown; + + result[i] = (flux_xup + flux_zup - flux_xdown - flux_zdown) / coord->cellvolume[i]; + + } + + return result; + } + + + + const Field3D Div_perp_flows(Field3D& f, Field3D& v, Field3D& vx, Field3D& vz, Field3D& spd) { + Mesh* mesh = f.getMesh(); + Coordinates* coord = mesh->getCoordinates(); + Field3D cellcross_x = coord->cellvolume / (coord->dx * sqrt(coord->g_11)); + Field3D cellcross_z = coord->cellvolume / (coord->dz * sqrt(coord->g_33)); + Field3D result{zeroFrom(f)}; + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + const auto ixp = i.xp(); + const auto ixm = i.xm(); + + const auto izp = i.zp(); + const auto izm = i.zm(); + + BoutReal cellarea_xup = 0.5 * (cellcross_x[i] + cellcross_x[ixp]); + BoutReal cellarea_xdown = 0.5 * (cellcross_x[i] + cellcross_x[ixm]); + + BoutReal cellarea_zup = 0.5 * (cellcross_z[i] + cellcross_z[izp]); + BoutReal cellarea_zdown = 0.5 * (cellcross_z[i] + cellcross_z[izm]); + + BoutReal f_xup = 0.5 * (f[i] + f[ixp]); + BoutReal f_xdown = 0.5 * (f[i] + f[ixm]); + + BoutReal f_zup = 0.5 * (f[i] + f[izp]); + BoutReal f_zdown = 0.5 * (f[i] + f[izm]); + + BoutReal vx_up = 0.5 * (vx[i] + vx[ixp]); + BoutReal vx_down = 0.5 * (vx[i] + vx[ixm]); + + BoutReal vz_up = 0.5 * (vz[i] + vz[izp]); + BoutReal vz_down = 0.5 * (vz[i] + vz[izm]); + + BoutReal v_xup = (0.5 * (v[i] + v[ixp])); + BoutReal v_xdown = (0.5 * (v[i] + v[ixm])); + BoutReal v_zup = (0.5 * (v[i] + v[izp])); + BoutReal v_zdown = (0.5 * (v[i] + v[izm])); + + BoutReal amax_xup = BOUTMAX(fabs(vx[i]),fabs(vx[ixp]),spd[i],spd[ixp]); + BoutReal amax_xdown = BOUTMAX(fabs(vx[i]),fabs(vx[ixm]),spd[i],spd[ixm]); + + BoutReal amax_zup = BOUTMAX(fabs(vz[i]),fabs(vz[izp]),spd[i],spd[izp]); + BoutReal amax_zdown = BOUTMAX(fabs(vz[i]),fabs(vz[izm]),spd[i],spd[izm]); + + + BoutReal flux_xup = (f_xup * vx_up * v_xup + 0.5 * amax_xup * f_xup * (v[i] - v[ixp]) ) * cellarea_xup; + BoutReal flux_xdown = (f_xdown * vx_down * v_xdown - 0.5 * amax_xdown * f_xdown * (v[i] - v[ixm])) * cellarea_xdown; + + BoutReal flux_zup = (f_zup * vz_up * v_zup + 0.5 * amax_zup * f_zup * (v[i] - v[izp]) ) * cellarea_zup; + BoutReal flux_zdown = (f_zdown * vz_down * v_zdown - 0.5 * amax_zdown * f_zdown * (v[i] - v[izm]) ) * cellarea_zdown; + + result[i] = (flux_xup + flux_zup - flux_xdown -flux_zdown) / coord->cellvolume[i]; + } + + return result; + + + } + + + Field3D Div_a_Grad_perp(Field3D a, Field3D b, bool use_finite=false) { + if (a.isFci()) { + if (use_finite) { + return Div_a_Grad_perp_curv(a,b); + } else { + return (*dagp)(a, b, false); + } + } + return FV::Div_a_Grad_perp(a, b); + } +}; + +namespace { +RegisterComponent registersolverneutralmixed("neutral_full_velocity_curv"); +} + +#endif // NEUTRAL_FULL_VELOCITY_CURV_H diff --git a/src/neutral_full_velocity_curv.cxx b/src/neutral_full_velocity_curv.cxx new file mode 100644 index 000000000..84fe57aed --- /dev/null +++ b/src/neutral_full_velocity_curv.cxx @@ -0,0 +1,687 @@ + +#include +#include +#include +#include +#include + +#include "../include/div_ops.hxx" +#include "../include/hermes_build_config.hxx" +#include "../include/neutral_full_velocity_curv.hxx" + +using bout::globals::mesh; + +using ParLimiter = FV::Upwind; + + +NeutralFullVelocityCurv::NeutralFullVelocityCurv(const std::string& name, Options& alloptions, Solver* solver) + : name(name) { + AUTO_TRACE(); + + // Normalisations + const Options& units = alloptions["units"]; + const BoutReal meters = units["meters"]; + const BoutReal seconds = units["seconds"]; + const BoutReal Nnorm = units["inv_meters_cubed"]; + const BoutReal Tnorm = units["eV"]; + const BoutReal Omega_ci = 1. / units["seconds"].as(); + + + // Need to take derivatives in X for cross-field diffusion terms + ASSERT0(mesh->xstart > 0); + + auto& options = alloptions[name]; + yboundary.init(options); + + AA = options["AA"].doc("Particle atomic mass. Proton = 1").withDefault(1.0); + + // Evolving variables e.g name is "h" or "h+" + solver->add(Nn, std::string("N") + name); + + + evolve_momentum = options["evolve_momentum"] + .doc("Evolve parallel neutral momentum?") + .withDefault(true); + + evolve_momentum_xz = options["evolve_momentum_xz"] + .doc("Evolve poloidal neutral momentum?") + .withDefault(true); + + momentum_advection = options["momentum_advection"] + .doc("Use non-linear advection of momentum?") + .withDefault(false); + + momentum_loss = options["momentum_loss"] + .doc("Use perpendicular momentu loss by CX with plasma, assuming plasma has no perpendicular velocity?") + .withDefault(false); + + evolve_pressure = options["evolve_pressure"] + .doc("Evolve neutral pressure?") + .withDefault(false); + + inherited_T = options["inherited_T"] + .doc("Inherit temperature from h+?") + .withDefault(false); + + isMMS = options["isMMS"] + .doc("Is this MMS? If yes, stop sources and sinks") + .withDefault(false); + + cross_terms = options["cross_terms"] + .doc("Use cross terms for momentum transport?") + .withDefault(true); + + if (evolve_momentum) { + solver->add(NVn, std::string("NV") + name); + } else { + output_warn.write( + "WARNING: Not evolving neutral parallel momentum. NVn and Vn set to zero\n"); + initial_Vn = options["initial_Vn"] + .doc("Initial neutral velocity in parallel direction?") + .withDefault(Field3D{0.0}) / (meters/seconds); + Vn = initial_Vn; + NVn = AA * Nn * Vn; + + } + + + + if (evolve_momentum_xz) { + solver->add(NVn_x, std::string("NVx") + name); + solver->add(NVn_z, std::string("NVz") + name); + } else { + output_warn.write( + "WARNING: Not evolving neutral parallelpoloidal momentum. Set to zero\n"); + initial_Vn_x = options["initial_Vn_x"] + .doc("Initial neutral velocity in x?") + .withDefault(Field3D{0.0}) / (meters/seconds); + + initial_Vn_z = options["initial_Vn_z"] + .doc("Initial neutral velocity in z?") + .withDefault(Field3D{0.0}) / (meters/seconds); + + + Vn_x = initial_Vn_x; + Vn_z = initial_Vn_z; + NVn_x = AA * Nn * Vn_x; + NVn_z = AA * Nn * Vn_z; + + } + + + + + + if (evolve_pressure) { + solver->add(Pn, std::string("P") + name); + } else { + output_warn.write( + "WARNING: Not evolving neutral pressure!"); + + initial_Tn = options["initial_Tn"] + .doc("Initial neutral temperature when pressure is not evolved?") + .withDefault(Field3D{1.0}) / Tnorm; + + Tn = initial_Tn; + Pn = Tn * Nn; + } + + + + + density_floor = options["density_floor"] + .doc("A minimum density used when dividing NVn by Nn. " + "Normalised units.") + .withDefault(1e14) / Nnorm; + + dissipative = options["dissipative"] + .doc("Use strong dissipation in parallel divergence?") + .withDefault(true); + + disable_dndt = options["disable_dndt"] + .doc("Disable dn/dt?") + .withDefault(false); + + use_finite_difference = options["use_finite_difference"] + .doc("Use finite difference for perpendicular diffusion?") + .withDefault(false); + + + parallel_dirichlet = options["parallel_dirichlet"] + .doc("Use parallel dirichlet boundary conditions for the plasma?") + .withDefault(true); + + + + n_lowsource = options["n_lowsource"].withDefault(-1.0) / Nnorm; + T_lowsource = options["T_lowsource"].withDefault(-1.0) / Tnorm; + lowsource_scale = options["lowsource_scale"].withDefault(1e-5) * Omega_ci; + + neutral_lmax = options["neutral_lmax"].doc("Largest distance to the target, limits diffusion").withDefault(0.1) / meters; + + 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; + + precondition = options["precondition"] + .doc("Enable preconditioning in neutral model?") + .withDefault(false); + + lax_flux = options["lax_flux"] + .doc("Enable stabilising lax flux?") + .withDefault(true); + + + neutral_viscosity = options["neutral_viscosity"] + .doc("Include neutral gas viscosity?") + .withDefault(false); + + neutral_conduction = options["neutral_conduction"] + .doc("Include neutral gas heat conduction?") + .withDefault(false); + + flux_limit = options["flux_limit"].doc("Limit the flux by taking harmonic mean of velocity and sound speed?").withDefault(-1.0); + + BoutReal diffusion_norm = meters * meters / seconds; + + include_D = options.isSet("anomalous_D"); + anomalous_D = options["anomalous_D"].withDefault(Field3D{0.0}) / diffusion_norm; + + include_nu = options.isSet("anomalous_nu"); + anomalous_nu = options["anomalous_nu"].withDefault(Field3D{0.0}) / diffusion_norm; + + anomalous_D.applyBoundary("neumann"); + anomalous_nu.applyBoundary("neumann"); + mesh->communicate(anomalous_D, anomalous_nu); + anomalous_D.applyParallelBoundary("parallel_neumann_o1"); + anomalous_nu.applyParallelBoundary("parallel_neumann_o1"); + + if (precondition) { + inv = Laplacian::create(&options["precon_laplace"]); + inv->setCoefA(1.0); + } + + // Optionally output time derivatives + output_ddt = + options["output_ddt"].doc("Save derivatives to output?").withDefault(false); + + diagnose = + options["diagnose"].doc("Save additional diagnostics?").withDefault(false); + + + // Try to read the density source from the mesh + // Units of particles per cubic meter per second + density_source = 0.0; + mesh->get(density_source, std::string("N") + name + "_src"); + // Allow the user to override the source + density_source = + alloptions[std::string("N") + name]["source"] + .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) + .withDefault(density_source) + / (Nnorm * Omega_ci); + + // Try to read the pressure source from the mesh + // Units of Pascals per second + pressure_source = 0.0; + mesh->get(pressure_source, std::string("P") + name + "_src"); + // Allow the user to override the source + pressure_source = alloptions[std::string("P") + name]["source"] + .doc(std::string("Source term in ddt(P") + name + + std::string("). Units [N/m^2/s]")) + .withDefault(pressure_source) + / (SI::qe * Nnorm * Tnorm * Omega_ci); + + + if (Nn.isFci()) { + dagp = FCI::getDagp_fv(alloptions, mesh); + } +} + +void NeutralFullVelocityCurv::transform(Options& state) { + AUTO_TRACE(); + + Nn.applyBoundary(); + mesh->communicate(Nn); + Nn.applyParallelBoundary(); + + if (!evolve_momentum) { + NVn = AA * initial_Vn * Nn; + } + + NVn.applyBoundary(); + mesh->communicate(NVn); + NVn.applyParallelBoundary(); + + if (!evolve_momentum_xz) { + NVn_x = AA * initial_Vn_x * Nn; + NVn_z = AA * initial_Vn_z * Nn; + } + + NVn_x.applyBoundary(); + NVn_z.applyBoundary(); + mesh->communicate(NVn_x, NVn_z); + NVn_x.applyParallelBoundary(); + NVn_z.applyParallelBoundary(); + + if (!evolve_pressure) { + if (inherited_T) { + Options& allspecies = state["species"]; + Options& donor_species = allspecies["h+"]; + const auto donor_T = GET_NOBOUNDARY(Field3D, donor_species["temperature"]); + Pn = donor_T * Nn; + } else { + Pn = initial_Tn * Nn; + } + } + + Pn.applyBoundary(); + mesh->communicate(Pn); + Pn.applyParallelBoundary(); + + Nn = floor(Nn, 1e-3 * density_floor); + Pn = floor(Pn, 1e-3 * pressure_floor); + + Nnlim = floor(Nn, density_floor); + Pnlim = floor(Pn, pressure_floor); + + Vn = NVn / (AA * Nnlim); + Vn_x = NVn_x / (AA * Nnlim); + Vn_z = NVn_z / (AA * Nnlim); + + + Vnlim = 1.0 * Vn; + Vn_xlim = 1.0 * Vn_x; + Vn_zlim = 1.0 * Vn_z; + + + + Tn = Pn / Nnlim; + + ///////////////////////////////////////////////////// + // Parallel boundary conditions + if (!isMMS && parallel_dirichlet) { + TRACE("Neutral boundary conditions"); + yboundary.iter_pnts([&](auto& pnt) { + // Free boundary (constant gradient) density + pnt.dirichlet_o2(Nn, pnt.extrapolate_sheath_o2(Nn)); + + // Zero gradient temperature, heat flux added later + pnt.neumann_o2(Tn,0.0); + + // Zero-gradient pressure + pnt.neumann_o1(Pn,0.0); + pnt.neumann_o1(Pnlim,0.0); + + // No flow into wall + pnt.dirichlet_o2(Vn,0.0); + pnt.dirichlet_o2(NVn,0.0); + + }); // end yboundary.iter_pnts() + } + + + // Set values in the state + auto& localstate = state["species"][name]; + set(localstate["density"], Nn); + set(localstate["AA"], AA); // Atomic mass + set(localstate["pressure"], Pn); + set(localstate["momentum"], NVn); + set(localstate["velocity"], Vn); + set(localstate["temperature"], Tn); + + set(localstate["momentum_x"], NVn_x); + set(localstate["momentum_z"], NVn_z); + set(localstate["velocity_x"], Vn_x); + set(localstate["velocity_z"], Vn_z); +} + +void NeutralFullVelocityCurv::finally(const Options& state) { + AUTO_TRACE(); + auto& localstate = state["species"][name]; + + logPnlim = log(Pnlim); + + Field3D Tnlim = floor(Tn, temperature_floor); + + sound_speed = 0; + if (lax_flux) { + sound_speed = sqrt(Tnlim * (5. / 3) / AA); + } + + if (isMMS) { + sound_speed = 0.0; + } + + if (flux_limit > 0.0) { + BOUT_FOR(i, Nn.getRegion("RGN_NOY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + BoutReal vmax = flux_limit * sound_speed[i]; + + if (Vn.yup()[iyp] > 0.0) { + Vnlim.yup()[iyp] = Vnlim.yup()[iyp] * vmax / (Vnlim.yup()[iyp] + vmax); + } else { + Vnlim.yup()[iyp] = -fabs(Vnlim.yup()[iyp]) * vmax / (fabs(Vnlim.yup()[iyp]) + vmax); + } + + if (Vn.ydown()[iym] > 0.0) { + Vnlim.ydown()[iym] = Vnlim.ydown()[iym] * vmax / (Vnlim.ydown()[iym] + vmax); + } else { + Vnlim.ydown()[iym] = -fabs(Vnlim.ydown()[iym]) * vmax / (fabs(Vnlim.ydown()[iym]) + vmax); + } + + if (Vn_x[i] > 0.0) { + Vn_xlim[i] = Vn_x[i] * vmax / (Vn_x[i] + vmax); + } else { + Vn_xlim[i] = -fabs(Vn_x[i]) * vmax / (fabs(Vn_x[i]) + vmax); + } + + if (Vn_z[i] > 0.0) { + Vn_zlim[i] = Vn_z[i] * vmax / (Vn_z[i] + vmax); + } else { + Vn_zlim[i] = -fabs(Vn_z[i]) * vmax / (fabs(Vn_z[i]) + vmax); + } + + + + } + + } + + + Field3D Rnn = sqrt(Tnlim / AA) / neutral_lmax; + + if (localstate.isSet("collision_frequency")) { + Dnn = (Tnlim / AA) / (get(localstate["collision_frequency"]) + Rnn); + } else { + Dnn = (Tnlim / AA) / Rnn; + } + + if (isMMS) { + Dnn.applyBoundary("free_o2"); + } else { + Dnn.applyBoundary("neumann"); + } + + mesh->communicate(Dnn); + Dnn.applyParallelBoundary("parallel_neumann_o1"); + + kappa_n = (5. / 2) * Dnn*Nn; + + eta_n = AA * (2. / 5) * kappa_n; + + + ///////////////////////////////////////////////////// + // Neutral density + TRACE("Neutral density"); + Field3D dummy; + + + // Parallel advection + ddt(Nn) = -FV::Div_par_mod(Nn, Vnlim, sound_speed, dummy, dissipative, false); + + + // Perpendicular advection + ddt(Nn) += -Div_perp(Nn, Vn_xlim, Vn_zlim, sound_speed, dissipative); + + + // Sources + Sn = density_source; // Save for possible output + if (localstate.isSet("density_source")) { + Sn += get(localstate["density_source"]); + } + if (!isMMS) { + ddt(Nn) += Sn; // Always add density_source + } + + + // Lowsource + if (n_lowsource > 0.0) { + ddt(Nn) += low_sourceterm(Nn, n_lowsource, lowsource_scale); + } + + + // Anomalous diffusion + if (include_D) { + ddt(Nn) += Div_a_Grad_perp(anomalous_D, Nn, use_finite_difference); + } + + + // Disabling time derivate for testing purposes + if (disable_dndt) { + ddt(Nn) = 0.0; + } + + + ///////////////////////////////////////////////////// + // Neutral pressure + + + TRACE("Neutral pressure"); + + if (evolve_pressure) { + Field3D dummy_Pn; + // Parallel advection + if (!isMMS) { + ddt(Pn) = -(5.0 / 3.0) * FV::Div_par_mod(Pn, Vnlim, sound_speed, dummy_Pn, dissipative); // Parallel advection + } else { + ddt(Pn) = -(5.0 / 3.0) * Div_par(Pn * Vnlim); + } + ddt(Pn) += (2. / 3) * Vnlim * Grad_par(Pn); + + + // Perpendicular advection + + ddt(Pn) -= (5.0 / 3.0) * Div_perp(Pn, Vn_xlim, Vn_zlim, sound_speed, dissipative); + + ddt(Pn) += (2.0 / 3.0) * (Vn_xlim * Grad_x(Pn) + Vn_z * Grad_z(Pn)); + + + // Thermal conduction + if (neutral_conduction) { + ddt(Pn) += (2.0/3.0) * Div_par_K_Grad_par_mod(kappa_n, Tn, dummy_Pn, false); + ddt(Pn) += (2.0/3.0) * Div_a_Grad_perp(kappa_n, Tn, use_finite_difference); + } + + Sp = pressure_source; + if (localstate.isSet("energy_source")) { + Sp += (2. / 3) * get(localstate["energy_source"]); + } + + if (!isMMS) { + ddt(Pn) += Sp; + } + + if (T_lowsource > 0.0) { + ddt(Pn) += low_sourceterm(Tn, T_lowsource, lowsource_scale); + } + + + } + + + ///////////////////////////////////////////////////// + // Neutral parallel momentum + + if (evolve_momentum) { + Field3D dummy_NVn; + + + // Parallel gradient + ddt(NVn) = -Grad_par(Pn); + + + // Parallel advection + if (momentum_advection) { + ddt(NVn) += -AA * FV::Div_par_fvv(Nnlim, Vnlim, sound_speed); + } + + + // Perpendicular advection + if (include_D) { + ddt(NVn) += Div_par_K_Grad_par_mod(AA * anomalous_D * Vn, Nn, dummy_NVn, false); + } + + if (include_nu) { + ddt(NVn) += Div_par_K_Grad_par_mod(AA * anomalous_nu * Nn, Vn, dummy_NVn, false); + } + + if (neutral_viscosity) { + Field3D viscosity_source = Div_par_K_Grad_par_mod(eta_n , Vn , mf_visc_par_ylow , false); + if (cross_terms) { + viscosity_source += Div_a_Grad_perp(eta_n , Vn, use_finite_difference); + } + ddt(NVn) += viscosity_source; + if (evolve_pressure) { + ddt(Pn) += -(2. /3) * Vn * viscosity_source; + } + } + + if (localstate.isSet("momentum_source")) { + Snv = get(localstate["momentum_source"]); + ddt(NVn) += Snv; + } else { + Snv = 0; + } + + } + + + + ///////////////////////////////////////////////////// + // Neutral perpendicular momentum + + if (evolve_momentum_xz) { + Field3D dummy_NVnxz; + + + // Perpendicular gradient + ddt(NVn_x) = -Grad_x(Pn); + ddt(NVn_z) = -Grad_z(Pn); + + + // Perpendicular momentum advection + if (momentum_advection) { + ddt(NVn_x) += -AA * Div_perp_flows(Nnlim, Vn_x, Vn_xlim, Vn_zlim, sound_speed); + ddt(NVn_z) += -AA * Div_perp_flows(Nnlim, Vn_z, Vn_xlim, Vn_zlim, sound_speed); + } + + + // Anomalous diffusion + if (include_D) { + ddt(NVn_x) += Div_a_Grad_perp(AA * anomalous_D * Vn_x, Nn, use_finite_difference); + ddt(NVn_z) += Div_a_Grad_perp(AA * anomalous_D * Vn_z, Nn, use_finite_difference); + } + + + // Anomalous viscosity + if (include_nu) { + ddt(NVn_x) += Div_a_Grad_perp(AA * anomalous_nu * Nn, Vn_x, use_finite_difference); + ddt(NVn_z) += Div_a_Grad_perp(AA * anomalous_nu * Nn, Vn_z, use_finite_difference); + } + + if (momentum_loss && localstate.isSet("collision_frequency") ) { + const Options& allspecies = state["species"]; + const Options& donor_species = allspecies["h+"]; + Field3D mom_lossfactor = - get(localstate["collision_frequency"]); + ddt(NVn_x) += mom_lossfactor * NVn_x; + ddt(NVn_z) += mom_lossfactor * NVn_z; + } + + + // Collisional viscosity + if (neutral_viscosity) { + Field3D viscosity_source_x = Div_a_Grad_perp(eta_n , Vn_x, use_finite_difference); + if (cross_terms) { + viscosity_source_x += Div_par_K_Grad_par_mod(eta_n , Vn_x , dummy_NVnxz, false); + } + ddt(NVn_x) += viscosity_source_x; + if (evolve_pressure) { + ddt(Pn) += -(2. /3) * Vn_x * viscosity_source_x; + } + + Field3D viscosity_source_z = Div_a_Grad_perp(eta_n , Vn_z, use_finite_difference); + if (cross_terms) { + viscosity_source_z += Div_par_K_Grad_par_mod(eta_n , Vn_z , dummy_NVnxz, false); + } + ddt(NVn_z) += viscosity_source_z; + if (evolve_pressure) { + ddt(Pn) += -(2. /3) * Vn_z * viscosity_source_z; + } + + } + } + + + + +} + +void NeutralFullVelocityCurv::outputVars(Options& state) { + // Normalisations + auto Nnorm = get(state["Nnorm"]); + auto Tnorm = get(state["Tnorm"]); + auto Omega_ci = get(state["Omega_ci"]); + auto Cs0 = get(state["Cs0"]); + auto rho_s0 = get(state["rho_s0"]); + const BoutReal Pnorm = SI::qe * Tnorm * Nnorm; + + state[std::string("N") + name].setAttributes({{"time_dimension", "t"}, + {"units", "m^-3"}, + {"conversion", Nnorm}, + {"standard_name", "density"}, + {"long_name", name + " number density"}, + {"species", name}, + {"source", "neutral_mixed"}}); + + state[std::string("P") + name].setAttributes({{"time_dimension", "t"}, + {"units", "Pa"}, + {"conversion", Pnorm}, + {"standard_name", "pressure"}, + {"long_name", name + " pressure"}, + {"species", name}, + {"source", "neutral_mixed"}}); + + state[std::string("NV") + name].setAttributes( + {{"time_dimension", "t"}, + {"units", "kg / m^2 / s"}, + {"conversion", SI::Mp * Nnorm * Cs0}, + {"standard_name", "momentum"}, + {"long_name", name + " parallel momentum"}, + {"species", name}, + {"source", "neutral_mixed"}}); + + if (diagnose) { + set_with_attrs(state[std::string("Vx") + name], Vn_x, + {{"source", "neutral_full_velocity_curv"}}); + set_with_attrs(state[std::string("Vz") + name], Vn_z, + {{"source", "neutral_full_velocity_curv"}}); + + + + + } + + +} + +void NeutralFullVelocityCurv::precon(const Options& state, BoutReal gamma) { + if (!precondition) { + return; + } + const auto& species = state["species"][name]; + const Field3D N = get(species["density"]); + + // Set the coefficient in Div_par( B * Grad_par ) + Field3D coef = - gamma * Dnn; + + inv->setCoefD(coef); + Field3D dT = ddt(Pn); + dT.applyBoundary("neumann"); + mesh->communicate(dT); + Field3D dummy = 0.0; + ddt(Pn) = inv->solve(dT, ddt(Pn)); + +} diff --git a/tests/integrated/2D-neutral-full-density-fci/data/BOUT.inp b/tests/integrated/2D-neutral-full-density-fci/data/BOUT.inp new file mode 100644 index 000000000..a75ec91aa --- /dev/null +++ b/tests/integrated/2D-neutral-full-density-fci/data/BOUT.inp @@ -0,0 +1,47 @@ +nout = 50 +timestep = 0.1 + + +[mesh] + +extrapolate_y = false + +[mesh:paralleltransform] +type = fci + + +[solver] +mxstep = 10000 +rtol = 1e-5 +atol = 1e-9 +mms = true # Run with MMS sources and output diagnostics +mms_initialise = true +[hermes] +components = h + + +Nnorm = 1e18 +Bnorm = 1 +Tnorm = 5 + +[h] +type = neutral_full_velocity_curv + +AA = 1.0 + +use_finite_difference = true + +isMMS = true + +evolve_pressure = false +evolve_momentum = false +evolve_momentum_xz = false + +anomalous_D = 0.2*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 2.50021345462797 +initial_Vn_x = 30.0*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) +initial_Vn_z = 20.0*sin(6.28318530717959*x)*sin(4.0*z + 2.34123135) +[Nh] +solution = 0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1 + +source = 3.29714527358133e-7*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.5)*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 5.24712167486348e-6*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(6.28318530717959*x)*cos(4.0*z + 2.34123135) + 3.93534125614761e-6*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(2.0*z + 1.33221)*cos(12.5663706143592*x) + 5.24712167486348e-8*sin(6.28318530717959*x)*sin(12.5663706143592*x)*sin(4.0*z + 2.34123135)*cos(2.0*z + 1.33221) - 1.31874551249943e-8*sin(6.28318530717959*x)*sin(12.5663706143592*x)*cos(2.0*z + 1.33221)*cos(4.0*z + 4.661) + 7.87068251229522e-8*sin(12.5663706143592*x)*sin(2.0*z + 1.33221)^2*cos(12.5663706143592*x) - 3.29686378124858e-9*sin(2.0*z + 1.33221)*sin(4.0*z + 4.661)*cos(6.28318530717959*x)*cos(12.5663706143592*x) +bndry_all = dirichlet_o2(`Nh:solution`) \ No newline at end of file diff --git a/tests/integrated/2D-neutral-full-density-fci/mms.py b/tests/integrated/2D-neutral-full-density-fci/mms.py new file mode 100644 index 000000000..4992d415a --- /dev/null +++ b/tests/integrated/2D-neutral-full-density-fci/mms.py @@ -0,0 +1,65 @@ +from __future__ import division +from __future__ import print_function + +from boutdata.mms import Metric, sin, cos, Div_par, Grad_par, exprToStr, diff, y, t, DDY, x, z, DDX, DDZ, Delp2, D2DX2, sqrt +from math import pi +from sympy import log +# Length of the y domain +Ly = 2.0 * pi + +# Atomic mass number +AA = 1.0 + +# metric tensor +metric = Metric() # Identity + +qe = 1.60217663e-19 +Me = 9.1093837e-31 +e0 = 8.85418781e-12 +Pi = pi + +Tnorm = 5 +Nnorm = 1e18 +Bnorm = 1.0 +Omega_ci = qe * Bnorm / (1836.0*Me); +seconds = 1.0 / Omega_ci; +rho_s = 0.00022847 +neutral_lmax = 0.02/ rho_s +# Define solution in terms of input x,y,z + + +Nn = (1e17 + 2e16* sin(4.0 * pi * x) * sin(2 * z + 1.33221) ) / Nnorm +Vn_x = 30 * sin(4.0 * pi * x) * sin(2 * z + 1.33221) / (rho_s/seconds) +Vn_z = 20 * sin(2.0 * pi * x) * sin(4 * z + 2.34123135) / (rho_s/seconds) + +anomalous_D = 0.5 + 0.2 * sin(2.0 * pi * x) * sin(4 * z + 4.661) / ((rho_s * rho_s) / seconds) + +replace = [(x, metric.x), (z, metric.z * 2.0 * pi ) ] + +Nn = Nn.subs(replace) +Vn_x = Vn_x.subs(replace) +Vn_z = Vn_z.subs(replace) +anomalous_D = anomalous_D.subs(replace) + + + +dNndt = ( -DDX(Nn*Vn_x) -DDZ(Nn*Vn_z)) * rho_s + (DDX(anomalous_D * DDX(Nn)) + DDZ(anomalous_D * DDZ(Nn))) * (rho_s * rho_s) + +SNn = diff(Nn, t) - dNndt + +replace = [(metric.x, x), (metric.z, z / (2.0 * pi) ) ] + +Nn = Nn.subs(replace) +Vn_x = Vn_x.subs(replace) +Vn_z = Vn_z.subs(replace) +SNn = SNn.subs(replace) +anomalous_D = anomalous_D.subs(replace) + +print("anomalous_D = " + exprToStr(anomalous_D * rho_s * rho_s / seconds)) +print("initial_Vn_x = " + exprToStr(Vn_x * (rho_s/seconds))) +print("initial_Vn_z = " + exprToStr(Vn_z * (rho_s/seconds))) + +print("[Nh]") +print("solution = " + exprToStr(Nn)) +print("\nsource = " + exprToStr(SNn)) +print("bndry_all = dirichlet_o2(`Nh:solution`)") diff --git a/tests/integrated/2D-neutral-full-density-fci/runtest b/tests/integrated/2D-neutral-full-density-fci/runtest new file mode 100755 index 000000000..f5bcdcc50 --- /dev/null +++ b/tests/integrated/2D-neutral-full-density-fci/runtest @@ -0,0 +1,140 @@ +#!/usr/bin/env python3 + +# Python script to run and analyse MMS test + +from __future__ import division +from __future__ import print_function +import numpy as np + +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 tarfile + +gridpath = "../../../../grids/mms_slab_fci_xz/" +archive_path = "MMS_Slab_fci_XZ.tar.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(nx): + return gridpath + f"MMS_straight_slab_{nx}_4_{nx-4}.fci.grid.nc" + + +# List of NY values to use +# nxlist = [36, 68, 132, 260] #, 640, 1280, 2560, 5120] +nxlist = [36, 68] +nout = 10 +timestep = 1e6 / nout + +nproc = 1 + +varnames = ["Nh"] + +results = {} +for var in varnames: + results[var] = {"l2": [], "inf": []} + +for ny in nxlist: + 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() + + # Collect data + for var in varnames: + E = collect("E_" + var, tind=[nout, nout], path="data", info=False) + print(str(np.array(E).shape)) + E = E[0, 2:-2, 0, :] + + l2 = sqrt(mean(E**2)) + linf = max(abs(E)) + results[var]["l2"].append(l2) + results[var]["inf"].append(linf) + + print("Error norm %s: l-2 %f l-inf %f" % (var, l2, linf)) + +# Calculate grid spacing +dy = 1.0 / array(nxlist) + +success = True + +for var in varnames: + l2 = results[var]["l2"] + order = log(l2[-1] / l2[-2]) / log(dy[-1] / dy[-2]) + print("Convergence order %s = %f" % (var, order)) + + if order < 1.8: + success = False + +# plot errors +try: + import matplotlib.pyplot as plt + + for var in varnames: + l2 = results[var]["l2"] + inf = results[var]["inf"] + order = log(l2[-1] / l2[-2]) / log(dy[-1] / dy[-2]) + + plt.plot(dy, l2, "-o", label=r"$l_2$ (" + var + ")") + plt.plot(dy, inf, "-x", label=r"$l_\infty$ (" + var + ")") + plt.plot( + dy, l2[-1] * (dy / dy[-1]) ** order, "--", label="Order %.1f" % (order) + ) + + plt.legend(loc="upper left") + plt.grid() + + plt.yscale("log") + plt.xscale("log") + + plt.xlabel(r"Mesh spacing $\delta y$") + plt.ylabel("Error norm") + + plt.savefig("fluid_norm.pdf") + plt.savefig("fluid_norm.png") + + # plt.show() + plt.close() +except: + # Plotting could fail for any number of reasons, and the actual + # error raised may depend on, among other things, the current + # matplotlib backend, so catch everything + pass + +if success: + print(" => Test passed") + exit(0) +else: + print(" => Test failed") + exit(1) diff --git a/tests/integrated/2D-neutral-full-momentum-fci/data/BOUT.inp b/tests/integrated/2D-neutral-full-momentum-fci/data/BOUT.inp new file mode 100644 index 000000000..ab6416d58 --- /dev/null +++ b/tests/integrated/2D-neutral-full-momentum-fci/data/BOUT.inp @@ -0,0 +1,67 @@ +nout = 50 +timestep = 0.1 + + +[mesh] + +extrapolate_y = false + +[mesh:paralleltransform] +type = fci + + +[solver] +type = cvode +mxstep = 10000 +rtol = 1e-5 +atol = 1e-9 +mms = true # Run with MMS sources and output diagnostics +mms_initialise = true +[hermes] +components = h + + +Nnorm = 1e18 +Bnorm = 1 +Tnorm = 5 + +[h] +type = neutral_full_velocity_curv + +AA = 1.0 + +use_finite_difference = true + +isMMS = true +neutral_lmax = 0.02 +evolve_pressure = false +evolve_momentum = false +evolve_momentum_xz = true + +neutral_viscosity = true + +disable_dndt = true + +initial_Tn = 5.0/(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1) +anomalous_D = 0.2*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 1.5 +anomalous_nu = 0.2*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 2.5 +[Nh] +solution = 0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1 + +source = 0.0 + + +bndry_all = dirichlet_o2(`Nh:solution`) + + + +[NVxh] +solution = 0.00137070296684326*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + +source = 4.51940680861117e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.299974387631475)*sin(12.5663706143592*x)^2*sin(2.0*z + 1.33221)^2 - 2.25970340430559e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.299974387631475)*sin(12.5663706143592*x)^2*cos(2.0*z + 1.33221)^2 - 2.25970340430559e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.299974387631475)*sin(2.0*z + 1.33221)^2*cos(12.5663706143592*x)^2 + 2.25970340430559e-8*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.499957312719125)*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) - 2.25970340430559e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.499957312719125)*sin(12.5663706143592*x)^2*cos(2.0*z + 1.33221)^2 - 2.25970340430559e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.499957312719125)*sin(2.0*z + 1.33221)^2*cos(12.5663706143592*x)^2 - 9.03804193247104e-10*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(6.28318530717959*x)*sin(12.5663706143592*x)*cos(2.0*z + 1.33221)*cos(4.0*z + 4.661) - 2.25951048311776e-10*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(2.0*z + 1.33221)*sin(4.0*z + 4.661)*cos(6.28318530717959*x)*cos(12.5663706143592*x) - 9.89059134374573e-7*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(12.5663706143592*x)^2*cos(2.0*z + 1.33221)^2/((0.2*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 1)^2*sqrt(1/(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1))) - 9.89059134374573e-7*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(2.0*z + 1.33221)^2*cos(12.5663706143592*x)^2/((0.2*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 1)^2*sqrt(1/(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1))) - 1.80760838649421e-11*sin(6.28318530717959*x)*sin(12.5663706143592*x)^2*sin(2.0*z + 1.33221)*cos(2.0*z + 1.33221)*cos(4.0*z + 4.661) - 4.51902096623552e-12*sin(12.5663706143592*x)*sin(2.0*z + 1.33221)^2*sin(4.0*z + 4.661)*cos(6.28318530717959*x)*cos(12.5663706143592*x) + 1.97811826874915e-6*sin(12.5663706143592*x)*sin(2.0*z + 1.33221)/sqrt(1/(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)) +bndry_all = dirichlet_o2(`NVxh:solution`) +[NVzh] +solution = 0.000913801977895509*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(6.28318530717959*x)*sin(4.0*z + 2.34123135) + +source = 3.01293787240745e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.299974387631475)*sin(6.28318530717959*x)*sin(12.5663706143592*x)*sin(2.0*z + 1.33221)*sin(4.0*z + 2.34123135) - 3.01293787240745e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.299974387631475)*sin(6.28318530717959*x)*sin(12.5663706143592*x)*cos(2.0*z + 1.33221)*cos(4.0*z + 2.34123135) - 7.53234468101862e-11*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.299974387631475)*sin(2.0*z + 1.33221)*sin(4.0*z + 2.34123135)*cos(6.28318530717959*x)*cos(12.5663706143592*x) + 3.20124648943291e-8*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.499957312719125)*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(6.28318530717959*x)*sin(4.0*z + 2.34123135) - 3.01293787240745e-10*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.499957312719125)*sin(6.28318530717959*x)*sin(12.5663706143592*x)*cos(2.0*z + 1.33221)*cos(4.0*z + 2.34123135) - 7.53234468101862e-11*(0.03999658501753*sin(6.28318530717959*x)*sin(4.0*z + 4.661) + 0.499957312719125)*sin(2.0*z + 1.33221)*sin(4.0*z + 2.34123135)*cos(6.28318530717959*x)*cos(12.5663706143592*x) - 1.20507225766281e-9*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(6.28318530717959*x)^2*cos(4.0*z + 2.34123135)*cos(4.0*z + 4.661) - 7.53170161039253e-11*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(4.0*z + 2.34123135)*sin(4.0*z + 4.661)*cos(6.28318530717959*x)^2 - 1.31874551249943e-6*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(6.28318530717959*x)*sin(12.5663706143592*x)*cos(2.0*z + 1.33221)*cos(4.0*z + 2.34123135)/((0.2*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 1)^2*sqrt(1/(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1))) - 3.29686378124858e-7*(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)*sin(2.0*z + 1.33221)*sin(4.0*z + 2.34123135)*cos(6.28318530717959*x)*cos(12.5663706143592*x)/((0.2*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 1)^2*sqrt(1/(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1))) - 1.20507225766281e-11*sin(6.28318530717959*x)^2*sin(12.5663706143592*x)*sin(4.0*z + 2.34123135)*cos(2.0*z + 1.33221)*cos(4.0*z + 4.661) - 3.01268064415701e-12*sin(6.28318530717959*x)*sin(2.0*z + 1.33221)*sin(4.0*z + 2.34123135)*sin(4.0*z + 4.661)*cos(6.28318530717959*x)*cos(12.5663706143592*x) + 2.80233421406129e-6*sin(6.28318530717959*x)*sin(4.0*z + 2.34123135)/sqrt(1/(0.02*sin(12.5663706143592*x)*sin(2.0*z + 1.33221) + 0.1)) +bndry_all = dirichlet_o2(`NVzh:solution`) diff --git a/tests/integrated/2D-neutral-full-momentum-fci/mms.py b/tests/integrated/2D-neutral-full-momentum-fci/mms.py new file mode 100644 index 000000000..54f80949b --- /dev/null +++ b/tests/integrated/2D-neutral-full-momentum-fci/mms.py @@ -0,0 +1,125 @@ +from __future__ import division +from __future__ import print_function + +from boutdata.mms import Metric, sin, cos, Div_par, Grad_par, exprToStr, diff, y, t, DDY, x, z, DDX, DDZ, Delp2, D2DX2, sqrt +from math import pi +from sympy import log +# Length of the y domain +Ly = 2.0 * pi + +# Atomic mass number +AA = 1.0 + +# metric tensor +metric = Metric() # Identity + +qe = 1.60217663e-19 +Me = 9.1093837e-31 +e0 = 8.85418781e-12 +Pi = pi + +Tnorm = 5 +Nnorm = 1e18 +Bnorm = 1.0 +Omega_ci = qe * Bnorm / (1836.0*Me); +seconds = 1.0 / Omega_ci; +rho_s = 0.00022847 +neutral_lmax = 0.02/ rho_s +# Define solution in terms of input x,y,z + + +Nn = (1e17 + 2e16* sin(4.0 * pi * x) * sin(2 * z + 1.33221) ) / Nnorm +Vn_x = 30 * sin(4.0 * pi * x) * sin(2 * z + 1.33221) / (rho_s/seconds) +Vn_z = 20 * sin(2.0 * pi * x) * sin(4 * z + 2.34123135) / (rho_s/seconds) +NVn_x = AA * Nn * Vn_x +NVn_z = AA * Nn * Vn_z +#Tn = 1.0 / Tnorm +#Pn = Nn * Tn + +Pn = 1.0 + 0*x +Tn = Pn / Nn + +anomalous_D = (1.5 + 0.2 * sin(2.0 * pi * x) * sin(4 * z + 4.661)) / ((rho_s * rho_s) / seconds) +anomalous_nu = (2.5 + 0.2 * sin(2.0 * pi * x) * sin(4 * z + 4.661)) / ((rho_s * rho_s) / seconds) +#anomalous_D = 0.0 + 0*x + +replace = [(x, metric.x), (z, metric.z * 2.0 * pi ) ] + +Nn = Nn.subs(replace) +Vn_x = Vn_x.subs(replace) +Vn_z = Vn_z.subs(replace) +anomalous_D = anomalous_D.subs(replace) +anomalous_nu = anomalous_nu.subs(replace) +NVn_x = NVn_x.subs(replace) +NVn_z = NVn_z.subs(replace) +Tn = Tn.subs(replace) +Pn = Pn.subs(replace) + +Rnn = sqrt(Tn / AA) / neutral_lmax +Dnn = (Tn / AA) / Rnn +kappa_n = (5.0 / 2.0) * Dnn * Nn +eta_n = AA * (2.0 / 5.0) * kappa_n + + + +dNndt = ( -DDX(Nn*Vn_x) -DDZ(Nn*Vn_z)) * rho_s + (DDX(anomalous_D * DDX(Nn)) + DDZ(anomalous_D * DDZ(Nn))) * (rho_s * rho_s) + +x_grad = -DDX(Pn) * rho_s +x_an_D = AA * ( DDX(anomalous_D * Vn_x * DDX(Nn)) + DDZ(anomalous_D * Vn_x * DDZ(Nn)) ) * (rho_s * rho_s) +x_an_nu = AA * ( DDX(anomalous_nu * Nn * DDX(Vn_x)) +DDZ(anomalous_nu * Nn * DDZ(Vn_x)) ) * (rho_s * rho_s) +x_visc = ( DDX(eta_n * DDX(Vn_x)) +DDZ(eta_n * DDZ(Vn_x)) ) * (rho_s * rho_s) + +dNVn_xdt = x_grad + x_an_D + x_an_nu + x_visc + +z_grad = -DDZ(Pn) * rho_s +z_an_D = AA * ( DDX(anomalous_D * Vn_z * DDX(Nn)) + DDZ(anomalous_D * Vn_z * DDZ(Nn)) ) * (rho_s * rho_s) +z_an_nu = AA * ( DDX(anomalous_nu * Nn * DDX(Vn_z)) + DDZ(anomalous_nu * Nn * DDZ(Vn_z)) ) * (rho_s * rho_s) +z_visc = ( DDX(eta_n * DDX(Vn_z)) +DDZ(eta_n * DDZ(Vn_z)) ) * (rho_s * rho_s) + + +dNVn_zdt = z_grad + z_an_D + z_an_nu + z_visc + + +SNn = diff(Nn, t) - dNndt +SNVn_x = diff(NVn_x, t) - dNVn_xdt +SNVn_z = diff(NVn_z, t) - dNVn_zdt + + +replace = [(metric.x, x), (metric.z, z / (2.0 * pi) ) ] + +Nn = Nn.subs(replace) +Vn_x = Vn_x.subs(replace) +Vn_z = Vn_z.subs(replace) +NVn_x = NVn_x.subs(replace) +NVn_z = NVn_z.subs(replace) +Tn = Tn.subs(replace) + + +SNn = SNn.subs(replace) +SNVn_x = SNVn_x.subs(replace) +SNVn_z = SNVn_z.subs(replace) + +anomalous_D = anomalous_D.subs(replace) +anomalous_nu = anomalous_nu.subs(replace) + + +print("initial_Tn = " + exprToStr(Tn * Tnorm)) +print("anomalous_D = " + exprToStr(anomalous_D * rho_s * rho_s / seconds)) +print("anomalous_nu = " + exprToStr(anomalous_nu * rho_s * rho_s / seconds)) + + + +print("[Nh]") +print("solution = " + exprToStr(Nn)) +print("\nsource = " + exprToStr(SNn)) +print("bndry_all = dirichlet_o2(`Nh:solution`)") + +print("[NVxh]") +print("solution = " + exprToStr(NVn_x)) +print("\nsource = " + exprToStr(SNVn_x)) +print("bndry_all = dirichlet_o2(`NVxh:solution`)") + +print("[NVzh]") +print("solution = " + exprToStr(NVn_z)) +print("\nsource = " + exprToStr(SNVn_z)) +print("bndry_all = dirichlet_o2(`NVzh:solution`)") diff --git a/tests/integrated/2D-neutral-full-momentum-fci/runtest b/tests/integrated/2D-neutral-full-momentum-fci/runtest new file mode 100755 index 000000000..2058206ea --- /dev/null +++ b/tests/integrated/2D-neutral-full-momentum-fci/runtest @@ -0,0 +1,140 @@ +#!/usr/bin/env python3 + +# Python script to run and analyse MMS test + +from __future__ import division +from __future__ import print_function +import numpy as np + +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 tarfile + +gridpath = "../../../../grids/mms_slab_fci_xz/" +archive_path = "MMS_Slab_fci_XZ.tar.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(nx): + return gridpath + f"MMS_straight_slab_{nx}_4_{nx-4}.fci.grid.nc" + + +# List of NY values to use +# nxlist = [36, 68, 132, 260] #, 640, 1280, 2560, 5120] +nxlist = [36, 68] +nout = 10 +timestep = 1e3 / nout + +nproc = 1 + +varnames = ["NVxh", "NVzh"] + +results = {} +for var in varnames: + results[var] = {"l2": [], "inf": []} + +for ny in nxlist: + 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() + + # Collect data + for var in varnames: + E = collect("E_" + var, tind=[nout, nout], path="data", info=False) + print(str(np.array(E).shape)) + E = E[0, 2:-2, 0, :] + + l2 = sqrt(mean(E**2)) + linf = max(abs(E)) + results[var]["l2"].append(l2) + results[var]["inf"].append(linf) + + print("Error norm %s: l-2 %f l-inf %f" % (var, l2, linf)) + +# Calculate grid spacing +dy = 1.0 / array(nxlist) + +success = True + +for var in varnames: + l2 = results[var]["l2"] + order = log(l2[-1] / l2[-2]) / log(dy[-1] / dy[-2]) + print("Convergence order %s = %f" % (var, order)) + + if order < 1.8: + success = False + +# plot errors +try: + import matplotlib.pyplot as plt + + for var in varnames: + l2 = results[var]["l2"] + inf = results[var]["inf"] + order = log(l2[-1] / l2[-2]) / log(dy[-1] / dy[-2]) + + plt.plot(dy, l2, "-o", label=r"$l_2$ (" + var + ")") + plt.plot(dy, inf, "-x", label=r"$l_\infty$ (" + var + ")") + plt.plot( + dy, l2[-1] * (dy / dy[-1]) ** order, "--", label="Order %.1f" % (order) + ) + + plt.legend(loc="upper left") + plt.grid() + + plt.yscale("log") + plt.xscale("log") + + plt.xlabel(r"Mesh spacing $\delta y$") + plt.ylabel("Error norm") + + plt.savefig("fluid_norm.pdf") + plt.savefig("fluid_norm.png") + + # plt.show() + plt.close() +except: + # Plotting could fail for any number of reasons, and the actual + # error raised may depend on, among other things, the current + # matplotlib backend, so catch everything + pass + +if success: + print(" => Test passed") + exit(0) +else: + print(" => Test failed") + exit(1)