diff --git a/hermes-3.cxx b/hermes-3.cxx index 57e1a74c..c5919631 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -248,12 +248,13 @@ int Hermes::init(bool restarting) { return 0; } -int Hermes::rhs(BoutReal time) { +int Hermes::rhs(BoutReal time, bool linear) { // Need to reset the state, since fields may be modified in transform steps state = Options(); - + set(state["time"], time); - state["units"] = units; + state["units"] = units; + set(state["solver"]["linear"], linear); // Call all the components scheduler->transform(state); diff --git a/hermes-3.hxx b/hermes-3.hxx index 3d273c57..b34a44ee 100644 --- a/hermes-3.hxx +++ b/hermes-3.hxx @@ -33,7 +33,7 @@ public: virtual ~Hermes() {} protected: int init(bool restarting) override; - int rhs(BoutReal t) override; + int rhs(BoutReal t, bool linear) override; int precon(BoutReal t, BoutReal gamma, BoutReal delta); /// Add variables to be written to the output file diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 4e28790b..f260a094 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -43,13 +43,27 @@ private: Field3D Dnn; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators + Field3D eta_n; ///< Viscosity + Field3D kappa_n; ///< Thermal conductivity bool sheath_ydown, sheath_yup; BoutReal nn_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. - BoutReal flux_limit; ///< Diffusive flux limit - BoutReal diffusion_limit; ///< Maximum diffusion coefficient + bool flux_limit; ///< Impose flux limiter? + bool particle_flux_limiter, heat_flux_limiter, momentum_flux_limiter; ///< Which limiters to impose + BoutReal maximum_mfp; ///< Maximum mean free path for diffusion. 0.1 by default, -1 is off. + BoutReal flux_limit_alpha; + BoutReal flux_limit_gamma; + Field3D particle_flux_factor; ///< Particle flux scaling factor + Field3D momentum_flux_factor; + Field3D heat_flux_factor; + + BoutReal flux_factor_timescale; ///< Timescale over which flux factors vary + BoutReal flux_factor_time {-1.0}; ///< Last time the flux factors were updated + Field3D particle_flux_avg; ///< Moving average flux factor + Field3D momentum_flux_avg; ///< Moving average flux factor + Field3D heat_flux_avg; ///< Moving average flux factor bool neutral_viscosity; ///< include viscosity? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 232ba75a..db64b818 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -51,18 +51,41 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .withDefault(true); flux_limit = options["flux_limit"] - .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") - .withDefault(0.2); + .doc("Use isotropic flux limiters?") + .withDefault(true); - diffusion_limit = options["diffusion_limit"] - .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") - .withDefault(-1.0) - / (meters * meters / seconds); // Normalise + flux_factor_timescale = options["flux_factor_timescale"] + .doc("Time constant in flux limitation factor variation [seconds]. < 0 is off.") + .withDefault(1e-6) / seconds; + + particle_flux_limiter = options["particle_flux_limiter"] + .doc("Enable particle flux limiter?") + .withDefault(true); + + heat_flux_limiter = options["heat_flux_limiter"] + .doc("Enable heat flux limiter?") + .withDefault(true); + + momentum_flux_limiter = options["momentum_flux_limiter"] + .doc("Enable momentum flux limiter?") + .withDefault(true); + + flux_limit_alpha = options["flux_limit_alpha"] + .doc("Scale flux limits") + .withDefault(1.0); + + flux_limit_gamma = options["flux_limit_gamma"] + .doc("Higher values increase sharpness of flux limiting") + .withDefault(2.0); neutral_viscosity = options["neutral_viscosity"] .doc("Include neutral gas viscosity?") .withDefault(true); + maximum_mfp = options["maximum_mfp"] + .doc("Optional maximum mean free path in [m] for diffusive processes. < 0 is off") + .withDefault(-1.0); + if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -156,6 +179,8 @@ void NeutralMixed::transform(Options& state) { Pnlim = floor(Nnlim * Tn, 1e-8); Pnlim.applyBoundary(); + Tnlim = Pnlim / Nnlim; + ///////////////////////////////////////////////////// // Parallel boundary conditions TRACE("Neutral boundary conditions"); @@ -245,69 +270,379 @@ void NeutralMixed::finally(const Options& state) { logPnlim = log(Pnlim); logPnlim.applyBoundary(); - /////////////////////////////////////////////////////// - // Calculate cross-field diffusion from collision frequency - // - // - BoutReal neutral_lmax = - 0.1 / get(state["units"]["meters"]); // Normalised length - - Field3D Rnn = sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] + // Note: If linear solve then don't recalculate coefficients + if (!get(state["solver"]["linear"])) { + /////////////////////////////////////////////////////// + // Calculate cross-field diffusion from collision frequency + // + // + if (localstate.isSet("collision_frequency")) { + // Dnn = Vth^2 / sigma + + if (maximum_mfp > 0) { // MFP limit enabled + Field3D Rnn = sqrt(Tn / AA) / (maximum_mfp / get(state["units"]["meters"])); + Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); + } else { // MFP limit disabled + Dnn = (Tn / AA) / (get(localstate["collision_frequency"])); + } - if (localstate.isSet("collision_frequency")) { - // Dnn = Vth^2 / sigma - Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); - } else { - Dnn = (Tn / AA) / Rnn; - } + } else { // If no collisions, hardcode max MFP to 0.1m + output_warn.write("No collisions set for the neutrals, limiting mean free path to 0.1m"); + Field3D Rnn = sqrt(Tn / AA) / (0.1 / get(state["units"]["meters"])); + Dnn = (Tn / AA) / Rnn; + } - if (flux_limit > 0.0) { - // Apply flux limit to diffusion, - // using the local thermal speed and pressure gradient magnitude - Field3D Dmax = flux_limit * sqrt(Tn / AA) / - (abs(Grad(logPnlim)) + 1. / neutral_lmax); - BOUT_FOR(i, Dmax.getRegion("RGN_NOBNDRY")) { - Dnn[i] = BOUTMIN(Dnn[i], Dmax[i]); + mesh->communicate(Dnn); + Dnn.clearParallelSlices(); + Dnn.applyBoundary(); + + // Neutral diffusion parameters have the same boundary condition as Dnn + DnnPn = Dnn * Pn; + DnnPn.applyBoundary(); + DnnNn = Dnn * Nn; + DnnNn.applyBoundary(); + DnnNVn = Dnn * NVn; + DnnNVn.applyBoundary(); + + if (sheath_ydown) { + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + Dnn(r.ind, mesh->ystart - 1, jz) = -Dnn(r.ind, mesh->ystart, jz); + DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); + DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); + DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); + } + } } - } - if (diffusion_limit > 0.0) { - // Impose an upper limit on the diffusion coefficient - BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { - Dnn[i] = BOUTMIN(Dnn[i], diffusion_limit); + if (sheath_yup) { + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + Dnn(r.ind, mesh->yend + 1, jz) = -Dnn(r.ind, mesh->yend, jz); + DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); + DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); + DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); + } + } } - } - mesh->communicate(Dnn); - Dnn.clearParallelSlices(); - Dnn.applyBoundary(); + // Heat conductivity + // Note: This is kappa_n = (5/2) * Pn / (m * nu) + // where nu is the collision frequency used in Dnn + kappa_n = (5. / 2) * DnnNn; + + // Viscosity + eta_n = AA * (2. / 5) * kappa_n; + + // Set factors that multiply the fluxes + particle_flux_factor = 1.0; + momentum_flux_factor = 1.0; + heat_flux_factor = 1.0; + if (flux_limit) { + // Apply flux limiters + // Note: Fluxes calculated here are cell centre, rather than cell edge + + // Cross-field velocity + Vector3D v_perp = -Dnn * Grad_perp(logPnlim); + + // Total velocity, perpendicular + parallel + Vector3D v_total = v_perp; + ASSERT2(v_total.covariant == true); + auto* coord = mesh->getCoordinates(); + v_total.y = Vn * (coord->J * coord->Bxy); // Set parallel component + + Field3D v_sq = v_total * v_total; // v dot v + Field3D v_abs = sqrt(v_sq); // |v| + + // Magnitude of the particle flux + Field3D particle_flux_abs = Nnlim * v_abs; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux from i to i+1 + + BoutReal fout = 0.5 * (DnnNn(i, j, k) + DnnNn(i + 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i + 1, j, k))) + * (logPnlim(i + 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + + // Calculate flux from i to i-1 + BoutReal fin = 0.5 * (DnnNn(i, j, k) + DnnNn(i - 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i - 1, j, k))) + * (logPnlim(i - 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i - 1, j, k)); + + // Flux from j to j+1 + BoutReal fup = 0.5 * (DnnNn(i, j, k) + DnnNn(i, j + 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j + 1, k))) + * (logPnlim(i, j + 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j + 1, k)); + + // Flux from j to j-1 + BoutReal fdown = 0.5 * (DnnNn(i, j, k) + DnnNn(i, j - 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j - 1, k))) + * (logPnlim(i, j - 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)); + + // Choose the maximum local particle flux, + // using cell center and edge values + particle_flux_abs(i, j, k) = BOUTMAX(particle_flux_abs(i, j, k), + fabs(fin), + fabs(fout), + fabs(fup), + fabs(fdown)); + } + } + } - // Neutral diffusion parameters have the same boundary condition as Dnn - DnnPn = Dnn * Pn; - DnnPn.applyBoundary(); - DnnNn = Dnn * Nn; - DnnNn.applyBoundary(); - Field3D DnnNVn = Dnn * NVn; - DnnNVn.applyBoundary(); + // Normalised particle flux limit + Field3D particle_limit = Nnlim * 0.25 * sqrt(8 * Tnlim / (PI * AA)); - if (sheath_ydown) { - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Dnn(r.ind, mesh->ystart - 1, jz) = -Dnn(r.ind, mesh->ystart, jz); - DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); - DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); - DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); + // Particle flux reduction factor + if (particle_flux_limiter) { + particle_flux_factor = pow(1. + pow(particle_flux_abs / (flux_limit_alpha * particle_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + particle_flux_factor = 1.0; } - } - } - if (sheath_yup) { - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Dnn(r.ind, mesh->yend + 1, jz) = -Dnn(r.ind, mesh->yend, jz); - DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); - DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); - DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); + // Flux of parallel momentum + // Note: Flux-limited particle flux is used in calculating the momentum flux before the + // momentum limiter is applied. + // + // The resulting momentum_flux_factor is applied to the momentum advection and viscosity terms, + // and so also scales the advection of kinetic energy + Vector3D momentum_flux = particle_flux_factor * NVn * v_total; + if (neutral_viscosity) { + momentum_flux -= eta_n * Grad_perp(Vn); + } + Field3D momentum_flux_abs = sqrt(momentum_flux * momentum_flux); + Field3D momentum_limit = Pnlim; + + if (momentum_flux_limiter) { + momentum_flux_factor = pow(1. + pow(momentum_flux_abs / (flux_limit_alpha * momentum_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + momentum_flux_factor = 1.0; + } + // Flux of heat + // Note: + // - Limiting the heat flux, not energy flux + // - Advection and heat conduction are both vectors, and e.g + // could be in opposite directions. + // - Flux doesn't include compression term, or kinetic energy transport + // that is in the momentum equation + Vector3D heat_flux = (3./2) * Pn * particle_flux_factor * v_total + + Pn * particle_flux_factor * v_perp + - kappa_n * Grad(Tn); + + Field3D heat_flux_abs = sqrt(heat_flux * heat_flux); + + mesh->communicate(particle_flux_factor); + particle_flux_factor.applyBoundary("neumann"); + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate energy flux from i to i+1 + + BoutReal fout = + // Convection with particle flux limiter + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i + 1, j, k) * DnnPn(i + 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i + 1, j, k))) + * (logPnlim(i + 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i + 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i + 1, j, k))) + * (Tn(i + 1, j, k) - Tn(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + ; + + // Calculate flux from i to i-1 + BoutReal fin = + // Convection + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i - 1, j, k) * DnnPn(i - 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i - 1, j, k))) + * (logPnlim(i - 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i - 1, j, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i - 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i - 1, j, k))) + * (Tn(i - 1, j, k) - Tn(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + ; + + // Flux from j to j+1 + BoutReal fup = + // Convection + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i, j + 1, k) * DnnPn(i, j + 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j + 1, k))) + * (logPnlim(i, j + 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i, j + 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j + 1, k))) + * (Tn(i, j + 1, k) - Tn(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) + ; + + // Flux from j to j-1 + BoutReal fdown = + // Convection + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i, j - 1, k) * DnnPn(i, j - 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j - 1, k))) + * (logPnlim(i, j - 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i, j - 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j - 1, k))) + * (Tn(i, j - 1, k) - Tn(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)) + ; + + heat_flux_abs(i, j, k) = BOUTMAX(heat_flux_abs(i, j, k), + fabs(fin), + fabs(fout), + fabs(fup), + fabs(fdown)); + } + } + } + + Field3D heat_limit = Pnlim * sqrt(2. * Tnlim / (PI * AA)); + + if (heat_flux_limiter) { + heat_flux_factor = pow(1. + pow(heat_flux_abs / (flux_limit_alpha * heat_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + heat_flux_factor = 1.0; + } + // Communicate guard cells and apply boundary conditions + // because the flux factors will be differentiated + mesh->communicate(particle_flux_factor, momentum_flux_factor, heat_flux_factor); + + // Fluxes at the boundary are not calculated correctly, so + // use limiting factors from one cell inwards + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + particle_flux_factor(r.ind, mesh->ystart, jz) = + BOUTMIN(particle_flux_factor(r.ind, mesh->ystart, jz), + particle_flux_factor(r.ind, mesh->ystart + 1, jz)); + momentum_flux_factor(r.ind, mesh->ystart, jz) = + BOUTMIN(momentum_flux_factor(r.ind, mesh->ystart, jz), + momentum_flux_factor(r.ind, mesh->ystart + 1, jz)); + heat_flux_factor(r.ind, mesh->ystart, jz) = + BOUTMIN(heat_flux_factor(r.ind, mesh->ystart, jz), + heat_flux_factor(r.ind, mesh->ystart + 1, jz)); + + // Neumann boundary conditions + particle_flux_factor(r.ind, mesh->ystart - 1, jz) = particle_flux_factor(r.ind, mesh->ystart, jz); + momentum_flux_factor(r.ind, mesh->ystart - 1, jz) = momentum_flux_factor(r.ind, mesh->ystart, jz); + heat_flux_factor(r.ind, mesh->ystart - 1, jz) = heat_flux_factor(r.ind, mesh->ystart, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + particle_flux_factor(r.ind, mesh->yend, jz) = + BOUTMIN(particle_flux_factor(r.ind, mesh->yend, jz), + particle_flux_factor(r.ind, mesh->yend - 1, jz)); + momentum_flux_factor(r.ind, mesh->yend, jz) = + BOUTMIN(momentum_flux_factor(r.ind, mesh->yend, jz), + momentum_flux_factor(r.ind, mesh->yend - 1, jz)); + heat_flux_factor(r.ind, mesh->yend, jz) = + BOUTMIN(heat_flux_factor(r.ind, mesh->yend, jz), + heat_flux_factor(r.ind, mesh->yend - 1, jz)); + + // Neumann boundary conditions + particle_flux_factor(r.ind, mesh->yend + 1, jz) = particle_flux_factor(r.ind, mesh->yend, jz); + momentum_flux_factor(r.ind, mesh->yend + 1, jz) = momentum_flux_factor(r.ind, mesh->yend, jz); + heat_flux_factor(r.ind, mesh->yend + 1, jz) = heat_flux_factor(r.ind, mesh->yend, jz); + } + } + + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + particle_flux_factor(mesh->xstart, j, k) = + BOUTMIN(particle_flux_factor(mesh->xstart, j, k), + particle_flux_factor(mesh->xstart + 1, j, k)); + momentum_flux_factor(mesh->xstart, j, k) = + BOUTMIN(momentum_flux_factor(mesh->xstart, j, k), + momentum_flux_factor(mesh->xstart + 1, j, k)); + heat_flux_factor(mesh->xstart, j, k) = + BOUTMIN(heat_flux_factor(mesh->xstart, j, k), + heat_flux_factor(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++) { + particle_flux_factor(mesh->xend, j, k) = + BOUTMIN(particle_flux_factor(mesh->xend, j, k), + particle_flux_factor(mesh->xend - 1, j, k)); + momentum_flux_factor(mesh->xend, j, k) = + BOUTMIN(momentum_flux_factor(mesh->xend, j, k), + momentum_flux_factor(mesh->xend - 1, j, k)); + heat_flux_factor(mesh->xend, j, k) = + BOUTMIN(heat_flux_factor(mesh->xend, j, k), + heat_flux_factor(mesh->xend - 1, j, k)); + } + } + } + + particle_flux_factor.applyBoundary("neumann"); + momentum_flux_factor.applyBoundary("neumann"); + heat_flux_factor.applyBoundary("neumann"); + + if (flux_factor_timescale > 0) { + // Smooth flux limitation factors over time + // Suppress rapid changes in flux limitation factors, while retaining + // the same steady-state solution. + BoutReal time = get(state["time"]); + if (flux_factor_time < 0.0) { + // First time, so just copy values + particle_flux_avg = particle_flux_factor; + momentum_flux_avg = momentum_flux_factor; + heat_flux_avg = heat_flux_factor; + } else { + // Weight by a factor that depends on the elapsed time. + BoutReal weight = exp(-(time - flux_factor_time) / flux_factor_timescale); + + particle_flux_avg = weight * particle_flux_avg + (1. - weight) * particle_flux_factor; + momentum_flux_avg = weight * momentum_flux_avg + (1. - weight) * momentum_flux_factor; + heat_flux_avg = weight * heat_flux_avg + (1. - weight) * heat_flux_factor; + + particle_flux_factor = particle_flux_avg; + momentum_flux_factor = momentum_flux_avg; + heat_flux_factor = heat_flux_avg; + } + flux_factor_time = time; } } } @@ -318,9 +653,11 @@ void NeutralMixed::finally(const Options& state) { ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed) // Advection - + FV::Div_a_Grad_perp(DnnNn, logPnlim) // Perpendicular diffusion - ; + + // Note: Parallel and perpendicular flux scaled by limiter + ddt(Nn) = -FV::Div_par_mod(Nn * particle_flux_factor, Vn, sound_speed) // Advection + + FV::Div_a_Grad_perp(DnnNn * particle_flux_factor, logPnlim) // Perpendicular diffusion + ; Sn = density_source; // Save for possible output if (localstate.isSet("density_source")) { @@ -333,26 +670,11 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral momentum"); ddt(NVn) = - -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - - Grad_par(Pn) // Pressure gradient - + FV::Div_a_Grad_perp(DnnNVn, logPnlim) // Perpendicular diffusion - ; - - if (neutral_viscosity) { - // NOTE: The following viscosity terms are are not (yet) balanced - // by a viscous heating term - - // Relationship between heat conduction and viscosity for neutral - // gas Chapman, Cowling "The Mathematical Theory of Non-Uniform - // Gases", CUP 1952 Ferziger, Kaper "Mathematical Theory of - // Transport Processes in Gases", 1972 - // eta_n = (2. / 5) * kappa_n; - // - - ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn) // Perpendicular viscosity - + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity - ; - } + // Note: The second argument (Vn) is squared in this operator, so the limiters are applied to the first argument + -AA * FV::Div_par_fvv(Nnlim * particle_flux_factor * momentum_flux_factor, Vn, sound_speed) // Momentum flow + - Grad_par(Pn) // Pressure gradient. Not included in flux limit. + + FV::Div_a_Grad_perp(DnnNVn * particle_flux_factor * momentum_flux_factor, logPnlim) // Perpendicular diffusion + ; if (localstate.isSet("momentum_source")) { Snv = get(localstate["momentum_source"]); @@ -363,11 +685,11 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Advection - - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp(DnnPn, logPnlim) // Perpendicular diffusion - + FV::Div_a_Grad_perp(DnnNn, Tn) // Conduction - + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction + ddt(Pn) = -FV::Div_par_mod(Pn * particle_flux_factor * heat_flux_factor, Vn, sound_speed) // Advection + - (2. / 3) * Pn * Div_par(Vn) // Compression + + FV::Div_a_Grad_perp((5. / 3) * DnnPn * particle_flux_factor * heat_flux_factor, logPnlim) // Perpendicular advection: q = 5/2 p u_perp + + (2. / 3) * (FV::Div_a_Grad_perp(kappa_n * heat_flux_factor, Tn) // Perpendicular Conduction + + FV::Div_par_K_Grad_par(kappa_n * heat_flux_factor, Tn)) // Parallel conduction ; Sp = pressure_source; @@ -376,6 +698,15 @@ void NeutralMixed::finally(const Options& state) { } ddt(Pn) += Sp; + if (neutral_viscosity) { + Field3D momentum_source = FV::Div_a_Grad_perp(eta_n * momentum_flux_factor, Vn) // Perpendicular viscosity + + FV::Div_par_K_Grad_par(eta_n * momentum_flux_factor, Vn) // Parallel viscosity + ; + + ddt(NVn) += momentum_source; // Viscosity + ddt(Pn) -= (2. / 3) * Vn * momentum_source; // Viscous heating + } + BOUT_FOR(i, Pn.getRegion("RGN_ALL")) { if ((Pn[i] < 1e-9) && (ddt(Pn)[i] < 0.0)) { ddt(Pn)[i] = 0.0; @@ -414,6 +745,7 @@ void NeutralMixed::outputVars(Options& state) { 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"}, @@ -468,6 +800,7 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "temperature"}, {"long_name", name + " temperature"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("Dnn") + name], Dnn, {{"time_dimension", "t"}, {"units", "m^2/s"}, @@ -475,6 +808,21 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "diffusion coefficient"}, {"long_name", name + " diffusion coefficient"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("eta_") + name], eta_n, + {{"time_dimension", "t"}, + {"units", "Pa s"}, + {"conversion", SQ(rho_s0) * Omega_ci * SI::Mp * Nnorm}, + {"standard_name", "viscosity"}, + {"long_name", name + " viscosity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("kappa_") + name], kappa_n, + {{"time_dimension", "t"}, + {"units", "W / m / eV"}, + {"conversion", SI::qe * Nnorm * rho_s0 * Cs0}, + {"standard_name", "heat conductivity"}, + {"long_name", name + " heat conductivity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("SN") + name], Sn, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, @@ -513,6 +861,30 @@ void NeutralMixed::outputVars(Options& state) { {"species", name}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("particle_flux_factor_") + name], particle_flux_factor, + {{"time_dimension", "t"}, + {"units", ""}, + {"conversion", 1.0}, + {"standard_name", "flux factor"}, + {"long_name", name + " particle flux factor"}, + {"species", name}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("momentum_flux_factor_") + name], momentum_flux_factor, + {{"time_dimension", "t"}, + {"units", ""}, + {"conversion", 1.0}, + {"standard_name", "flux factor"}, + {"long_name", name + " momentum flux factor"}, + {"species", name}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("heat_flux_factor_") + name], heat_flux_factor, + {{"time_dimension", "t"}, + {"units", ""}, + {"conversion", 1.0}, + {"standard_name", "flux factor"}, + {"long_name", name + " heat flux factor"}, + {"species", name}, + {"source", "neutral_mixed"}}); } } @@ -524,7 +896,7 @@ void NeutralMixed::precon(const Options& state, BoutReal gamma) { // Neutral gas diffusion // Solve (1 - gamma*Dnn*Delp2)^{-1} - Field3D coef = -gamma * Dnn; + Field3D coef = -gamma * Dnn * particle_flux_factor; if (state.isSet("scale_timederivs")) { coef *= get(state["scale_timederivs"]);