diff --git a/docs/sphinx/code_structure.rst b/docs/sphinx/code_structure.rst index a0d49347..09dd8a21 100644 --- a/docs/sphinx/code_structure.rst +++ b/docs/sphinx/code_structure.rst @@ -54,7 +54,7 @@ for new variables which are added to the state: * `charge` Charge, in units of proton charge (i.e. electron=-1) * `density` - * `momentum` + * `momentum` Parallel momentum * `pressure` * `velocity` Parallel velocity * `temperature` @@ -74,12 +74,14 @@ for new variables which are added to the state: * `fields` * `vorticity` - * `phi` Electrostatic potential - * `DivJdia` Divergence of diamagnetic current - * `DivJcol` Divergence of collisional current - * `DivJextra` Divergence of current, including 2D parallel current - closures. Not including diamagnetic, parallel current due to - flows, or polarisation currents + * `phi` Electrostatic potential + * `Apar` Electromagnetic potential b dot A in induction terms + * `Apar_flutter` The electromagnetic potential (b dot A) in flutter terms + * `DivJdia` Divergence of diamagnetic current + * `DivJcol` Divergence of collisional current + * `DivJextra` Divergence of current, including 2D parallel current + closures. Not including diamagnetic, parallel current due to + flows, or polarisation currents For example to get the electron density:: diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index 30fe024d..d9046a50 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -2044,11 +2044,19 @@ electromagnetic ~~~~~~~~~~~~~~~ **Notes**: When using this module, + 1. Set ``sound_speed:alfven_wave=true`` so that the shear Alfven wave speed is included in the calculation of the fastest parallel wave speed for numerical dissipation. -2. For tokamak simulations use zero-Laplacian boundary conditions - by setting ``electromagnetic:apar_boundary_neumann=false``. +2. For tokamak simulations use Neumann boundary condition on the core + and Dirichlet on SOL and PF boundaries by setting + ``electromagnetic:apar_core_neumann=true`` (this is the default). +3. Set the potential core boundary to be constant in Y by setting + ``vorticity:phi_core_averagey = true`` +4. Magnetic flutter terms must be enabled to be active + (``electromagnetic:magnetic_flutter=true``). They use an + ``Apar_flutter`` field, not the ``Apar`` field that is calculated + from the induction terms. This component modifies the definition of momentum of all species, to include the contribution from the electromagnetic potential @@ -2062,8 +2070,17 @@ Assumes that "momentum" :math:`p_s` calculated for all species p_s = m_s n_s v_{||s} + Z_s e n_s A_{||} which arises once the electromagnetic contribution to the force on -each species is included in the momentum equation. This is normalised -so that in dimensionless quantities +each species is included in the momentum equation. This requires +an additional term in the momentum equation: + +.. math:: + + \frac{\partial p_s}{\partial t} = \cdots + Z_s e A_{||} \frac{\partial n_s}{\partial t} + +This is implemented so that the density time-derivative is calculated using the lowest order +terms (parallel flow, ExB drift and a low density numerical diffusion term). + +The above equations are normalised so that in dimensionless quantities: .. math:: @@ -2098,5 +2115,22 @@ To convert the species momenta into a current, we take the sum of - \frac{1}{\beta_{em}} \nabla_\perp^2 A_{||} + \sum_s \frac{Z^2 n_s}{A}A_{||} = \sum_s \frac{Z}{A} p_s +The toroidal variation of density :math:`n_s` must be kept in this +equation. By default the iterative "Naulin" solver is used to do +this: A fast FFT-based method is used in a fixed point iteration, +correcting for the density variation. + +Magnetic flutter terms are disabled by default, and can be enabled by setting + +.. code-block:: ini + + [electromagnetic] + magnetic_flutter = true + +This writes an ``Apar_flutter`` field to the state, which then enables perturbed +parallel derivative terms in the ``evolve_density``, ``evolve_pressure``, ``evolve_energy`` and +``evolve_momentum`` components. Parallel flow terms are modified, and parallel heat +conduction. + .. doxygenstruct:: Electromagnetic :members: diff --git a/include/div_ops.hxx b/include/div_ops.hxx index f5a0204f..5c97adf8 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -23,8 +23,8 @@ */ -#ifndef HERMES_DIV_OPS_H -#define HERMES_DIV_OPS_H +#ifndef DIV_OPS_H +#define DIV_OPS_H #include #include @@ -44,6 +44,11 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flux = true, bool poloidal = false, bool positive = false); +/// This version has an extra coefficient 'g' that is linearly interpolated +/// onto cell faces +const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D &n, const Field3D &g, const Field3D &f, + bool bndry_flux = true, bool positive = false); + const Field3D Div_Perp_Lap_FV_Index(const Field3D& a, const Field3D& f, bool xflux); const Field3D Div_Z_FV_Index(const Field3D& a, const Field3D& f); @@ -216,9 +221,10 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, flux = bndryval * vpar * vpar; } else { // Add flux due to difference in boundary values - flux = s.R * vpar * sv.R - + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k))) - * (s.R * sv.R - bndryval * vpar); + flux = + s.R * sv.R * sv.R + + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k))) + * (s.R * sv.R - bndryval * vpar); } } else { // Maximum wave speed in the two cells @@ -244,9 +250,10 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, flux = bndryval * vpar * vpar; } else { // Add flux due to difference in boundary values - flux = s.L * vpar * sv.L - - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k))) - * (s.L * sv.L - bndryval * vpar); + flux = + s.L * sv.L * sv.L + - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k))) + * (s.L * sv.L - bndryval * vpar); } } else { // Maximum wave speed in the two cells @@ -690,4 +697,4 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Fi } // namespace FV -#endif // HERMES_DIV_OPS_H +#endif // DIV_OPS_H diff --git a/include/electromagnetic.hxx b/include/electromagnetic.hxx index eb43b70e..15a3588e 100644 --- a/include/electromagnetic.hxx +++ b/include/electromagnetic.hxx @@ -64,6 +64,13 @@ private: std::unique_ptr aparSolver; // Laplacian solver in X-Z + 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; + bool diagnose; ///< Output additional diagnostics? }; diff --git a/include/vorticity.hxx b/include/vorticity.hxx index 18fd8b94..acf2c281 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -37,6 +37,8 @@ struct Vorticity : public Component { /// Relax radial phi boundaries towards zero-gradient? /// - phi_boundary_timescale: float, 1e-4 /// Timescale for phi boundary relaxation [seconds] + /// - phi_core_averagey: bool, default false + /// Average phi core boundary in Y? (if phi_boundary_relax) /// - phi_dissipation: bool, default true /// Parallel dissipation of potential (Recommended) /// - poloidal_flows: bool, default true @@ -128,7 +130,8 @@ private: bool phi_boundary_relax; ///< Relax boundary to zero-gradient 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 split_n0; // Split phi into n=0 and n!=0 components LaplaceXY* laplacexy; // Laplacian solver in X-Y (n=0) diff --git a/src/div_ops.cxx b/src/div_ops.cxx index 1ff0eab3..08a422b4 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -1372,3 +1372,160 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, return result; } +const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D &n, const Field3D &g, const Field3D &f, + bool bndry_flux, bool positive) { + Field3D result{0.0}; + + Coordinates *coord = mesh->getCoordinates(); + + ////////////////////////////////////////// + // X-Z advection. + // + // Z + // | + // + // fmp --- vU --- fpp + // | nU | + // | | + // vL nL nR vR -> X + // | | + // | nD | + // fmm --- vD --- fpm + // + + int nz = mesh->LocalNz; + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < nz; k++) { + int kp = (k + 1) % nz; + int kpp = (kp + 1) % nz; + int km = (k - 1 + nz) % nz; + int kmm = (km - 1 + nz) % nz; + + // 1) Interpolate stream function f onto corners fmp, fpp, fpm + + BoutReal fmm = 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + + f(i - 1, j, km)); + BoutReal fmp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) + + f(i - 1, j, kp)); // 2nd order accurate + BoutReal fpp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + + f(i + 1, j, kp)); + BoutReal fpm = 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + + f(i + 1, j, km)); + + // 2) Calculate velocities on cell faces + + BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx + BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx + + BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) / + coord->dz(i, j); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) / + coord->dz(i, j); // J*df/dz + + // 3) Calculate g on cell faces + + BoutReal gU = 0.5 * (g(i, j, kp) + g(i, j, k)); + BoutReal gD = 0.5 * (g(i, j, km) + g(i, j, k)); + BoutReal gR = 0.5 * (g(i + 1, j, k) + g(i, j, k)); + BoutReal gL = 0.5 * (g(i - 1, j, k) + g(i, j, k)); + + // 4) Calculate n on the cell faces. The sign of the + // velocity determines which side is used. + + // X direction + Stencil1D s; + s.c = n(i, j, k); + s.m = n(i - 1, j, k); + s.mm = n(i - 2, j, k); + s.p = n(i + 1, j, k); + s.pp = n(i + 2, j, k); + + MC(s, coord->dx(i, j)); + + // Right side + if ((i == mesh->xend) && (mesh->lastX())) { + // At right boundary in X + + if (bndry_flux) { + BoutReal flux; + if (vR > 0.0) { + // Flux to boundary + flux = vR * s.R * gR; + } else { + // Flux in from boundary + flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)) * gR; + } + + result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + } + } else { + // Not at a boundary + if (vR > 0.0) { + // Flux out into next cell + BoutReal flux = vR * s.R * gR; + result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + } + } + + // Left side + + if ((i == mesh->xstart) && (mesh->firstX())) { + // At left boundary in X + + if (bndry_flux) { + BoutReal flux; + + if (vL < 0.0) { + // Flux to boundary + flux = vL * s.L * gL; + + } else { + // Flux in from boundary + flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)) * gL; + } + result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + } + } else { + // Not at a boundary + + if (vL < 0.0) { + BoutReal flux = vL * s.L * gL; + result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + } + } + + /// NOTE: Need to communicate fluxes + + // Z direction + s.m = n(i, j, km); + s.mm = n(i, j, kmm); + s.p = n(i, j, kp); + s.pp = n(i, j, kpp); + + MC(s, coord->dz(i, j)); + + if (vU > 0.0) { + BoutReal flux = vU * s.R * gU/ (coord->J(i, j) * coord->dz(i, j)); + result(i, j, k) += flux; + result(i, j, kp) -= flux; + } + if (vD < 0.0) { + BoutReal flux = vD * s.L * gD / (coord->J(i, j) * coord->dz(i, j)); + result(i, j, k) -= flux; + result(i, j, km) += flux; + } + } + } + } + FV::communicateFluxes(result); + return result; +} diff --git a/src/electromagnetic.cxx b/src/electromagnetic.cxx index de228b95..7048b47d 100644 --- a/src/electromagnetic.cxx +++ b/src/electromagnetic.cxx @@ -20,23 +20,54 @@ Electromagnetic::Electromagnetic(std::string name, Options &alloptions, Solver*) auto& options = alloptions[name]; + // Use the "Naulin" solver because we need to include toroidal + // variations of the density (A coefficient) + if (!options["laplacian"].isSet("type")) { + options["laplacian"]["type"] = "naulin"; + } aparSolver = Laplacian::create(&options["laplacian"]); - if (options["apar_boundary_neumann"] - .doc("Neumann radial boundaries? False => Zero Laplace") + const_gradient = options["const_gradient"] + .doc("Extrapolate gradient of Apar into all radial boundaries?") + .withDefault(false); + + // Give Apar an initial value because we solve Apar by iteration + // starting from the previous solution + Apar = 0.0; + + if (const_gradient) { + // Set flags to take the gradient from the RHS + aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_RHS); + aparSolver->setOuterBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_RHS); + last_time = 0.0; + + apar_boundary_timescale = options["apar_boundary_timescale"] + .doc("Timescale for Apar boundary relaxation [seconds]") + .withDefault(1e-8) + / get(alloptions["units"]["seconds"]); + + } else if (options["apar_boundary_neumann"] + .doc("Neumann on all radial boundaries?") .withDefault(false)) { - // Set zero-gradient (neumann) boundary conditions + // Set zero-gradient (neumann) boundary condition DC on the core aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); aparSolver->setOuterBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); - } else { - // Laplacian = 0 boundary conditions - aparSolver->setInnerBoundaryFlags(INVERT_DC_LAP + INVERT_AC_LAP); - aparSolver->setOuterBoundaryFlags(INVERT_DC_LAP + INVERT_AC_LAP); + + } else if (options["apar_core_neumann"] + .doc("Neumann radial boundary in the core? False => Dirichlet") + .withDefault(true) + and bout::globals::mesh->periodicY(bout::globals::mesh->xstart)) { + // Set zero-gradient (neumann) boundary condition DC on the core + aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD); } diagnose = options["diagnose"] .doc("Output additional diagnostics?") .withDefault(false); + + magnetic_flutter = options["magnetic_flutter"] + .doc("Set magnetic flutter terms (Apar_flutter)?") + .withDefault(false); } void Electromagnetic::transform(Options &state) { @@ -75,7 +106,46 @@ void Electromagnetic::transform(Options &state) { // Invert Helmholtz equation for Apar aparSolver->setCoefA((-beta_em) * alpha_em); - Apar = aparSolver->solve((-beta_em) * Ajpar); + + if (const_gradient) { + // Set gradient boundary condition from gradient inside boundary + Field3D rhs = (-beta_em) * Ajpar; + + const auto* mesh = Apar.getMesh(); + const auto* coords = Apar.getCoordinates(); + + BoutReal time = get(state["time"]); + BoutReal weight = 1.0; + if (time > last_time) { + weight = exp((last_time - time) / apar_boundary_timescale); + } + last_time = time; + + if (mesh->firstX()) { + const int x = mesh->xstart - 1; + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { + rhs(x, y, z) = (weight * (Apar(x + 1, y, z) - Apar(x, y, z)) + + (1 - weight) * (Apar(x + 2, y, z) - Apar(x + 1, y, z))) / + (sqrt(coords->g_11(x, y)) * coords->dx(x, y)); + } + } + } + if (mesh->lastX()) { + const int x = mesh->xend + 1; + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { + rhs(x, y, z) = (weight * (Apar(x, y, z) - Apar(x - 1, y, z)) + + (1 - weight) * (Apar(x - 1, y, z) - Apar(x - 2, y, z))) / + sqrt(coords->g_11(x, y)) / coords->dx(x, y); + } + } + } + // Use previous value of Apar as initial guess + Apar = aparSolver->solve(rhs, Apar); + } else { + Apar = aparSolver->solve((-beta_em) * Ajpar, Apar); + } // Save in the state set(state["fields"]["Apar"], Apar); @@ -95,10 +165,10 @@ void Electromagnetic::transform(Options &state) { const Field3D N = GET_NOBOUNDARY(Field3D, species["density"]); Field3D nv = getNonFinal(species["momentum"]); - nv -= Z * DC(N) * Apar; + nv -= Z * N * Apar; // Note: velocity is momentum / (A * N) Field3D v = getNonFinal(species["velocity"]); - v -= (Z / A) * DC(N) * Apar / floor(N, 1e-5); + v -= (Z / A) * N * Apar / floor(N, 1e-5); // Need to update the guard cells bout::globals::mesh->communicate(nv, v); v.applyBoundary("dirichlet"); @@ -107,6 +177,36 @@ void Electromagnetic::transform(Options &state) { set(species["momentum"], nv); set(species["velocity"], v); } + + if (magnetic_flutter) { + // Magnetic flutter terms + Apar_flutter = Apar - DC(Apar); + + // Ensure that guard cells are communicated + Apar.getMesh()->communicate(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) + 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) + Vector3D delta_B = Curl(A); + + // 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 + } } void Electromagnetic::outputVars(Options &state) { @@ -126,6 +226,16 @@ void Electromagnetic::outputVars(Options &state) { {"long_name", "Parallel component of vector potential A"} }); + if (magnetic_flutter) { + set_with_attrs(state["Apar_flutter"], Apar_flutter, { + {"time_dimension", "t"}, + {"units", "T m"}, + {"conversion", Bnorm * rho_s0}, + {"standard_name", "b dot A"}, + {"long_name", "Vector potential A|| used in flutter terms"} + }); + } + if (diagnose) { set_with_attrs(state["Ajpar"], Ajpar, { {"time_dimension", "t"}, diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index 892627b4..33ed120b 100644 --- a/src/evolve_density.cxx +++ b/src/evolve_density.cxx @@ -124,8 +124,6 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv } } - - neumann_boundary_average_z = alloptions[std::string("N") + name]["neumann_boundary_average_z"] .doc("Apply neumann boundary with Z average?") .withDefault(false); @@ -232,6 +230,14 @@ void EvolveDensity::finally(const Options& state) { } ddt(N) -= FV::Div_par_mod(N, V, fastest_wave); + + if (state.isSection("fields") and state["fields"].isSet("Apar_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); + } } if (low_n_diffuse) { diff --git a/src/evolve_energy.cxx b/src/evolve_energy.cxx index 2432528e..a6ce0e64 100644 --- a/src/evolve_energy.cxx +++ b/src/evolve_energy.cxx @@ -257,6 +257,12 @@ void EvolveEnergy::finally(const Options& state) { } ddt(E) -= FV::Div_par_mod(E + P, V, fastest_wave); + + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + // Magnetic flutter term + const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); + ddt(E) -= Div_n_g_bxGrad_f_B_XZ(E + P, V, -Apar_flutter); + } } if (species.isSet("low_n_coeff")) { @@ -322,6 +328,20 @@ void EvolveEnergy::finally(const Options& state) { // Note: Flux through boundary turned off, because sheath heat flux // is calculated and removed separately ddt(E) += FV::Div_par_K_Grad_par(kappa_par, T, false); + + if (state.isSection("fields") and state["fields"].isSet("Apar_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)) + // The first term is already calculated above. + // Here we add the terms containing db + 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); + ddt(E) += Div_par(kappa_par * db_dot_T) + - Div_n_g_bxGrad_f_B_XZ(kappa_par, db_dot_T + b0_dot_T, Apar_flutter); + } } if (hyper_z > 0.) { diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index b693f2ed..9aa1289a 100644 --- a/src/evolve_momentum.cxx +++ b/src/evolve_momentum.cxx @@ -100,15 +100,27 @@ void EvolveMomentum::finally(const Options &state) { // Apply a floor to the density Field3D Nlim = floor(N, density_floor); + // Typical wave speed used for numerical diffusion + Field3D fastest_wave; + if (state.isSet("fastest_wave")) { + fastest_wave = get(state["fastest_wave"]); + } else { + Field3D T = get(species["temperature"]); + fastest_wave = sqrt(T / AA); + } + + // Parallel flow + V = get(species["velocity"]); + if (state.isSection("fields") and state["fields"].isSet("phi") and species.isSet("charge")) { - BoutReal Z = get(species["charge"]); + const BoutReal Z = get(species["charge"]); if (Z != 0.0) { // Electrostatic potential set and species has charge // -> include ExB flow and parallel force - Field3D phi = get(state["fields"]["phi"]); + const Field3D phi = get(state["fields"]["phi"]); ddt(NV) = -Div_n_bxGrad_f_B_XPPM(NV, phi, bndry_flux, poloidal_flows, true); // ExB drift @@ -116,6 +128,26 @@ void EvolveMomentum::finally(const Options &state) { // Parallel electric field // Force density = - Z N ∇ϕ ddt(NV) -= Z * N * Grad_par(phi); + + if (state["fields"].isSet("Apar")) { + // Include a correction term for electromagnetic simulations + const Field3D Apar = get(state["fields"]["Apar"]); + + const Field3D density_source = species.isSet("density_source") ? + get(species["density_source"]) + : zeroFrom(Apar); + + // This is Z * Apar * dn/dt, keeping just leading order terms + Field3D dndt = density_source + - FV::Div_par_mod(N, V, fastest_wave) + - Div_n_bxGrad_f_B_XPPM(N, phi, bndry_flux, poloidal_flows, true) + ; + if (low_n_diffuse_perp) { + dndt += Div_Perp_Lap_FV_Index(density_floor / floor(N, 1e-3 * density_floor), N, + bndry_flux); + } + ddt(NV) += Z * Apar * dndt; + } } else { ddt(NV) = 0.0; } @@ -123,30 +155,29 @@ void EvolveMomentum::finally(const Options &state) { ddt(NV) = 0.0; } - // Parallel flow - V = get(species["velocity"]); - - // Typical wave speed used for numerical diffusion - Field3D fastest_wave; - if (state.isSet("fastest_wave")) { - fastest_wave = get(state["fastest_wave"]); - } else { - Field3D T = get(species["temperature"]); - fastest_wave = sqrt(T / AA); - } - // Note: // - Density floor should be consistent with calculation of V // otherwise energy conservation is affected // - using the same operator as in density and pressure equations doesn't work ddt(NV) -= AA * FV::Div_par_fvv(Nlim, V, fastest_wave, fix_momentum_boundary_flux); - + // Parallel pressure gradient if (species.isSet("pressure")) { Field3D P = get(species["pressure"]); ddt(NV) -= Grad_par(P); } + if (state.isSection("fields") and state["fields"].isSet("Apar_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); + } + } + if (species.isSet("low_n_coeff")) { // Low density parallel diffusion Field3D low_n_coeff = get(species["low_n_coeff"]); @@ -175,6 +206,8 @@ void EvolveMomentum::finally(const Options &state) { // If N < density_floor then NV and NV_solver may differ // -> Add term to force NV_solver towards NV + // Note: This correction is calculated in transform() + // because NV may be modified by electromagnetic terms ddt(NV) += NV_err; // Scale time derivatives @@ -202,6 +235,8 @@ void EvolveMomentum::finally(const Options &state) { } // Restore NV to the value returned by the solver // so that restart files contain the correct values + // Note: Copy boundary condition so dump file has correct boundary. + NV_solver.setBoundaryTo(NV); NV = NV_solver; } diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index 37f5da43..390d2b4c 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -268,6 +268,13 @@ void EvolvePressure::finally(const Options& state) { ddt(P) += (2. / 3) * V * Grad_par(P); } + + if (state.isSection("fields") and state["fields"].isSet("Apar_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); + } } if (species.isSet("low_n_coeff")) { @@ -419,6 +426,20 @@ void EvolvePressure::finally(const Options& state) { // Note: Flux through boundary turned off, because sheath heat flux // is calculated and removed separately ddt(P) += (2. / 3) * FV::Div_par_K_Grad_par(kappa_par, T, false); + + if (state.isSection("fields") and state["fields"].isSet("Apar_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)) + // The first term is already calculated above. + // Here we add the terms containing db + 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); + 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)); + } } if (hyper_z > 0.) { diff --git a/src/sheath_boundary.cxx b/src/sheath_boundary.cxx index 74a5050d..9d416c61 100644 --- a/src/sheath_boundary.cxx +++ b/src/sheath_boundary.cxx @@ -348,6 +348,10 @@ void SheathBoundary::transform(Options &state) { BoutReal q = ((gamma_e - 1 - 1 / (electron_adiabatic - 1)) * tesheath - 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; + if (q > 0.0) { + // Note: This could happen if tesheath > 2*phisheath + q = 0.0; + } // Multiply by cell area to get power BoutReal flux = q * (coord->J[i] + coord->J[im]) @@ -410,7 +414,10 @@ void SheathBoundary::transform(Options &state) { BoutReal q = ((gamma_e - 1 - 1 / (electron_adiabatic - 1)) * tesheath - 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; - + if (q < 0.0) { + // Note: This could happen if tesheath > 2*phisheath + q = 0.0; + } // Multiply by cell area to get power BoutReal flux = q * (coord->J[i] + coord->J[ip]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 936f71e1..496f7296 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -161,6 +161,10 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { // the first time RHS function is called phi_boundary_last_update = -1.; + phi_core_averagey = options["phi_core_averagey"] + .doc("Average phi core boundary in Y?") + .withDefault(false) and mesh->periodicY(mesh->xstart); + phi_boundary_timescale = options["phi_boundary_timescale"] .doc("Timescale for phi boundary relaxation [seconds]") .withDefault(1e-4) @@ -217,6 +221,7 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { void Vorticity::transform(Options& state) { AUTO_TRACE(); + phi.name = "phi"; auto& fields = state["fields"]; // Set the boundary of phi. Both 2D and 3D fields are kept, though the 3D field @@ -270,12 +275,32 @@ void Vorticity::transform(Options& state) { 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++) { - BoutReal phivalue = 0.0; - for (int k = 0; k < mesh->LocalNz; k++) { - phivalue += phi(mesh->xstart, j, k); + 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 } - phivalue /= mesh->LocalNz; // Average in Z of point next to boundary // Old value of phi at boundary BoutReal oldvalue =