From bcc57fd66561e76eb3b9857970860feaf1da6a5a Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 2 Aug 2023 18:09:51 +0100 Subject: [PATCH 01/18] Add neutral pump to recycling component - Pump specified in grid anywhere along SOL/PFR edge - Only one recycle multiplier per species for all pumps --- include/recycling.hxx | 14 ++++--- src/recycling.cxx | 98 +++++++++++++++++++++++++++++++++++-------- 2 files changed, 90 insertions(+), 22 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index bf7b2111..aedd5741 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -46,20 +46,24 @@ private: /// Flux multiplier (recycling fraction). /// Combination of recycling fraction and species change e.g h+ -> h2 results in 0.5 multiplier - BoutReal target_multiplier, sol_multiplier, pfr_multiplier; + BoutReal target_multiplier, sol_multiplier, pfr_multiplier, pump_multiplier; BoutReal target_energy, sol_energy, pfr_energy; ///< Energy of recycled particle (normalised to Tnorm) }; std::vector channels; // Recycling channels - bool target_recycle, sol_recycle, pfr_recycle; ///< Flags for enabling recycling in different regions + bool target_recycle, sol_recycle, pfr_recycle, neutral_pump; ///< Flags for enabling recycling in different regions bool diagnose; ///< Save additional post-processing variables? Field3D density_source, energy_source; ///< Recycling particle and energy sources for all locations + Field3D is_pump; ///< 1 = pump, 0 = no pump. Works only in SOL/PFR - Field3D target_recycle_density_source, target_recycle_energy_source; ///< Recycling particle and energy sources for target recycling only - Field3D sol_recycle_density_source, sol_recycle_energy_source; ///< Recycling particle and energy sources for edge recycling only - Field3D pfr_recycle_density_source, pfr_recycle_energy_source; ///< Recycling particle and energy sources for edge recycling only + // Recycling particle and energy sources for the different sources of recycling + // Note that SOL, PFR and pump are not applicable to 1D + Field3D target_recycle_density_source, target_recycle_energy_source; + Field3D sol_recycle_density_source, sol_recycle_energy_source; + Field3D pfr_recycle_density_source, pfr_recycle_energy_source; + Field3D pump_recycle_density_source, pump_recycle_energy_source; Field3D radial_particle_outflow, radial_energy_outflow; ///< Radial fluxes coming from evolve_density and evolve_pressure used in recycling calc diff --git a/src/recycling.cxx b/src/recycling.cxx index 51bccff7..2d48d5f1 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -21,7 +21,13 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { .as(), ','); - + + // Neutral pump + // Mark cells as having a pump by setting the Field3D is_pump to 1 in the grid file + // Works only on SOL and PFR edges, where it locally modifies the recycle multiplier to the pump albedo + is_pump = 0.0; + mesh->get(is_pump, std::string("is_pump")); + for (const auto& species : species_list) { std::string from = trim(species, " \t\r()"); // The species name in the list @@ -52,6 +58,12 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { .doc("Multiply the pfr recycled flux by this factor. Should be >=0 and <= 1") .withDefault(1.0); + BoutReal pump_recycle_multiplier = + from_options["pump_recycle_multiplier"] + .doc("Multiply the pump boundary recycling flux by this factor (like albedo). Should be >=0 and <= 1") + .withDefault(1.0); + + BoutReal target_recycle_energy = from_options["target_recycle_energy"] .doc("Fixed energy of the recycled particles at target [eV]") .withDefault(3.0) @@ -69,14 +81,15 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { if ((target_recycle_multiplier < 0.0) or (target_recycle_multiplier > 1.0) or (sol_recycle_multiplier < 0.0) or (sol_recycle_multiplier > 1.0) - or (pfr_recycle_multiplier < 0.0) or (pfr_recycle_multiplier > 1.0)) { - throw BoutException("recycle_fraction must be betweeen 0 and 1"); + or (pfr_recycle_multiplier < 0.0) or (pfr_recycle_multiplier > 1.0) + or (pump_recycle_multiplier < 0.0) or (pump_recycle_multiplier > 1.0)) { + throw BoutException("All recycle multipliers must be betweeen 0 and 1"); } // Populate recycling channel vector channels.push_back({ from, to, - target_recycle_multiplier, sol_recycle_multiplier, pfr_recycle_multiplier, + target_recycle_multiplier, sol_recycle_multiplier, pfr_recycle_multiplier, pump_recycle_multiplier, target_recycle_energy, sol_recycle_energy, pfr_recycle_energy}); // Boolean flags for enabling recycling in different regions @@ -91,6 +104,10 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { pfr_recycle = from_options["pfr_recycle"] .doc("Recycling in the PFR edge?") .withDefault(false); + + neutral_pump = from_options["neutral_pump"] + .doc("Neutral pump enabled? Note, need location in grid file") + .withDefault(false); } } @@ -215,6 +232,12 @@ void Recycling::transform(Options& state) { BoutReal volume = J(mesh->xend, iy) * dx(mesh->xend, iy) * dy(mesh->xend, iy) * dz(mesh->xend, iy); + // If cell is a pump, overwrite multiplier with pump multiplier + BoutReal multiplier = channel.pfr_multiplier; + if ((is_pump(mesh->xend, iy, iz) == 1.0) and (neutral_pump)) { + multiplier = channel.pump_multiplier; + } + // Flow of recycled species back from the edge // SOL edge = LHS flow of inner guard cells on the high X side (mesh->xend+1) // Recycling source is 0 for each cell where the flow goes into instead of out of the domain @@ -224,14 +247,21 @@ void Recycling::transform(Options& state) { } // Divide by volume to get source - sol_recycle_density_source(mesh->xend, iy, iz) += recycle_particle_flow / volume; - density_source(mesh->xend, iy, iz) += sol_recycle_density_source(mesh->xend, iy, iz); - - // For now, this is a fixed temperature - sol_recycle_energy_source(mesh->xend, iy, iz) += channel.sol_energy * recycle_particle_flow / volume; - energy_source(mesh->xend, iy, iz) += sol_recycle_energy_source(mesh->xend, iy, iz); + BoutReal recycle_source = recycle_particle_flow / volume; + + // Add to appropriate diagnostic field depending if pump or not + if ((is_pump(mesh->xend+1, iy, iz) == 1.0) and (neutral_pump)) { + pump_recycle_density_source(mesh->xend+1, iy, iz) += recycle_source; + pump_recycle_energy_source(mesh->xend+1, iy, iz) += recycle_source * channel.pfr_energy; + } else { + pfr_recycle_density_source(mesh->xend+1, iy, iz) += recycle_source; + pfr_recycle_energy_source(mesh->xend+1, iy, iz) += recycle_source * channel.pfr_energy; + } - + // Add to density source which will be picked up by evolve_density.cxx + // Add to energy source which will be picked up by evolve_pressure.cxx + density_source(mesh->xend+1, iy, iz) += recycle_source; + energy_source(mesh->xend+1, iy, iz) += recycle_source * channel.pfr_energy; } } @@ -257,21 +287,36 @@ void Recycling::transform(Options& state) { BoutReal volume = J(mesh->xstart, iy) * dx(mesh->xstart, iy) * dy(mesh->xstart, iy) * dz(mesh->xstart, iy); + // If cell is a pump, overwrite multiplier with pump multiplier + BoutReal multiplier = channel.pfr_multiplier; + if ((is_pump(mesh->xstart, iy, iz) == 1.0) and (neutral_pump)) { + multiplier = channel.pump_multiplier; + } + // Flow of recycled species back from the edge // PFR edge = LHS flow of the first domain cell on the low X side (mesh->xstart) // Recycling source is 0 for each cell where the flow goes into instead of out of the domain BoutReal recycle_particle_flow = 0; if (radial_particle_outflow(mesh->xstart, iy, iz) > 0) { - recycle_particle_flow = channel.pfr_multiplier * radial_particle_outflow(mesh->xstart, iy, iz); + recycle_particle_flow = multiplier * radial_particle_outflow(mesh->xstart, iy, iz); } // Divide by volume to get source - pfr_recycle_density_source(mesh->xstart, iy, iz) += recycle_particle_flow / volume; - density_source(mesh->xstart, iy, iz) += pfr_recycle_density_source(mesh->xstart, iy, iz); + BoutReal recycle_source = recycle_particle_flow / volume; + + // Add to appropriate diagnostic field depending if pump or not + if ((is_pump(mesh->xstart, iy, iz) == 1.0) and (neutral_pump)) { + pump_recycle_density_source(mesh->xstart, iy, iz) += recycle_source; + pump_recycle_energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; + } else { + pfr_recycle_density_source(mesh->xstart, iy, iz) += recycle_source; + pfr_recycle_energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; + } - // For now, this is a fixed temperature - pfr_recycle_energy_source(mesh->xstart, iy, iz) += channel.pfr_energy * recycle_particle_flow / volume; - energy_source(mesh->xstart, iy, iz) += pfr_recycle_energy_source(mesh->xstart, iy, iz); + // Add to density source which will be picked up by evolve_density.cxx + // Add to energy source which will be picked up by evolve_pressure.cxx + density_source(mesh->xstart, iy, iz) += recycle_source; + energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; } } @@ -357,6 +402,25 @@ void Recycling::outputVars(Options& state) { {"long_name", std::string("PFR recycling energy source of ") + channel.to}, {"source", "recycling"}}); } + + // Neutral pump + if (neutral_pump) { + set_with_attrs(state[{std::string("S") + channel.to + std::string("_pump")}], pump_recycle_density_source, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "particle source"}, + {"long_name", std::string("Pump recycling particle source of ") + channel.to}, + {"source", "recycling"}}); + + set_with_attrs(state[{std::string("E") + channel.to + std::string("_pump")}], pump_recycle_energy_source, + {{"time_dimension", "t"}, + {"units", "W m^-3"}, + {"conversion", Pnorm * Omega_ci}, + {"standard_name", "energy source"}, + {"long_name", std::string("Pump recycling energy source of ") + channel.to}, + {"source", "recycling"}}); + } } } From 06099fbd610baea72d3d7e0502a6d50c2733877a Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 7 Aug 2023 14:41:31 +0100 Subject: [PATCH 02/18] Several neutral pump fixes --- include/recycling.hxx | 2 +- src/recycling.cxx | 55 +++++++++++++++++++++++-------------------- 2 files changed, 31 insertions(+), 26 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index aedd5741..b56bbf82 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -56,7 +56,7 @@ private: bool diagnose; ///< Save additional post-processing variables? Field3D density_source, energy_source; ///< Recycling particle and energy sources for all locations - Field3D is_pump; ///< 1 = pump, 0 = no pump. Works only in SOL/PFR + Field2D is_pump; ///< 1 = pump, 0 = no pump. Works only in SOL/PFR // Recycling particle and energy sources for the different sources of recycling // Note that SOL, PFR and pump are not applicable to 1D diff --git a/src/recycling.cxx b/src/recycling.cxx index 2d48d5f1..d50516ad 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -215,6 +215,10 @@ void Recycling::transform(Options& state) { } } + // Initialise counters of pump recycling fluxes + pump_recycle_density_source = 0; + pump_recycle_energy_source = 0; + // Recycling at the SOL edge (2D/3D only) if (sol_recycle) { @@ -223,6 +227,7 @@ void Recycling::transform(Options& state) { sol_recycle_density_source = 0; sol_recycle_energy_source = 0; + if(mesh->lastX()){ // Only do this for the processor which has the edge region for(int iy=0; iy < mesh->LocalNy ; iy++){ @@ -250,12 +255,12 @@ void Recycling::transform(Options& state) { BoutReal recycle_source = recycle_particle_flow / volume; // Add to appropriate diagnostic field depending if pump or not - if ((is_pump(mesh->xend+1, iy, iz) == 1.0) and (neutral_pump)) { + if ((is_pump(mesh->xend+1, iy) == 1.0) and (neutral_pump)) { pump_recycle_density_source(mesh->xend+1, iy, iz) += recycle_source; pump_recycle_energy_source(mesh->xend+1, iy, iz) += recycle_source * channel.pfr_energy; } else { - pfr_recycle_density_source(mesh->xend+1, iy, iz) += recycle_source; - pfr_recycle_energy_source(mesh->xend+1, iy, iz) += recycle_source * channel.pfr_energy; + sol_recycle_density_source(mesh->xend+1, iy, iz) += recycle_source; + sol_recycle_energy_source(mesh->xend+1, iy, iz) += recycle_source * channel.pfr_energy; } // Add to density source which will be picked up by evolve_density.cxx @@ -305,14 +310,14 @@ void Recycling::transform(Options& state) { BoutReal recycle_source = recycle_particle_flow / volume; // Add to appropriate diagnostic field depending if pump or not - if ((is_pump(mesh->xstart, iy, iz) == 1.0) and (neutral_pump)) { + if ((is_pump(mesh->xstart, iy) == 1.0) and (neutral_pump)) { pump_recycle_density_source(mesh->xstart, iy, iz) += recycle_source; pump_recycle_energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; } else { pfr_recycle_density_source(mesh->xstart, iy, iz) += recycle_source; pfr_recycle_energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; } - + // Add to density source which will be picked up by evolve_density.cxx // Add to energy source which will be picked up by evolve_pressure.cxx density_source(mesh->xstart, iy, iz) += recycle_source; @@ -331,7 +336,7 @@ void Recycling::transform(Options& state) { } void Recycling::outputVars(Options& state) { - + AUTO_TRACE(); // Normalisations auto Nnorm = get(state["Nnorm"]); @@ -401,26 +406,26 @@ void Recycling::outputVars(Options& state) { {"standard_name", "energy source"}, {"long_name", std::string("PFR recycling energy source of ") + channel.to}, {"source", "recycling"}}); - } + } - // Neutral pump - if (neutral_pump) { - set_with_attrs(state[{std::string("S") + channel.to + std::string("_pump")}], pump_recycle_density_source, - {{"time_dimension", "t"}, - {"units", "m^-3 s^-1"}, - {"conversion", Nnorm * Omega_ci}, - {"standard_name", "particle source"}, - {"long_name", std::string("Pump recycling particle source of ") + channel.to}, - {"source", "recycling"}}); - - set_with_attrs(state[{std::string("E") + channel.to + std::string("_pump")}], pump_recycle_energy_source, - {{"time_dimension", "t"}, - {"units", "W m^-3"}, - {"conversion", Pnorm * Omega_ci}, - {"standard_name", "energy source"}, - {"long_name", std::string("Pump recycling energy source of ") + channel.to}, - {"source", "recycling"}}); - } + // Neutral pump + if (neutral_pump) { + set_with_attrs(state[{std::string("S") + channel.to + std::string("_pump_recycle")}], pump_recycle_density_source, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "particle source"}, + {"long_name", std::string("Pump recycling particle source of ") + channel.to}, + {"source", "recycling"}}); + + set_with_attrs(state[{std::string("E") + channel.to + std::string("_pump_recycle")}], pump_recycle_energy_source, + {{"time_dimension", "t"}, + {"units", "W m^-3"}, + {"conversion", Pnorm * Omega_ci}, + {"standard_name", "energy source"}, + {"long_name", std::string("Pump recycling energy source of ") + channel.to}, + {"source", "recycling"}}); + } } } From ca94b6d87a859b1222363598745cdd7138a09770 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 16 Aug 2023 15:36:23 -0700 Subject: [PATCH 03/18] ion_viscosity: Only keep axisymmetric part of perpendicular viscosity Removing non-axisymmetric components, so we just include the poloidal flow damping. --- src/ion_viscosity.cxx | 74 ++++++++++++++++++++----------------------- 1 file changed, 35 insertions(+), 39 deletions(-) diff --git a/src/ion_viscosity.cxx b/src/ion_viscosity.cxx index fbb773d1..64232cd8 100644 --- a/src/ion_viscosity.cxx +++ b/src/ion_viscosity.cxx @@ -67,9 +67,9 @@ void IonViscosity::transform(Options &state) { Options& allspecies = state["species"]; auto coord = mesh->getCoordinates(); - const Field3D Bxy = coord->Bxy; - const Field3D sqrtB = sqrt(Bxy); - const Field3D Grad_par_logB = Grad_par(log(Bxy)); + const Field2D Bxy = coord->Bxy; + const Field2D sqrtB = sqrt(Bxy); + const Field2D Grad_par_logB = Grad_par(log(Bxy)); // Loop through all species for (auto& kv : allspecies.getChildren()) { @@ -123,35 +123,40 @@ void IonViscosity::transform(Options &state) { continue; // Skip perpendicular flow parts below } + // The following calculation is performed on the axisymmetric components + // Need electrostatic potential for ExB flow - const Field3D phi = get(state["fields"]["phi"]); + const Field2D phi_av = DC(get(state["fields"]["phi"])); // Density and temperature needed for diamagnetic flows - const Field3D N = get(species["density"]); - const Field3D T = get(species["temperature"]); + const Field2D N_av = DC(get(species["density"])); + const Field2D T_av = DC(get(species["temperature"])); + const Field2D P_av = DC(P); + const Field2D tau_av = DC(tau); + const Field2D V_av = DC(V); // Parallel ion stress tensor component - Field3D Pi_cipar = -0.96 * P * tau * - (2. * Grad_par(V) + V * Grad_par_logB); + Field2D Pi_cipar = -0.96 * P_av * tau_av * + (2. * Grad_par(V_av) + V_av * Grad_par_logB); // Could also be written as: // Pi_cipar = -0.96*Pi*tau*2.*Grad_par(sqrt(Bxy)*Vi)/sqrt(Bxy); // Perpendicular ion stress tensor // 0.96 P tau kappa * (V_E + V_di + 1.61 b x Grad(T)/B ) // Note: Heat flux terms are neglected for now - Field3D Pi_ciperp = -0.5 * 0.96 * P * tau * - (Curlb_B * Grad(phi + 1.61 * T) - Curlb_B * Grad(P) / N); + Field2D Pi_ciperp = -0.5 * 0.96 * P_av * tau_av * + (Curlb_B * Grad(phi_av + 1.61 * T_av) - Curlb_B * Grad(P_av) / N_av); // Limit size of stress tensor components // If the off-diagonal components of the pressure tensor are large compared // to the scalar pressure, then the model is probably breaking down. // This can happen in very low density regions BOUT_FOR(i, Pi_ciperp.getRegion("RGN_NOBNDRY")) { - if (fabs(Pi_ciperp[i]) > P[i]) { - Pi_ciperp[i] = SIGN(Pi_ciperp[i]) * P[i]; + if (fabs(Pi_ciperp[i]) > P_av[i]) { + Pi_ciperp[i] = SIGN(Pi_ciperp[i]) * P_av[i]; } - if (fabs(Pi_cipar[i]) > P[i]) { - Pi_cipar[i] = SIGN(Pi_cipar[i]) * P[i]; + if (fabs(Pi_cipar[i]) > P_av[i]) { + Pi_cipar[i] = SIGN(Pi_cipar[i]) * P_av[i]; } } @@ -160,26 +165,16 @@ void IonViscosity::transform(Options &state) { Pi_cipar.applyBoundary("neumann"); // Apply parallel boundary conditions - Pi_ciperp = toFieldAligned(Pi_ciperp); - Pi_cipar = toFieldAligned(Pi_cipar); - { - int jy = 1; - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Pi_ciperp(r.ind, jy, jz) = Pi_ciperp(r.ind, jy + 1, jz); - Pi_cipar(r.ind, jy, jz) = Pi_cipar(r.ind, jy + 1, jz); - } - } - jy = mesh->yend + 1; - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Pi_ciperp(r.ind, jy, jz) = Pi_ciperp(r.ind, jy - 1, jz); - Pi_cipar(r.ind, jy, jz) = Pi_cipar(r.ind, jy - 1, jz); - } - } + int jy = 1; + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + Pi_ciperp(r.ind, jy) = Pi_ciperp(r.ind, jy + 1); + Pi_cipar(r.ind, jy) = Pi_cipar(r.ind, jy + 1); + } + jy = mesh->yend + 1; + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + Pi_ciperp(r.ind, jy) = Pi_ciperp(r.ind, jy - 1); + Pi_cipar(r.ind, jy) = Pi_cipar(r.ind, jy - 1); } - Pi_ciperp = fromFieldAligned(Pi_ciperp); - Pi_cipar = fromFieldAligned(Pi_cipar); mesh->communicate(Pi_ciperp, Pi_cipar); @@ -187,13 +182,13 @@ void IonViscosity::transform(Options &state) { //const Field3D div_Pi_ciperp = - (2. / 3) * B32 * Grad_par(Pi_ciperp / B32); add(species["momentum_source"], div_Pi_ciperp); - subtract(species["energy_source"], V * div_Pi_ciperp); + subtract(species["energy_source"], V_av * div_Pi_ciperp); // Total scalar ion viscous pressure - Field3D Pi_ci = Pi_cipar + Pi_ciperp; + Field2D Pi_ci = Pi_cipar + Pi_ciperp; #if CHECKLEVEL >= 1 - for (auto& i : N.getRegion("RGN_NOBNDRY")) { + for (auto& i : N_av.getRegion("RGN_NOBNDRY")) { if (!std::isfinite(Pi_cipar[i])) { throw BoutException("{} Pi_cipar non-finite at {}.\n", species_name, i); } @@ -205,12 +200,13 @@ void IonViscosity::transform(Options &state) { // Divergence of current in vorticity equation Field3D DivJ = Div(0.5 * Pi_ci * Curlb_B) - - Div_n_bxGrad_f_B_XPPM(1. / 3, Pi_ci, false); + Div_n_bxGrad_f_B_XPPM(1. / 3, Pi_ci, false, true); add(state["fields"]["DivJextra"], DivJ); // Transfer of energy between ion internal energy and ExB flow - subtract(species["energy_source"], 0.5 * Pi_ci * Curlb_B * Grad(phi + P) - - (1 / 3) * bracket(Pi_ci, phi + P, BRACKET_ARAKAWA)); + subtract(species["energy_source"], + Field3D(0.5 * Pi_ci * Curlb_B * Grad(phi_av + P_av) + - (1 / 3) * bracket(Pi_ci, phi_av + P_av, BRACKET_STD))); if (diagnose) { // Find the diagnostics struct for this species From 60486b7331fa8a1ef11a3066d8173bc9220539ff Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 20 Aug 2023 14:37:03 -0700 Subject: [PATCH 04/18] vorticity: Fix diamagnetic_polarisation energy transfer Should include both parallel and diamagnetic current, but only diamagnetic term was included in the transfer of energy to/from ion pressure. Unfortunately calculating this term depends on current, that depends on sheath boundaries, that depends on potential. It therefore can't be calculated in `vorticity` unless assumptions are made on parallel current at the boundary (e.g. zero-gradient). Instead, these terms are now calculated in the `polarisation_drift` component, because they are connected to the divergence of the polarisation drift velocity. Fix only affects simulations with vorticity equation and hot ions. --- include/polarisation_drift.hxx | 2 + src/polarisation_drift.cxx | 104 ++++++++++++++++++++++++--------- src/vorticity.cxx | 26 --------- 3 files changed, 79 insertions(+), 53 deletions(-) diff --git a/include/polarisation_drift.hxx b/include/polarisation_drift.hxx index 3ee2b8ea..0663ad9f 100644 --- a/include/polarisation_drift.hxx +++ b/include/polarisation_drift.hxx @@ -70,6 +70,8 @@ private: bool boussinesq; // If true, assume a constant mass density in Jpol BoutReal average_atomic_mass; // If boussinesq=true, mass density to use BoutReal density_floor; // Minimum mass density if boussinesq=false + bool advection; // Advect fluids by an approximate polarisation velocity? + bool diamagnetic_polarisation; // Calculate compression terms? }; namespace { diff --git a/src/polarisation_drift.cxx b/src/polarisation_drift.cxx index 4b11c640..5c49b450 100644 --- a/src/polarisation_drift.cxx +++ b/src/polarisation_drift.cxx @@ -44,7 +44,16 @@ PolarisationDrift::PolarisationDrift(std::string name, // Use a density floor to prevent divide-by-zero errors density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-5); } - + + advection = options["advection"] + .doc("Include advection by polarisation drift, using potential flow approximation?") + .withDefault(true); + + diamagnetic_polarisation = + options["diamagnetic_polarisation"] + .doc("Include diamagnetic drift in polarisation current?") + .withDefault(true); + diagnose = options["diagnose"] .doc("Output additional diagnostics?") .withDefault(false); @@ -56,6 +65,72 @@ void PolarisationDrift::transform(Options &state) { // Iterate through all subsections Options& allspecies = state["species"]; + // Calculate divergence of all currents except the polarisation current + + if (IS_SET(state["fields"]["DivJdia"])) { + DivJ = get(state["fields"]["DivJdia"]); + } else { + DivJ = 0.0; + } + + // Parallel current due to species parallel flow + for (auto& kv : allspecies.getChildren()) { + const Options& species = kv.second; + + if (!species.isSet("charge") or !species.isSet("momentum")) { + continue; // Not charged, or no parallel flow + } + const BoutReal Z = get(species["charge"]); + if (fabs(Z) < 1e-5) { + continue; // Not charged + } + + const Field3D NV = GET_VALUE(Field3D, species["momentum"]); + const BoutReal A = get(species["AA"]); + + // Note: Using NV rather than N*V so that the cell boundary flux is correct + DivJ += Div_par((Z / A) * NV); + } + + if (diamagnetic_polarisation) { + // Compression of ion diamagnetic contribution to polarisation velocity + // Note: For now this ONLY includes the parallel and diamagnetic current terms + // Other terms e.g. ion viscous current, are in their separate components + if (!boussinesq) { + throw BoutException("diamagnetic_polarisation not implemented for non-Boussinesq"); + } + + // Calculate energy exchange term nonlinear in pressure + // (3 / 2) ddt(Pi) += (Pi / n0) * Div((Pe + Pi) * Curlb_B + Jpar); + for (auto& kv : allspecies.getChildren()) { + Options& species = allspecies[kv.first]; // Note: need non-const + + if (!(IS_SET_NOBOUNDARY(species["pressure"]) and species.isSet("charge") + and species.isSet("AA"))) { + // No pressure, charge or mass -> no polarisation current due to + // diamagnetic flow + continue; + } + + const auto charge = get(species["charge"]); + if (fabs(charge) < 1e-5) { + // No charge + continue; + } + + const auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); + const auto AA = get(species["AA"]); + + add(species["energy_source"], + P * (AA / average_atomic_mass / charge) * DivJ); + } + } + + if (!advection) { + return; + } + // Calculate advection terms using a potential-flow approximation + // Calculate the total mass density of species // which contribute to polarisation current Field3D mass_density; @@ -81,39 +156,14 @@ void PolarisationDrift::transform(Options &state) { // Apply a floor to prevent divide-by-zero errors mass_density = floor(mass_density, density_floor); } - - // Calculate divergence of all currents except the polarisation current - DivJ = 0.0; - + if (IS_SET(state["fields"]["DivJextra"])) { DivJ += get(state["fields"]["DivJextra"]); } - if (IS_SET(state["fields"]["DivJdia"])) { - DivJ += get(state["fields"]["DivJdia"]); - } if (IS_SET(state["fields"]["DivJcol"])) { DivJ += get(state["fields"]["DivJcol"]); } - // Parallel current due to species parallel flow - for (auto& kv : allspecies.getChildren()) { - const Options& species = kv.second; - - if (!species.isSet("charge") or !species.isSet("momentum")) { - continue; // Not charged, or no parallel flow - } - const BoutReal Z = get(species["charge"]); - if (fabs(Z) < 1e-5) { - continue; // Not charged - } - - const Field3D NV = GET_VALUE(Field3D, species["momentum"]); - const BoutReal A = get(species["AA"]); - - // Note: Using NV rather than N*V so that the cell boundary flux is correct - DivJ += Div_par((Z / A) * NV); - } - // Solve for time derivative of potential // Using Div(mass_density / B^2 Grad_perp(dphi/dt)) = DivJ diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 5224f189..6b7951ae 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -500,32 +500,6 @@ void Vorticity::transform(Options& state) { DivJdia = Div(Jdia); ddt(Vort) += DivJdia; - if (diamagnetic_polarisation) { - // Calculate energy exchange term nonlinear in pressure - // (3 / 2) ddt(Pi) += Pi * Div((Pe + Pi) * Curlb_B); - for (auto& kv : allspecies.getChildren()) { - Options& species = allspecies[kv.first]; // Note: need non-const - - if (!(IS_SET_NOBOUNDARY(species["pressure"]) and species.isSet("charge") - and species.isSet("AA"))) { - continue; // No pressure, charge or mass -> no polarisation current due to - // diamagnetic flow - } - - const auto charge = get(species["charge"]); - if (fabs(charge) < 1e-5) { - // No charge - continue; - } - - const auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); - const auto AA = get(species["AA"]); - - add(species["energy_source"], - P * (AA / average_atomic_mass / charge) * DivJdia); - } - } - set(fields["DivJdia"], DivJdia); } From a2119abe31ec2ce5bde0f95704a3b5c9210ea27d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 27 Aug 2023 09:36:29 -0700 Subject: [PATCH 05/18] vorticity: Add option to damp vorticity near core Core boundary poloidal flow. By damping vorticity this is damping towards a Neumann condition on the ion flow. --- include/vorticity.hxx | 5 ++++- src/vorticity.cxx | 20 ++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/include/vorticity.hxx b/include/vorticity.hxx index 7731b187..18fd8b94 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -49,7 +49,9 @@ struct Vorticity : public Component { /// Kinematic viscosity [m^2/s] /// - vort_dissipation: bool, default false /// Parallel dissipation of vorticity? - /// + /// - damp_core_vorticity: bool, default false + /// Damp axisymmetric component of vorticity in cell next to core boundary + /// Vorticity(std::string name, Options &options, Solver *solver); /// Optional inputs @@ -121,6 +123,7 @@ private: bool vort_dissipation; ///< Parallel dissipation of vorticity bool phi_dissipation; ///< Parallel dissipation of potential bool phi_sheath_dissipation; ///< Dissipation at the sheath if phi < 0 + bool damp_core_vorticity; ///< Damp axisymmetric component of vorticity bool phi_boundary_relax; ///< Relax boundary to zero-gradient BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised] diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 6b7951ae..d2185d99 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -107,6 +107,10 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { .doc("Add dissipation when phi < 0.0 at the sheath") .withDefault(false); + damp_core_vorticity = options["damp_core_vorticity"] + .doc("Damp vorticity at the core boundary?") + .withDefault(false); + // Add phi to restart files so that the value in the boundaries // is restored on restart. This is done even when phi is not evolving, // so that phi can be saved and re-loaded @@ -655,6 +659,22 @@ void Vorticity::finally(const Options& state) { } ddt(Vort) += fromFieldAligned(dissipation); } + + if (damp_core_vorticity) { + // Damp axisymmetric vorticity near core boundary + if (mesh->firstX() and mesh->periodicY(mesh->xstart)) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal vort_avg = 0.0; // Average Vort in Z + for (int k = 0; k < mesh->LocalNz; k++) { + vort_avg += Vort(mesh->xstart, j, k); + } + vort_avg /= mesh->LocalNz; + for (int k = 0; k < mesh->LocalNz; k++) { + ddt(Vort)(mesh->xstart, j, k) -= 0.01 * vort_avg; + } + } + } + } } void Vorticity::outputVars(Options& state) { From 5a6bc51173256b5d7f6873eed1a90e70c62fee45 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 27 Aug 2023 20:06:50 +0100 Subject: [PATCH 06/18] Fix typo in neutral pump - Was not picking up correct recycle multiplier --- src/recycling.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index 06a468ac..185e79e5 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -246,7 +246,7 @@ void Recycling::transform(Options& state) { // Recycling source is 0 for each cell where the flow goes into instead of out of the domain BoutReal recycle_particle_flow = 0; if (radial_particle_outflow(mesh->xend+1, iy, iz) > 0) { - recycle_particle_flow = channel.sol_multiplier * radial_particle_outflow(mesh->xend+1, iy, iz); + recycle_particle_flow = multiplier * radial_particle_outflow(mesh->xend+1, iy, iz); } // Divide by volume to get source From 2c9695d9f865ed87ef9bcd3a5dd4d088abe2f8f7 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 27 Aug 2023 22:31:28 -0700 Subject: [PATCH 07/18] ion_viscosity: Add boundary condition and diagnostic - Set Neumann boundary condition on Pi_ci before taking Div() - Add a DivJvis output diagnostic, to inspect the divergence of current passed to the vorticity equation. --- include/ion_viscosity.hxx | 1 + src/ion_viscosity.cxx | 14 +++++++++++++- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/include/ion_viscosity.hxx b/include/ion_viscosity.hxx index 06fa0ad8..118caf42 100644 --- a/include/ion_viscosity.hxx +++ b/include/ion_viscosity.hxx @@ -61,6 +61,7 @@ private: struct Diagnostics { Field3D Pi_ciperp; ///< Perpendicular part of Pi scalar Field3D Pi_cipar; ///< Parallel part of Pi scalar + Field3D DivJ; ///< Divergence of current in vorticity equation }; /// Store diagnostics for each species diff --git a/src/ion_viscosity.cxx b/src/ion_viscosity.cxx index 64232cd8..6f253023 100644 --- a/src/ion_viscosity.cxx +++ b/src/ion_viscosity.cxx @@ -198,6 +198,8 @@ void IonViscosity::transform(Options &state) { } #endif + Pi_ci.applyBoundary("neumann"); + // Divergence of current in vorticity equation Field3D DivJ = Div(0.5 * Pi_ci * Curlb_B) - Div_n_bxGrad_f_B_XPPM(1. / 3, Pi_ci, false, true); @@ -213,12 +215,13 @@ void IonViscosity::transform(Options &state) { auto search = diagnostics.find(species_name); if (search == diagnostics.end()) { // First time, create diagnostic - diagnostics.emplace(species_name, Diagnostics {Pi_ciperp, Pi_cipar}); + diagnostics.emplace(species_name, Diagnostics {Pi_ciperp, Pi_cipar, DivJ}); } else { // Update diagnostic values auto& d = search->second; d.Pi_ciperp = Pi_ciperp; d.Pi_cipar = Pi_cipar; + d.DivJ = DivJ; } } } @@ -231,6 +234,7 @@ void IonViscosity::outputVars(Options &state) { // Normalisations auto Nnorm = get(state["Nnorm"]); auto Tnorm = get(state["Tnorm"]); + auto Omega_ci = get(state["Omega_ci"]); BoutReal Pnorm = SI::qe * Tnorm * Nnorm; // Pressure normalisation for (const auto& it : diagnostics) { @@ -254,6 +258,14 @@ void IonViscosity::outputVars(Options &state) { {"long_name", species_name + " parallel collisional viscous pressure"}, {"species", species_name}, {"source", "ion_viscosity"}}); + + set_with_attrs(state[std::string("DivJvis_") + species_name], d.DivJ, + {{"time_dimension", "t"}, + {"units", "A m^-3"}, + {"conversion", SI::qe * Nnorm * Omega_ci}, + {"long_name", std::string("Divergence of viscous current due to species") + species_name}, + {"species", species_name}, + {"source", "ion_viscosity"}}); } } } From 09ff0dc3efdd3bf38d96064bba17ca1fb16dabfb Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 27 Aug 2023 22:33:57 -0700 Subject: [PATCH 08/18] vorticity: Add pressure boundary when calculating DivJdia Diamagnetic current needs a boundary on P, but that is set in the sheath boundary condition. For now DivJdia is calculated before the sheath boundary, so have to apply the P boundary here too. --- src/vorticity.cxx | 58 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/src/vorticity.cxx b/src/vorticity.cxx index d2185d99..936f71e1 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -23,6 +23,32 @@ Ind3D indexAt(const Field3D& f, int x, int y, int z) { int nz = f.getNz(); return Ind3D{(x * ny + y) * nz + z, ny, nz}; } + +/// Limited free gradient of log of a quantity +/// This ensures that the guard cell values remain positive +/// while also ensuring that the quantity never increases +/// +/// fm fc | fp +/// ^ boundary +/// +/// exp( 2*log(fc) - log(fm) ) +/// +BoutReal limitFree(BoutReal fm, BoutReal fc) { + if (fm < fc) { + return fc; // Neumann rather than increasing into boundary + } + if (fm < 1e-10) { + return fc; // Low / no density condition + } + BoutReal fp = SQ(fc) / fm; +#if CHECKLEVEL >= 2 + if (!std::isfinite(fp)) { + throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp); + } +#endif + + return fp; +} } Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { @@ -460,6 +486,38 @@ void Vorticity::transform(Options& state) { auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); + // Note: We need boundary conditions on P, so apply the same + // free boundary condition as sheath_boundary. + if (P.hasParallelSlices()) { + Field3D &P_ydown = P.ydown(); + Field3D &P_yup = P.yup(); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + P_ydown(r.ind, mesh->ystart - 1, jz) = 2 * P(r.ind, mesh->ystart, jz) - P_yup(r.ind, mesh->ystart + 1, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + P_yup(r.ind, mesh->yend + 1, jz) = 2 * P(r.ind, mesh->yend, jz) - P_ydown(r.ind, mesh->yend - 1, jz); + } + } + } else { + Field3D P_fa = toFieldAligned(P); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + auto i = indexAt(P_fa, r.ind, mesh->ystart, jz); + P_fa[i.ym()] = limitFree(P_fa[i.yp()], P_fa[i]); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + auto i = indexAt(P_fa, r.ind, mesh->yend, jz); + P_fa[i.yp()] = limitFree(P_fa[i.ym()], P_fa[i]); + } + } + P = fromFieldAligned(P_fa); + } + // Note: This calculation requires phi derivatives at the Y boundaries // Setting to free boundaries if (phi.hasParallelSlices()) { From 8ac3efb40e7a5ba02b31ff3d18a3b820ba005376 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 30 Aug 2023 12:45:36 -0700 Subject: [PATCH 09/18] Adding TCV-X21 input and analysis scripts Python scripts to gather data, convert to TCV-X21 dataset, and make plots for comparison to experiment. --- examples/tcv-x21/convert_to_tcvx21.py | 198 ++++++++++++++++++++++++ examples/tcv-x21/data/BOUT.inp | 138 +++++++++++++++++ examples/tcv-x21/gather_data.py | 215 ++++++++++++++++++++++++++ examples/tcv-x21/make_tcvx21_plots.py | 101 ++++++++++++ 4 files changed, 652 insertions(+) create mode 100644 examples/tcv-x21/convert_to_tcvx21.py create mode 100644 examples/tcv-x21/data/BOUT.inp create mode 100644 examples/tcv-x21/gather_data.py create mode 100644 examples/tcv-x21/make_tcvx21_plots.py diff --git a/examples/tcv-x21/convert_to_tcvx21.py b/examples/tcv-x21/convert_to_tcvx21.py new file mode 100644 index 00000000..5f847c65 --- /dev/null +++ b/examples/tcv-x21/convert_to_tcvx21.py @@ -0,0 +1,198 @@ +from pathlib import Path +import tcvx21 +from tcvx21.record_c.record_writer_m import RecordWriter + + +def write_x21_dataset(hermes_data: dict, output_file: Path): + """ + Write data to NetCDF, in the format of the TCV-X21 datasets + """ + + result = { + "LFS-LP": { + "name": "Low-field-side target Langmuir probes", + "hermes_location": "lfs", + "observables": { + "density": { + "name": "Plasma density", + "hermes_name": "Ne", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "electron_temp": { + "name": "Electron temperature", + "hermes_name": "Te", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "ion_temp": { + "name": "Ion temperature", + "hermes_name": "Ti", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "potential": { + "name": "Plasma potential", + "hermes_name": "phi", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "current": { + "name": "Parallel current", + "hermes_name": "Jpar", + "experimental_hierarchy": 1, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "vfloat": { + "name": "Floating potential", + "hermes_name": "Vfl", + "experimental_hierarchy": 1, + "dimensionality": 1, + "simulation_hierarchy": 2, + }, + }, + }, + "HFS-LP": { + "name": "High-field-side target Langmuir probes", + "hermes_location": "hfs", + "observables": { + "density": { + "name": "Plasma density", + "hermes_name": "Ne", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "electron_temp": { + "name": "Electron temperature", + "hermes_name": "Te", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "ion_temp": { + "name": "Ion temperature", + "hermes_name": "Ti", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "potential": { + "name": "Plasma potential", + "hermes_name": "phi", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "current": { + "name": "Parallel current", + "hermes_name": "Jpar", + "experimental_hierarchy": 1, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "vfloat": { + "name": "Floating potential", + "hermes_name": "Vfl", + "experimental_hierarchy": 1, + "dimensionality": 1, + "simulation_hierarchy": 2, + }, + }, + }, + "FHRP": { + "name": "Outboard midplane reciprocating probe", + "hermes_location": "omp", + "observables": { + "density": { + "name": "Plasma density", + "hermes_name": "Ne", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "electron_temp": { + "name": "Electron temperature", + "hermes_name": "Te", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "ion_temp": { + "name": "Ion temperature", + "hermes_name": "Ti", + "experimental_hierarchy": -1, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "potential": { + "name": "Plasma potential", + "hermes_name": "phi", + "experimental_hierarchy": 2, + "dimensionality": 1, + "simulation_hierarchy": 1, + }, + "vfloat": { + "name": "Floating potential", + "hermes_name": "Vfl", + "experimental_hierarchy": 1, + "dimensionality": 1, + "simulation_hierarchy": 2, + }, + }, + }, + } + + # Add Hermes-3 data + for dname, diagnostic in result.items(): + observables = diagnostic["observables"] + location = diagnostic["hermes_location"] + for oname, observable in observables.items(): + hermes_obs = hermes_data[observable["hermes_name"]][location] + observable["units"] = hermes_obs["units"] + observable["values"] = hermes_obs["mean"] + observable["errors"] = hermes_obs["std"] + observable["Ru"] = hermes_obs["Ru"] + observable["Ru_units"] = hermes_obs["Ru_units"] + + additional_attributes = {} + + writer = RecordWriter( + file_path=output_file, + descriptor="Hermes-3", + description="Hermes-3 simulation dataset", + allow_overwrite=True, + ) + writer.write_data_dict(result, additional_attributes) + + +if __name__ == "__main__": + import argparse + import pickle + + parser = argparse.ArgumentParser( + description="Convert Hermes-3 pickle file into TCV-X21 NetCDF file" + ) + + parser.add_argument( + "pickle_file_path", type=str, help="Pickle file containing Hermes-3 data" + ) + + parser.add_argument( + "-o", + "--output", + default="hermes_tcvx21_data.nc", + help="Output NetCDF file in TCV-X21 format", + ) + + args = parser.parse_args() + + with open(args.pickle_file_path, "rb") as f: + data = pickle.load(f) + + write_x21_dataset(data, Path(args.output)) diff --git a/examples/tcv-x21/data/BOUT.inp b/examples/tcv-x21/data/BOUT.inp new file mode 100644 index 00000000..ce026a50 --- /dev/null +++ b/examples/tcv-x21/data/BOUT.inp @@ -0,0 +1,138 @@ +# Turbulence simulation + +nout = 200 +timestep = 100 + +MZ = 81 + +zperiod = 5 + +[mesh] + +file = "65402_68x32_revIp_wide_fixBp_curv.nc" + +extrapolate_y = false # Can result in negative Jacobians in guard cells + +[mesh:paralleltransform] +type = shifted + +[solver] + +use_precon = true +mxstep = 100000 +cvode_max_order = 3 + +[hermes] +components = (e, i, sound_speed, vorticity, + sheath_boundary, collisions, + diamagnetic_drift, classical_diffusion, ion_viscosity, + polarisation_drift + ) + +Nnorm = 1e19 # Reference density [m^-3] +Bnorm = 1 # Reference magnetic field [T] +Tnorm = 50 # Reference temperature [eV] + +loadmetric = false # Recalculate metric tensor? (false -> use grid values) +normalise_metric = true + +[polarisation_drift] +advection = false +diamagnetic_polarisation = true +average_atomic_mass = vorticity:average_atomic_mass + +[ion_viscosity] +perpendicular = true +diagnose = true + +[diamagnetic_drift] +diamag_form = x * (1 - x) # 0 = gradient; 1 = divergence + +[vorticity] + +exb_advection_simplified = false +diamagnetic = true # Include diamagnetic current? +diamagnetic_polarisation = true # Include diamagnetic drift in polarisation current? +average_atomic_mass = i:AA # Weighted average atomic mass, for polarisaion current +poloidal_flows = true # Include poloidal ExB flow +split_n0 = false # Split phi into n=0 and n!=0 components + +vort_dissipation = false +phi_dissipation = true +phi_sheath_dissipation = true +damp_core_vorticity = true + +phi_boundary_relax = true +phi_boundary_timescale = 1e-6 + +################################################################ +# Electrons + +[e] +# Evolve the electron density, parallel momentum, and fix Te +type = evolve_density, evolve_pressure, evolve_momentum + +AA = 1 / 1836 +charge = -1 + +poloidal_flows = true + +diagnose = true + +low_n_diffuse_perp = true + +[Ne] +neumann_boundary_average_z = true # Neumann average Z boundaries in X + +function = 3.0 - 2.7*x + 1e-3*mixmode(x - z) + +flux = 3e21 # /s + +# sum( Se_src[x,y] * J*dx*dy*2pi ) +# Note: Depends on source shape and mesh +shape_factor = 14.325989540821935 + +source = flux * shape_factor * exp(-((x - 0.05)/0.05)^2) +source_only_in_core = true + +[Pe] +#bndry_core = dirichlet(3.0) +#bndry_all = neumann +neumann_boundary_average_z = true # Neumnn average Z boundaries in X + +function = 3*(1.0 - 0.9*x)^2 + +heating = 60e3 # Power input per species [W] + +T_source = heating / (Ne:flux * 1.602e-19 * 1.5) + +source = Ne:source * T_source * 1.602e-19 +source_only_in_core = true + +# Post-process source: +# float((bd['Pe_src'][-1,:,:,0] * bd['J'] * bd['dx'] * bd['dy'] * 2 * 3.14159).sum()) + +################################################################ +# Ions +[i] +# Set ion density from quasineutrality, evolve parallel flow +type = quasineutral, evolve_pressure, evolve_momentum + +AA = 2 # Atomic mass +charge = 1 + +poloidal_flows = true + +low_n_diffuse_perp = true +#hyper_z_T = 1e-5 + +[Pi] +#bndry_core = dirichlet(3.0) +#bndry_all = neumann +neumann_boundary_average_z = true # Neumann average boundaries in X + +function = 3*(1.0 - 0.9*x)^2 + +source = Pe:source +source_only_in_core = true + diff --git a/examples/tcv-x21/gather_data.py b/examples/tcv-x21/gather_data.py new file mode 100644 index 00000000..720e0f0e --- /dev/null +++ b/examples/tcv-x21/gather_data.py @@ -0,0 +1,215 @@ +import numpy as np +from boutdata import collect +from boututils.datafile import DataFile +import pickle + +qe = 1.602e-19 +AA_me = 1.0 / 1836 # Electron mass divided by proton mass + + +def extract_data(path, gridfilepath, ymid=18): + """ + Read data from a Hermes-3/BOUT++ output directory. + Returned as a dictionary. + + Note: Here we should use xHermes, but the grid file does not + include the Y guard cell values. xHermes/xBOUT therefore fails + to load the dataset due to array size mismatches if Y boundaries + are requested. + + Since extrapolate_y = false in this case, the zShift angle is + constant into the boundary. Hence no mapping to field-aligned + coordinates is needed to reconstruct sheath entrance values. + + NOTE: The code below is specific to this case! + """ + + # Many diagnostics are mapped in flux space to R at midplane + with DataFile(gridfilepath) as grid: + psixy = grid["psixy"] + Rxy = grid["Rxy"] + ixsep = grid["ixseps1"] + + psi_mid = psixy[:, ymid] + R_mid = Rxy[:, ymid] + Rsep = 0.5 * (R_mid[ixsep - 1] + R_mid[ixsep]) + + from scipy.interpolate import interp1d + + R_u = interp1d( + psi_mid, R_mid - Rsep, fill_value="extrapolate" + ) # Linearly interpolate R as a function of psi + + Nnorm = collect("Nnorm", path=path) + Tnorm = collect("Tnorm", path=path) + wci = collect("Omega_ci", path=path) + Cs0 = collect("Cs0", path=path) + time = collect("t_array", path=path) / wci # Seconds + run_id = collect("run_id", path=path) + AA = 2 # Atomic species + + # Dictionary to be populated by the add_var function + result = {} + + def add_var(name, units, data_txyz): + assert len(data_txyz.shape) == 4 + assert data_txyz.shape[1] == psixy.shape[0] # X axis + assert data_txyz.shape[2] == psixy.shape[1] + 4 # Y axis. + # Data includes 4 Y guards that are not in the grid file + + result[name] = {} + + def add_location(location, name, Rx, data_txz): + assert len(Rx.shape) == 1 + assert len(data_txz.shape) == 3 + assert len(Rx) == data_txz.shape[1] + + result[name][location] = { + "name": name, + "units": units, + "nt": data_txz.shape[0], + "tmin": time[0], + "tmax": time[-1], + "duration": time[-1] - time[0], + "nz": data_txz.shape[-1], + "run_ids": [run_id], + "paths": [path], + "gridfilepath": gridfilepath, + "mean": np.mean(data_txz, axis=(0, -1)), # Average in time and Z + "meansq": np.mean(data_txz**2, axis=(0, -1)), + "std": np.std(data_txz, axis=(0, -1)), + "Ru": Rx * 100.0, # Major radius in cm + "Ru_units": "cm", + } + + # Outboard midplane. data includes 2 Y guard cells + add_location("omp", name, R_mid - Rsep, data_txyz[:, :, ymid + 2, :]) + + # Interpolate at high field side, mapping to outboard midplane R + add_location( + "hfs", + name, + R_u(psixy[:, 0]), + 0.5 * (data_txyz[:, :, 1, :] + data_txyz[:, :, 2, :]), + ) + + # Interpolate at low field side + add_location( + "lfs", + name, + R_u(psixy[:, -1]), + 0.5 * (data_txyz[:, :, -2, :] + data_txyz[:, :, -3, :]), + ) + + Ne = collect("Ne", path=path, yguards=True) + add_var("Ne", "1/m^3", Ne * Nnorm) + + Ne_floor = np.clip(Ne, 1e-5, None) + + Pe = collect("Pe", path=path, yguards=True) + add_var("Pe", "Pa", Pe * Nnorm * Tnorm * qe) + + Te = Pe / Ne_floor + add_var("Te", "eV", Te * Tnorm) + + Pi = collect("Pi", path=path, yguards=True) + add_var("Pi", "Pa", Pi * Nnorm * Tnorm * qe) + + Ti = Pi / Ne_floor + add_var("Ti", "eV", Ti * Tnorm) + + phi = collect("phi", path=path, yguards=True) + add_var("phi", "V", phi * Tnorm) + + Vfl = phi - 3 * Te + add_var("Vfl", "V", Vfl * Tnorm) + + # Note: NVi and NVe contain mass factors + NVi = collect("NVi", path=path, yguards=True) + add_var("Vi", "m/s", NVi / (AA * Ne_floor) * Nnorm * Cs0) + + NVe = collect("NVe", path=path, yguards=True) + Jpar = NVi / AA - NVe / AA_me + add_var("Jpar", "A/m^2", Jpar * qe * Nnorm * Cs0) + + return result + + +def combine_data(dataset1, dataset2): + """ + Combine two datasets together. Returns a new dataset. + """ + result = {} + for name in dataset1.keys(): + result[name] = {} + for location in dataset1[name].keys(): + data1 = dataset1[name][location] + data2 = dataset2[name][location] + + # Check that the datasets are compatible + assert data1["name"] == data2["name"] + assert data1["units"] == data2["units"] + assert data1["nz"] == data2["nz"] + assert data1["gridfilepath"] == data2["gridfilepath"] + assert data1["Ru_units"] == data2["Ru_units"] + + # Weighting factors for combining means + nt_tot = data1["nt"] + data2["nt"] + w1 = data1["nt"] / nt_tot + w2 = data2["nt"] / nt_tot + mean = data1["mean"] * w1 + data2["mean"] * w2 + meansq = data1["meansq"] * w1 + data2["meansq"] * w2 + # Re-calculate standard deviation + std = np.sqrt(meansq - mean**2) + + result[name][location] = { + "name": data1["name"], + "units": data1["units"], + "nt": nt_tot, + "tmin": min([data1["tmin"], data2["tmin"]]), + "tmax": max([data1["tmax"], data2["tmax"]]), + "duration": data1["duration"] + data2["duration"], + "nz": data1["nz"], + "run_ids": data1["run_ids"] + data2["run_ids"], + "paths": data1["paths"] + data2["paths"], + "gridfilepath": data1["gridfilepath"], + "mean": mean, + "meansq": meansq, + "std": std, + "Ru": data1["Ru"], + "Ru_units": data1["Ru_units"], + } + return result + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser( + description="Gather Hermes-3 data from one or more directories" + ) + + parser.add_argument( + "paths", type=str, nargs="+", help="Paths containing datasets to concatenate" + ) + + parser.add_argument( + "-o", "--output", default="data.pickle", help="Pickle file to write to" + ) + + parser.add_argument( + "-g", "--grid", default="bout.grd.nc", type=str, help="Grid file" + ) + + args = parser.parse_args() + + print(f"Got {len(args.paths)} files: {args.paths}") + print(f"Outputting to '{args.output}'") + + data = extract_data(args.paths[0], args.grid, ymid=18) + for path in args.paths[1:]: + data2 = extract_data(path, args.grid, ymid=18) + data = combine_data(data, data2) + + with open(args.output, "wb") as f: + pickle.dump(data, f) diff --git a/examples/tcv-x21/make_tcvx21_plots.py b/examples/tcv-x21/make_tcvx21_plots.py new file mode 100644 index 00000000..68006d9e --- /dev/null +++ b/examples/tcv-x21/make_tcvx21_plots.py @@ -0,0 +1,101 @@ +import tcvx21 +from tcvx21 import Quantity +import matplotlib.pyplot as plt +from tcvx21.record_c import Record +import numpy as np +from scipy.stats import linregress +from scipy.optimize import curve_fit + +plt.style.use(tcvx21.style_sheet) + +tcv = Record("TCV-X21/1.experimental_data/TCV_forward_field.nc", color="C0") + +hermes = Record("hermes_tcvx21_data.nc", color="C1") +hermes_highres = Record("hermes_tcvx21_highres.nc", color="C2") + +# Get Te, ne profiles + +ne_tcv_ts = tcv.get_observable("TS", "density") +ne_tcv_fhrp = tcv.get_observable("FHRP", "density") +ne_lowres = hermes.get_observable("FHRP", "density") +ne_highres = hermes_highres.get_observable("FHRP", "density") + + +def make_mask(positions): + return np.logical_and( + positions > Quantity(0.0, "cm"), + positions < Quantity(1.5, "cm"), + ) + + +def length_scale(dataset): + mask = make_mask(dataset.positions) + x = dataset.positions[mask].magnitude + y = dataset.values[mask].magnitude + y_err = dataset.errors[mask].magnitude + + # Get initial guess + result = linregress(x, np.log(y)) + L = -1.0 / result.slope + y0 = np.exp(result.intercept) + + # Use curve_fit to calculate errors + def logfit(x, y0, L): + return y0 * np.exp(-x / L) + + result = curve_fit(logfit, x, y, p0=(y0, L), sigma=y_err, absolute_sigma=True) + + y0 = result[0][0] + L = result[0][1] + y0_err, L_err = np.sqrt(np.diag(result[1])) + + return ( + ( + Quantity(L, dataset.positions.units), + Quantity(L_err, dataset.positions.units), + ), + (Quantity(y0, dataset.values.units), Quantity(y0_err, dataset.values.units)), + ) + + +for label, dataset in [ + ("lambda_n (TS)", ne_tcv_ts), + ("lambda_n (FHRP)", ne_tcv_fhrp), + ("lambda_n (Hermes-3 lowres)", ne_lowres), + ("lambda_n (Hermes-3 highres)", ne_highres), + ("lambda_T (TS)", tcv.get_observable("TS", "electron_temp")), + ("lambda_T (FHRP)", tcv.get_observable("FHRP", "electron_temp")), + ("lambda_T (Hermes-3 lowres)", hermes.get_observable("FHRP", "electron_temp")), + ( + "lambda_T (Hermes-3 highres)", + hermes_highres.get_observable("FHRP", "electron_temp"), + ), +]: + result = length_scale(dataset) + L, L_err = result[0] + y0, y0_err = result[1] + print(f"{label}: {L} +/- {L_err}\n\tIntercept: {y0} +/- {y0_err}") + +# Make plots of midplane and target profiles +for region, measurement in [ + ("LFS-LP", "density"), + ("LFS-LP", "electron_temp"), + ("LFS-LP", "potential"), + ("LFS-LP", "current"), + ("LFS-LP", "vfloat"), + ("FHRP", "density"), + ("FHRP", "electron_temp"), + ("FHRP", "potential"), + ("HFS-LP", "density"), + ("HFS-LP", "electron_temp"), + ("HFS-LP", "potential"), + ("HFS-LP", "current"), + ("HFS-LP", "vfloat"), +]: + fig, ax = plt.subplots() + tcv.get_observable(region, measurement).plot(ax=ax) + hermes.get_observable(region, measurement).plot(ax=ax) + hermes_highres.get_observable(region, measurement).plot(ax=ax) + ax.set_title(f"{region}:{measurement}") + fig.legend() + plt.show() From a942af7828f591caecba6b6d49547ffb9747888a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 31 Aug 2023 08:20:56 -0700 Subject: [PATCH 10/18] Capture and save slope limiter name HERMES_SLOPE_LIMITER is captured as a const string `hermes::limiter_typename`. That is then saved to the options and output dump files. --- hermes-3.cxx | 6 ++++++ include/hermes_build_config.hxx.in | 1 + 2 files changed, 7 insertions(+) diff --git a/hermes-3.cxx b/hermes-3.cxx index 4bf08419..57e1a74c 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -162,6 +162,10 @@ int Hermes::init(bool restarting) { options["revision"] = hermes::version::revision; options["revision"].setConditionallyUsed(); + output.write("Slope limiter: {}\n", hermes::limiter_typename); + options["slope_limiter"] = hermes::limiter_typename; + options["slope_limiter"].setConditionallyUsed(); + // Choose normalisations Tnorm = options["Tnorm"].doc("Reference temperature [eV]").withDefault(100.); Nnorm = options["Nnorm"].doc("Reference density [m^-3]").withDefault(1e19); @@ -277,6 +281,8 @@ void Hermes::outputVars(Options& options) { // Save the Hermes version in the output dump files options["HERMES_REVISION"].force(hermes::version::revision); + options["HERMES_SLOPE_LIMITER"].force(hermes::limiter_typename); + // Save normalisation quantities. These may be used by components // to calculate conversion factors to SI units diff --git a/include/hermes_build_config.hxx.in b/include/hermes_build_config.hxx.in index e529c58b..e084258e 100644 --- a/include/hermes_build_config.hxx.in +++ b/include/hermes_build_config.hxx.in @@ -9,6 +9,7 @@ namespace hermes { /// Slope limiter to use in advection operators using Limiter=FV::@HERMES_SLOPE_LIMITER@; + const char* const limiter_typename = "@HERMES_SLOPE_LIMITER@"; } #endif // HERMES_BUILD_CONFIG_HXX From aecbafe5278ead6d00b427eb05900d865d0fa685 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 3 Sep 2023 14:17:40 +0100 Subject: [PATCH 11/18] Add neutral particle/energy sink for pump in wall - Using v_th as velocity for flux like neutral_boundary - Calculating free BC for neutrals but not propagating to state --- src/recycling.cxx | 54 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 51 insertions(+), 3 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index 185e79e5..90516bba 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -2,6 +2,7 @@ #include "../include/recycling.hxx" #include // for trim, strsplit +#include "../include/hermes_utils.hxx" // For indexAt #include #include #include @@ -131,6 +132,11 @@ void Recycling::transform(Options& state) { Options& species_to = state["species"][channel.to]; + const Field3D Nn = get(species_to["density"]); + const Field3D Pn = get(species_to["pressure"]); + const Field3D Tn = get(species_to["temperature"]); + const BoutReal AAn = get(species_to["AA"]); + // Recycling particle and energy sources will be added to these global sources // which are then passed to the density and pressure equations density_source = species_to.isSet("density_source") @@ -252,10 +258,52 @@ void Recycling::transform(Options& state) { // Divide by volume to get source BoutReal recycle_source = recycle_particle_flow / volume; - // Add to appropriate diagnostic field depending if pump or not + // Compute the neutral loss sink and combine with the recycled source to compute overall neutral sources if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { - pump_recycle_density_source(mesh->xend, iy, iz) += recycle_source; - pump_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy; + + auto i = indexAt(Nn, mesh->lastX(), iy, iz); // Final domain cell + auto is = i.xm(); // Second to final domain cell + auto ig = i.xp(); // First guard cell + + // Free boundary condition on Nn, Pn, Tn + // These are NOT communicated back into state and will exist only in this component + // This will prevent neutrals leaking through cross-field transport from neutral_mixed or other components + // While enabling us to still calculate radial wall fluxes separately here + const BoutReal Nn_guard = SQ(Nn[i]) / Nn[is]; + const BoutReal Pn_guard = SQ(Pn[i]) / Pn[is]; + const BoutReal Tn_guard = SQ(Tn[i]) / Tn[is]; + + // Calculate wall conditions + const BoutReal nnsheath = 0.5 * (Nn[i] + Nn_guard); + const BoutReal tnsheath = 0.5 * (Tn[i] + Tn_guard); + const BoutReal v_th = sqrt(tnsheath / AAn); + + // Convert dy to poloidal length: dl = dy * sqrt(g22) = dy * h_theta + // Convert dz to toroidal length: = dz*sqrt(g_33) = dz * R = 2piR + // Calculate radial wall area in [m^2] + // Calculate final cell volume [m^3] + BoutReal dpolsheath = 0.5*(coord->dy[i] + coord->dy[ig]) * 1/( 0.5*(sqrt(coord->g22[i]) + sqrt(coord->g22[ig])) ); + BoutReal dtorsheath = 0.5*(coord->dz[i] + coord->dz[ig]) * 0.5*(sqrt(coord->g_33[i]) + sqrt(coord->g_33[ig])); + BoutReal dasheath = dpolsheath * dtorsheath; // [m^2] + BoutReal dv = coord->J[i] * coord->dx[i] * coord->dy[i] * coord->dz[i]; + + // Calculate particle and energy fluxes of neutrals hitting the pump + // Assume thermal velocity greater than perpendicular velocity and use it for flux calc + BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] + BoutReal psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + + // Use gamma=3.5 as per Stangeby p.69, total energy of drifting Maxwellian + BoutReal eflow = 3.5 * tnsheath * v_th * nnsheath * dasheath; // [W] + BoutReal esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] + + // Pump puts neutral particle and energy source in final domain cell + // Source accounts for recycled ions and the particle sink due to neutrals hitting the pump + // Note that the ions are explicitly transported out of the domain by anomalous_diffusion.cxx + // the pump picks this up and adds a recycling source based on this, but doesn't need an ion sink. + // Neutrals are not explicitly transported out by any component and must be taken out by the sink below. + // Pump multiplier controls both fraction of recycled ions and fraction of returned neutrals + pump_recycle_density_source(mesh->xend, iy, iz) += recycle_source - psink; + pump_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy - esink; } else { wall_recycle_density_source(mesh->xend, iy, iz) += recycle_source; wall_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy; From d5d5cdfba6c3625f8cc8997279c2600241697ebb Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 3 Sep 2023 14:22:24 +0100 Subject: [PATCH 12/18] Add diagnostic on pump location --- src/recycling.cxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index 90516bba..4781ec06 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -24,7 +24,7 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { // Neutral pump - // Mark cells as having a pump by setting the Field3D is_pump to 1 in the grid file + // Mark cells as having a pump by setting the Field2D is_pump to 1 in the grid file // Works only on SOL and PFR edges, where it locally modifies the recycle multiplier to the pump albedo is_pump = 0.0; mesh->get(is_pump, std::string("is_pump")); @@ -449,6 +449,11 @@ void Recycling::outputVars(Options& state) { {"standard_name", "energy source"}, {"long_name", std::string("Pump recycling energy source of ") + channel.to}, {"source", "recycling"}}); + + set_with_attrs(state[{std::string("is_pump")}], is_pump, + {{"standard_name", "neutral pump location"}, + {"long_name", std::string("Neutral pump location")}, + {"source", "recycling"}}); } } From ef81db6bbcb705ce94a7a57e24f2d1d1479d6075 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 3 Sep 2023 18:09:41 +0100 Subject: [PATCH 13/18] Add pump neutral sink at PFR --- src/recycling.cxx | 61 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 51 insertions(+), 10 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index 4781ec06..a253c5bd 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -243,7 +243,7 @@ void Recycling::transform(Options& state) { // If cell is a pump, overwrite multiplier with pump multiplier BoutReal multiplier = channel.pfr_multiplier; - if ((is_pump(mesh->xend, iy, iz) == 1.0) and (neutral_pump)) { + if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { multiplier = channel.pump_multiplier; } @@ -261,7 +261,7 @@ void Recycling::transform(Options& state) { // Compute the neutral loss sink and combine with the recycled source to compute overall neutral sources if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { - auto i = indexAt(Nn, mesh->lastX(), iy, iz); // Final domain cell + auto i = indexAt(Nn, mesh->xend, iy, iz); // Final domain cell auto is = i.xm(); // Second to final domain cell auto ig = i.xp(); // First guard cell @@ -269,14 +269,14 @@ void Recycling::transform(Options& state) { // These are NOT communicated back into state and will exist only in this component // This will prevent neutrals leaking through cross-field transport from neutral_mixed or other components // While enabling us to still calculate radial wall fluxes separately here - const BoutReal Nn_guard = SQ(Nn[i]) / Nn[is]; - const BoutReal Pn_guard = SQ(Pn[i]) / Pn[is]; - const BoutReal Tn_guard = SQ(Tn[i]) / Tn[is]; + BoutReal nnguard = SQ(Nn[i]) / Nn[is]; + BoutReal pnguard = SQ(Pn[i]) / Pn[is]; + BoutReal tnguard = SQ(Tn[i]) / Tn[is]; // Calculate wall conditions - const BoutReal nnsheath = 0.5 * (Nn[i] + Nn_guard); - const BoutReal tnsheath = 0.5 * (Tn[i] + Tn_guard); - const BoutReal v_th = sqrt(tnsheath / AAn); + BoutReal nnsheath = 0.5 * (Nn[i] + nnguard); + BoutReal tnsheath = 0.5 * (Tn[i] + tnguard); + BoutReal v_th = sqrt(tnsheath / AAn); // Convert dy to poloidal length: dl = dy * sqrt(g22) = dy * h_theta // Convert dz to toroidal length: = dz*sqrt(g_33) = dz * R = 2piR @@ -354,8 +354,49 @@ void Recycling::transform(Options& state) { // Add to appropriate diagnostic field depending if pump or not if ((is_pump(mesh->xstart, iy) == 1.0) and (neutral_pump)) { - pump_recycle_density_source(mesh->xstart, iy, iz) += recycle_source; - pump_recycle_energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; + auto i = indexAt(Nn, mesh->xstart, iy, iz); // Final domain cell + auto is = i.xp(); // Second to final domain cell + auto ig = i.xm(); // First guard cell + + // Free boundary condition on Nn, Pn, Tn + // These are NOT communicated back into state and will exist only in this component + // This will prevent neutrals leaking through cross-field transport from neutral_mixed or other components + // While enabling us to still calculate radial wall fluxes separately here + BoutReal nnguard = SQ(Nn[i]) / Nn[is]; + BoutReal pnguard = SQ(Pn[i]) / Pn[is]; + BoutReal tnguard = SQ(Tn[i]) / Tn[is]; + + // Calculate wall conditions + BoutReal nnsheath = 0.5 * (Nn[i] + nnguard); + BoutReal tnsheath = 0.5 * (Tn[i] + tnguard); + BoutReal v_th = sqrt(tnsheath / AAn); + + // Convert dy to poloidal length: dl = dy * sqrt(g22) = dy * h_theta + // Convert dz to toroidal length: = dz*sqrt(g_33) = dz * R = 2piR + // Calculate radial wall area in [m^2] + // Calculate final cell volume [m^3] + BoutReal dpolsheath = 0.5*(coord->dy[i] + coord->dy[ig]) * 1/( 0.5*(sqrt(coord->g22[i]) + sqrt(coord->g22[ig])) ); + BoutReal dtorsheath = 0.5*(coord->dz[i] + coord->dz[ig]) * 0.5*(sqrt(coord->g_33[i]) + sqrt(coord->g_33[ig])); + BoutReal dasheath = dpolsheath * dtorsheath; // [m^2] + BoutReal dv = coord->J[i] * coord->dx[i] * coord->dy[i] * coord->dz[i]; + + // Calculate particle and energy fluxes of neutrals hitting the pump + // Assume thermal velocity greater than perpendicular velocity and use it for flux calc + BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] + BoutReal psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + + // Use gamma=3.5 as per Stangeby p.69, total energy of drifting Maxwellian + BoutReal eflow = 3.5 * tnsheath * v_th * nnsheath * dasheath; // [W] + BoutReal esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] + + // Pump puts neutral particle and energy source in final domain cell + // Source accounts for recycled ions and the particle sink due to neutrals hitting the pump + // Note that the ions are explicitly transported out of the domain by anomalous_diffusion.cxx + // the pump picks this up and adds a recycling source based on this, but doesn't need an ion sink. + // Neutrals are not explicitly transported out by any component and must be taken out by the sink below. + // Pump multiplier controls both fraction of recycled ions and fraction of returned neutrals + pump_recycle_density_source(mesh->xend, iy, iz) += recycle_source - psink; + pump_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy - esink; } else { wall_recycle_density_source(mesh->xstart, iy, iz) += recycle_source; wall_recycle_energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; From 7a3ee544966ad7199098f066118425b61f9cfe0b Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 3 Sep 2023 18:10:45 +0100 Subject: [PATCH 14/18] Change pump diagnostic name to "Sd_pump" - This is for clarity - now it's not just the neutral recycling flow but also the neutral sink. --- src/recycling.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index a253c5bd..fb00f040 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -475,7 +475,7 @@ void Recycling::outputVars(Options& state) { // Neutral pump if (neutral_pump) { - set_with_attrs(state[{std::string("S") + channel.to + std::string("_pump_recycle")}], pump_recycle_density_source, + set_with_attrs(state[{std::string("S") + channel.to + std::string("_pump")}], pump_recycle_density_source, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, {"conversion", Nnorm * Omega_ci}, @@ -483,7 +483,7 @@ void Recycling::outputVars(Options& state) { {"long_name", std::string("Pump recycling particle source of ") + channel.to}, {"source", "recycling"}}); - set_with_attrs(state[{std::string("E") + channel.to + std::string("_pump_recycle")}], pump_recycle_energy_source, + set_with_attrs(state[{std::string("E") + channel.to + std::string("_pump")}], pump_recycle_energy_source, {{"time_dimension", "t"}, {"units", "W m^-3"}, {"conversion", Pnorm * Omega_ci}, From e62296cd897448c89bb0eb21bc629b1495291c44 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 4 Sep 2023 16:29:18 +0100 Subject: [PATCH 15/18] Fix neutral pump sinks not added to state --- src/recycling.cxx | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index fb00f040..92c1d6f7 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -258,6 +258,10 @@ void Recycling::transform(Options& state) { // Divide by volume to get source BoutReal recycle_source = recycle_particle_flow / volume; + // Pump source terms + BoutReal esink = 0; + BoutReal psink = 0; + // Compute the neutral loss sink and combine with the recycled source to compute overall neutral sources if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { @@ -290,11 +294,11 @@ void Recycling::transform(Options& state) { // Calculate particle and energy fluxes of neutrals hitting the pump // Assume thermal velocity greater than perpendicular velocity and use it for flux calc BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] - BoutReal psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] // Use gamma=3.5 as per Stangeby p.69, total energy of drifting Maxwellian BoutReal eflow = 3.5 * tnsheath * v_th * nnsheath * dasheath; // [W] - BoutReal esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] + esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] // Pump puts neutral particle and energy source in final domain cell // Source accounts for recycled ions and the particle sink due to neutrals hitting the pump @@ -304,6 +308,7 @@ void Recycling::transform(Options& state) { // Pump multiplier controls both fraction of recycled ions and fraction of returned neutrals pump_recycle_density_source(mesh->xend, iy, iz) += recycle_source - psink; pump_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy - esink; + } else { wall_recycle_density_source(mesh->xend, iy, iz) += recycle_source; wall_recycle_energy_source(mesh->xend, iy, iz) += recycle_source * channel.sol_energy; @@ -311,8 +316,9 @@ void Recycling::transform(Options& state) { // Add to density source which will be picked up by evolve_density.cxx // Add to energy source which will be picked up by evolve_pressure.cxx - density_source(mesh->xend, iy, iz) += recycle_source; - energy_source(mesh->xend, iy, iz) += recycle_source * channel.pfr_energy; + // psink and esink are additional sinks from the neutral pump which are 0 if disabled + density_source(mesh->xend, iy, iz) += recycle_source - psink; + energy_source(mesh->xend, iy, iz) += recycle_source * channel.pfr_energy - esink; } } @@ -352,6 +358,10 @@ void Recycling::transform(Options& state) { // Divide by volume to get source BoutReal recycle_source = recycle_particle_flow / volume; + // Pump source terms + BoutReal esink = 0; + BoutReal psink = 0; + // Add to appropriate diagnostic field depending if pump or not if ((is_pump(mesh->xstart, iy) == 1.0) and (neutral_pump)) { auto i = indexAt(Nn, mesh->xstart, iy, iz); // Final domain cell @@ -383,11 +393,11 @@ void Recycling::transform(Options& state) { // Calculate particle and energy fluxes of neutrals hitting the pump // Assume thermal velocity greater than perpendicular velocity and use it for flux calc BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] - BoutReal psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] // Use gamma=3.5 as per Stangeby p.69, total energy of drifting Maxwellian BoutReal eflow = 3.5 * tnsheath * v_th * nnsheath * dasheath; // [W] - BoutReal esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] + esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] // Pump puts neutral particle and energy source in final domain cell // Source accounts for recycled ions and the particle sink due to neutrals hitting the pump @@ -404,8 +414,9 @@ void Recycling::transform(Options& state) { // Add to density source which will be picked up by evolve_density.cxx // Add to energy source which will be picked up by evolve_pressure.cxx - density_source(mesh->xstart, iy, iz) += recycle_source; - energy_source(mesh->xstart, iy, iz) += recycle_source * channel.pfr_energy; + // psink and esink are additional sinks from the neutral pump which are 0 if disabled + density_source(mesh->xend, iy, iz) += recycle_source - psink; + energy_source(mesh->xend, iy, iz) += recycle_source * channel.pfr_energy - esink; } } From 34661f2e370b3a25b988c587a4e32a5ef9c318ac Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 5 Sep 2023 17:19:49 +0100 Subject: [PATCH 16/18] Fix v_th, pump hflux gamma in neutral_bndry/pump The thermal velocity was naively calculated and was wrong. Replaced with thermal velocity in one direction for a static Maxwellian from Stangeby. Pump heat flux gamma changed to 2 to represent a static, not drifting Maxwellian. --- src/neutral_boundary.cxx | 12 ++++++------ src/recycling.cxx | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/neutral_boundary.cxx b/src/neutral_boundary.cxx index 10190a55..2a2dca30 100644 --- a/src/neutral_boundary.cxx +++ b/src/neutral_boundary.cxx @@ -99,7 +99,7 @@ void NeutralBoundary::transform(Options& state) { const BoutReal tnsheath = 0.5 * (Tn[im] + Tn[i]); // Thermal speed - const BoutReal v_th = sqrt(tnsheath / AA); + const BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AA) ); // Stangeby p.69 eqns. 2.21, 2.24 // Calculate effective gamma from particle and energy reflection coefficients BoutReal target_gamma_heat = 1 - target_energy_refl_factor * target_fast_refl_fraction @@ -142,7 +142,7 @@ void NeutralBoundary::transform(Options& state) { const BoutReal tnsheath = 0.5 * (Tn[ip] + Tn[i]); // Thermal speed - const BoutReal v_th = sqrt(tnsheath / AA); + const BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AA) ); // Stangeby p.69 eqns. 2.21, 2.24 // Calculate effective gamma from particle and energy reflection coefficients BoutReal target_gamma_heat = 1 - target_energy_refl_factor * target_fast_refl_fraction @@ -177,8 +177,8 @@ void NeutralBoundary::transform(Options& state) { const BoutReal nnsheath = 0.5 * (Nn[ig] + Nn[i]); const BoutReal tnsheath = 0.5 * (Tn[ig] + Tn[i]); - // Thermal speed in one direction only (?) - const BoutReal v_th = sqrt(tnsheath / AA); + // Thermal speed of static Maxwellian in one direction + const BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AA) ); // Stangeby p.69 eqns. 2.21, 2.24 // Calculate effective gamma from particle and energy reflection coefficients BoutReal sol_gamma_heat = 1 - sol_energy_refl_factor * sol_fast_refl_fraction @@ -224,8 +224,8 @@ void NeutralBoundary::transform(Options& state) { const BoutReal nnsheath = 0.5 * (Nn[ig] + Nn[i]); const BoutReal tnsheath = 0.5 * (Tn[ig] + Tn[i]); - // Thermal speed in one direction only (?) - const BoutReal v_th = sqrt(tnsheath / AA); + // Thermal speed of static Maxwellian in one direction + const BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AA) ); // Stangeby p.69 eqns. 2.21, 2.24 // Calculate effective gamma from particle and energy reflection coefficients diff --git a/src/recycling.cxx b/src/recycling.cxx index 92c1d6f7..e81bf377 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -280,7 +280,7 @@ void Recycling::transform(Options& state) { // Calculate wall conditions BoutReal nnsheath = 0.5 * (Nn[i] + nnguard); BoutReal tnsheath = 0.5 * (Tn[i] + tnguard); - BoutReal v_th = sqrt(tnsheath / AAn); + BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AAn) ); // Stangeby p.69 eqns. 2.21, 2.24 // Convert dy to poloidal length: dl = dy * sqrt(g22) = dy * h_theta // Convert dz to toroidal length: = dz*sqrt(g_33) = dz * R = 2piR @@ -296,8 +296,8 @@ void Recycling::transform(Options& state) { BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] - // Use gamma=3.5 as per Stangeby p.69, total energy of drifting Maxwellian - BoutReal eflow = 3.5 * tnsheath * v_th * nnsheath * dasheath; // [W] + // Use gamma=2.0 as per Stangeby p.69, total energy of static Maxwellian + BoutReal eflow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] // Pump puts neutral particle and energy source in final domain cell @@ -379,7 +379,7 @@ void Recycling::transform(Options& state) { // Calculate wall conditions BoutReal nnsheath = 0.5 * (Nn[i] + nnguard); BoutReal tnsheath = 0.5 * (Tn[i] + tnguard); - BoutReal v_th = sqrt(tnsheath / AAn); + BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AAn) ); // Stangeby p.69 eqns. 2.21, 2.24 // Convert dy to poloidal length: dl = dy * sqrt(g22) = dy * h_theta // Convert dz to toroidal length: = dz*sqrt(g_33) = dz * R = 2piR @@ -395,8 +395,8 @@ void Recycling::transform(Options& state) { BoutReal pflow = v_th * nnsheath * dasheath; // [s^-1] psink = pflow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] - // Use gamma=3.5 as per Stangeby p.69, total energy of drifting Maxwellian - BoutReal eflow = 3.5 * tnsheath * v_th * nnsheath * dasheath; // [W] + // Use gamma=2.0 as per Stangeby p.69, total energy of static Maxwellian + BoutReal eflow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] esink = eflow / dv * (1 - multiplier); // heatsink [W m^-3] // Pump puts neutral particle and energy source in final domain cell From 33ed3cfcf00db780987e2a91e33ac62861db881f Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 25 Sep 2023 15:34:32 +0200 Subject: [PATCH 17/18] Add neutral pump documentation --- docs/sphinx/components.rst | 52 +++++++++++++++++++++++++++++++++++--- 1 file changed, 49 insertions(+), 3 deletions(-) diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index 1a287961..e5a35533 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -1189,9 +1189,55 @@ Each returning atom has an energy of 3.5eV: sol_recycle_multiplier = 1 # Recycling fraction sol_recycle_energy = 3.5 # Energy of recycled particles [eV] - sol_recycle = true - sol_recycle_multiplier = 1 # Recycling fraction - sol_recycle_energy = 3.5 # Energy of recycled particles [eV] + pfr_recycle = true + pfr_recycle_multiplier = 1 # Recycling fraction + pfr_recycle_energy = 3.5 # Energy of recycled particles [eV] + +Neutral pump +^^^^^^^^^^^^^^^ + +The recycling component also features a neutral pump, which is currently implemented for +the SOL and PFR edges only, and so is not available in 1D. The pump is a region of the wall +which facilitates particle loss by incomplete recycling and neutral absorption. + +The particle loss rate :math:`\Gamma_{N_{n}}` is the sum of the incident ions that are not recycled and the +incident neutrals which are not reflected, both of which are controlled by the pump multiplier :math:`M_{p}` +which is set by the `pump_multiplier` option in the input file. The unrecycled ion flux :math:`\Gamma_{N_{i}}^{unrecycled}` is calculated exactly the same +as for edge recycling but with `pump_multiplier` replacing the `recycle_multiplier`. + +.. math:: + + \begin{aligned} + \Gamma_{N_{n}} &= \Gamma_{N_{i}}^{unrecycled} + M_{p} \times \Gamma_{N_{n}}^{incident} \\ + \Gamma_{N_{n}}^{incident} &= N_{n} v_{th} = N_{n} \frac{1}{4} \sqrt{\frac{8 T_{n}}{\pi m_{n}}} \\ + \end{aligned} + +Where the thermal velocity formulation is for a static maxwellian in 1D (see Stangeby p.64, eqns 2.21, 2.24) +and the temperature is in `eV`. + +The heat loss rate :math:`\Gamma_{E_{n}}` is calculated as: + +.. math:: + + \begin{aligned} + \Gamma_{E_{n}} &= \Gamma_{E_{i}}^{unrecycled} + M_{p} \times \Gamma_{E_{n}}^{incident} \\ + \Gamma_{E_{n}}^{incident} &= \gamma T_{n} N_{n} v_{th} = 2 T_{n} N_{n} \frac{1}{4} \sqrt{\frac{8 T_{n}}{\pi m_{n}}} \\ + \end{aligned} + +Where the incident heat flux is for a static maxwellian in 1D (see Stangeby p.69, eqn 2.30). + +The pump will be placed in any cell that + 1. Is the final domain cell before the guard cells + 2. Is on the SOL or PFR edge + 3. Has a `is_pump` value of 1 + +The field `is_pump` must be created by the user and added to the grid file as a `Field2D`. + +Diagnostic variables +^^^^^^^^^^^^^^^ +Diagnostic variables for the recycled particle and energy fluxes are provided separately for the targets, the pump as well as the SOL and PFR which are grouped together as `wall`. +as well as the pump. In addition, the field `is_pump` is saved to help in plotting the pump location. + .. doxygenstruct:: Recycling :members: From 26d1c9c15efd1f64e23b7f986a6eb7fc2b1a9602 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 25 Sep 2023 15:36:28 +0200 Subject: [PATCH 18/18] Add requirement of edge recycling for pump in docs --- docs/sphinx/components.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index e5a35533..0cb370a7 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -1196,9 +1196,11 @@ Each returning atom has an energy of 3.5eV: Neutral pump ^^^^^^^^^^^^^^^ -The recycling component also features a neutral pump, which is currently implemented for +The recycling component also features a neutral pump which is currently implemented for the SOL and PFR edges only, and so is not available in 1D. The pump is a region of the wall -which facilitates particle loss by incomplete recycling and neutral absorption. +which facilitates particle loss by incomplete recycling and neutral absorption. + +The pump requires wall recycling to be enabled on the relevant wall region. The particle loss rate :math:`\Gamma_{N_{n}}` is the sum of the incident ions that are not recycled and the incident neutrals which are not reflected, both of which are controlled by the pump multiplier :math:`M_{p}`