From de9a01364e47c1ba440e57a687a06ace2d905096 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 3 Sep 2023 20:06:29 +0100 Subject: [PATCH 01/19] Add fast recycling at target - Very simple expression at the moment, need to check later whether there are any factors missing from a detailed derivation - Unsure whether neutral energy source is energy or pressure. Previously neutral energy was T_FC * N and added to energy_source, which seems inconsistent. --- include/recycling.hxx | 3 +++ src/recycling.cxx | 56 ++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 56 insertions(+), 3 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index eb30d25e..3892f4ab 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -48,6 +48,8 @@ private: /// 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_energy, sol_energy, pfr_energy; ///< Energy of recycled particle (normalised to Tnorm) + BoutReal target_fast_recycle_fraction, pfr_fast_recycle_fraction, sol_fast_recycle_fraction; ///< Fraction of ions undergoing fast reflection + BoutReal target_fast_recycle_energy_factor, sol_fast_recycle_energy_factor, pfr_fast_recycle_energy_factor; ///< Fraction of energy retained by fast recycled neutrals }; std::vector channels; // Recycling channels @@ -62,6 +64,7 @@ private: Field3D radial_particle_outflow, radial_energy_outflow; ///< Radial fluxes coming from evolve_density and evolve_pressure used in recycling calc + }; namespace { diff --git a/src/recycling.cxx b/src/recycling.cxx index 8906dd3f..fa94339e 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -67,6 +67,36 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { .withDefault(3.0) / Tnorm; // Normalise from eV + BoutReal target_fast_recycle_fraction = + from_options["target_fast_recycle_fraction"] + .doc("Fraction of ions undergoing fast reflection at target") + .withDefault(0); + + BoutReal pfr_fast_recycle_fraction = + from_options["pfr_fast_recycle_fraction"] + .doc("Fraction of ions undergoing fast reflection at pfr") + .withDefault(0); + + BoutReal sol_fast_recycle_fraction = + from_options["sol_fast_recycle_fraction"] + .doc("Fraction of ions undergoing fast reflection at sol") + .withDefault(0); + + BoutReal target_fast_recycle_energy_factor = + from_options["target_fast_recycle_energy_factor"] + .doc("Fraction of energy retained by fast recycled neutrals at target") + .withDefault(0); + + BoutReal sol_fast_recycle_energy_factor = + from_options["sol_fast_recycle_energy_factor"] + .doc("Fraction of energy retained by fast recycled neutrals at sol") + .withDefault(0); + + BoutReal pfr_fast_recycle_energy_factor = + from_options["pfr_fast_recycle_energy_factor"] + .doc("Fraction of energy retained by fast recycled neutrals at pfr") + .withDefault(0); + 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)) { @@ -77,7 +107,9 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { channels.push_back({ from, to, target_recycle_multiplier, sol_recycle_multiplier, pfr_recycle_multiplier, - target_recycle_energy, sol_recycle_energy, pfr_recycle_energy}); + target_recycle_energy, sol_recycle_energy, pfr_recycle_energy, + target_fast_recycle_fraction, pfr_fast_recycle_fraction, sol_fast_recycle_fraction, + target_fast_recycle_energy_factor, sol_fast_recycle_energy_factor, pfr_fast_recycle_energy_factor}); // Boolean flags for enabling recycling in different regions target_recycle = from_options["target_recycle"] @@ -111,6 +143,7 @@ void Recycling::transform(Options& state) { const Field3D N = get(species_from["density"]); const Field3D V = get(species_from["velocity"]); // Parallel flow velocity + const Field3D T = get(species_from["temperature"]); // Ion temperature Options& species_to = state["species"][channel.to]; @@ -123,6 +156,8 @@ void Recycling::transform(Options& state) { ? getNonFinal(species_to["energy_source"]) : 0.0; + + // Recycling at the divertor target plates if (target_recycle) { @@ -156,7 +191,14 @@ void Recycling::transform(Options& state) { / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); // energy of recycled particles - target_recycle_energy_source(r.ind, mesh->ystart, jz) += channel.target_energy * flow + // VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV + // TODO: Is the energy here 3eV * NV or 3/2 * 3eV * NV? It was just 3eV * NV before but it's adding to energy source + BoutReal tisheath = (T(r.ind, mesh->ystart, jz) + T(r.ind, mesh->ystart-1, jz)) * 0.5; // Ion temp at wall + BoutReal neutral_energy = flow * ( + channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction * tisheath // Fast recycling part + + (1 - channel.target_fast_recycle_fraction) * channel.target_energy); // Thermal recycling part + + target_recycle_energy_source(r.ind, mesh->ystart, jz) += neutral_energy / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); energy_source(r.ind, mesh->ystart, jz) += channel.target_energy * flow / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); @@ -190,7 +232,15 @@ void Recycling::transform(Options& state) { density_source(r.ind, mesh->yend, jz) += flow / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - target_recycle_energy_source(r.ind, mesh->yend, jz) += channel.target_energy * flow + // energy of recycled particles + // VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV + // TODO: Is the energy here 3eV * NV or 3/2 * 3eV * NV? It was just 3eV * NV before but it's adding to energy source + BoutReal tisheath = (T(r.ind, mesh->yend, jz) + T(r.ind, mesh->yend+1, jz)) * 0.5; // Ion temp at wall + BoutReal neutral_energy = flow * ( + channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction * tisheath // Fast recycling part + + (1 - channel.target_fast_recycle_fraction) * channel.target_energy); // Thermal recycling part + + target_recycle_energy_source(r.ind, mesh->yend, jz) += neutral_energy / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); energy_source(r.ind, mesh->yend, jz) += channel.target_energy * flow / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); From 799570d0384b606905c1563d8f94ea23a010d277 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 4 Sep 2023 14:24:49 +0100 Subject: [PATCH 02/19] Fix typo in fast recycling The fast recycling energy wasn't being added to the actual energy source, only to the diagnostic energy source --- src/recycling.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index fa94339e..545abbf4 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -200,7 +200,7 @@ void Recycling::transform(Options& state) { target_recycle_energy_source(r.ind, mesh->ystart, jz) += neutral_energy / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - energy_source(r.ind, mesh->ystart, jz) += channel.target_energy * flow + energy_source(r.ind, mesh->ystart, jz) += neutral_energy / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); } } @@ -242,7 +242,7 @@ void Recycling::transform(Options& state) { target_recycle_energy_source(r.ind, mesh->yend, jz) += neutral_energy / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - energy_source(r.ind, mesh->yend, jz) += channel.target_energy * flow + energy_source(r.ind, mesh->yend, jz) += neutral_energy / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); } } From 6de14e39c355ea0e8cc15f45be402bf04762062a Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 27 Sep 2023 15:36:07 +0100 Subject: [PATCH 03/19] Add sheath_bndry_simple heat flux to edge flux tallies - This allows the heat flow through the sheath to be used by other parts of the code --- src/sheath_boundary_simple.cxx | 45 ++++++++++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index 62513b86..14d9d94b 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -279,7 +279,10 @@ void SheathBoundarySimple::transform(Options& state) { } } } - + + // Field to capture total sheath heat flux for diagnostics + Field3D electron_sheath_power_ylow = zeroFrom(Ne); + ////////////////////////////////////////////////////////////////// // Electrons @@ -334,6 +337,14 @@ void SheathBoundarySimple::transform(Options& state) { BoutReal power = flux / (coord->dy[i] * coord->J[i]); electron_energy_source[i] += power; + + // Total heat flux for diagnostic purposes + q = (gamma_e * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; + flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + power = flux / (coord->dy[i] * coord->J[i]); + + electron_sheath_power_ylow[i] += power; // lower Y, so power placed in final domain cell + } } } @@ -386,6 +397,13 @@ void SheathBoundarySimple::transform(Options& state) { BoutReal power = flux / (coord->dy[i] * coord->J[i]); electron_energy_source[i] -= power; + + // Total heat flux for diagnostic purposes + q = (gamma_e * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; + flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + power = flux / (coord->dy[i] * coord->J[i]); + + electron_sheath_power_ylow[ip] -= power; // upper Y, so power placed in first guard cell } } } @@ -399,6 +417,9 @@ void SheathBoundarySimple::transform(Options& state) { // Note: electron_energy_source includes any sources previously set in other components set(electrons["energy_source"], fromFieldAligned(electron_energy_source)); + // Add the total sheath power flux to the tracker of y power flows + add(electrons["energy_flow_ylow"], electron_sheath_power_ylow); + if (IS_SET_NOBOUNDARY(electrons["velocity"])) { setBoundary(electrons["velocity"], fromFieldAligned(Ve)); } @@ -452,6 +473,9 @@ void SheathBoundarySimple::transform(Options& state) { ? toFieldAligned(getNonFinal(species["energy_source"])) : zeroFrom(Ni); + // Field to capture total sheath heat flux for diagnostics + Field3D ion_sheath_power_ylow = zeroFrom(Ne); + if (lower_y) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { @@ -503,6 +527,13 @@ void SheathBoundarySimple::transform(Options& state) { BoutReal power = flux / (coord->dy[i] * coord->J[i]); energy_source[i] += power; + + // Calculation of total heat flux for diagnostic purposes + q = (gamma_i * tisheath + 0.5 * Mi * C_i_sq) * nisheath * visheath; + flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + power = flux / (coord->dy[i] * coord->J[i]); + + ion_sheath_power_ylow[i] += power; // lower Y, so power placed in final domain cell } } } @@ -561,10 +592,17 @@ void SheathBoundarySimple::transform(Options& state) { ASSERT2(std::isfinite(power)); energy_source[i] -= power; // Note: Sign negative because power > 0 + + // Calculation of total heat flux for diagnostic purposes + q = (gamma_i * tisheath + 0.5 * Mi * C_i_sq) * nisheath * visheath; + flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + power = flux / (coord->dy[i] * coord->J[i]); + + ion_sheath_power_ylow[ip] += power; // Upper Y, so power placed in first guard cell } } + } - // Finished boundary conditions for this species // Put the modified fields back into the state. setBoundary(species["density"], fromFieldAligned(Ni)); @@ -582,5 +620,8 @@ void SheathBoundarySimple::transform(Options& state) { // Additional loss of energy through sheath // Note: energy_source already includes previously set values set(species["energy_source"], fromFieldAligned(energy_source)); + + // Add the total sheath power flux to the tracker of y power flows + add(species["energy_flow_ylow"], ion_sheath_power_ylow); } } From e9f224b7978a15af9ec49590aba4ff6bae92ae16 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 27 Sep 2023 18:39:07 +0100 Subject: [PATCH 04/19] Fix flow cross-sectional area and clarify notation --- src/sheath_boundary_simple.cxx | 56 +++++++++++++++++----------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index 14d9d94b..d67395fb 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -330,20 +330,20 @@ void SheathBoundarySimple::transform(Options& state) { * nesheath * vesheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + BoutReal heatflow = q * (coord->J[i] + coord->J[im]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); electron_energy_source[i] += power; // Total heat flux for diagnostic purposes - q = (gamma_e * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; - flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); - power = flux / (coord->dy[i] * coord->J[i]); + q = gamma_e * tesheath * nesheath * vesheath; // Wm-2 + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W - electron_sheath_power_ylow[i] += power; // lower Y, so power placed in final domain cell + electron_sheath_power_ylow[i] += heatflow; // lower Y, so power placed in final domain cell } } @@ -390,20 +390,20 @@ void SheathBoundarySimple::transform(Options& state) { * nesheath * vesheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + BoutReal heatflow = q * (coord->J[i] + coord->J[ip]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); electron_energy_source[i] -= power; // Total heat flux for diagnostic purposes - q = (gamma_e * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; - flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); - power = flux / (coord->dy[i] * coord->J[i]); + q = gamma_e * tesheath * nesheath * vesheath; // Wm-2 + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W - electron_sheath_power_ylow[ip] -= power; // upper Y, so power placed in first guard cell + electron_sheath_power_ylow[ip] -= heatflow; // upper Y, so power placed in first guard cell } } } @@ -520,20 +520,20 @@ void SheathBoundarySimple::transform(Options& state) { * nisheath * visheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); + BoutReal heatflow = q * (coord->J[i] + coord->J[im]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); energy_source[i] += power; // Calculation of total heat flux for diagnostic purposes - q = (gamma_i * tisheath + 0.5 * Mi * C_i_sq) * nisheath * visheath; - flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); - power = flux / (coord->dy[i] * coord->J[i]); + q = gamma_i * tisheath * nisheath * visheath; // Wm-2 + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W - ion_sheath_power_ylow[i] += power; // lower Y, so power placed in final domain cell + ion_sheath_power_ylow[i] += heatflow; // lower Y, so power placed in final domain cell } } } @@ -584,21 +584,21 @@ void SheathBoundarySimple::transform(Options& state) { * nisheath * visheath; // Multiply by cell area to get power - BoutReal flux = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); + BoutReal heatflow = q * (coord->J[i] + coord->J[ip]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); // This omits dx*dz because we divide by dx*dz next // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = flux / (coord->dy[i] * coord->J[i]); + BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); ASSERT2(std::isfinite(power)); energy_source[i] -= power; // Note: Sign negative because power > 0 // Calculation of total heat flux for diagnostic purposes - q = (gamma_i * tisheath + 0.5 * Mi * C_i_sq) * nisheath * visheath; - flux = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); - power = flux / (coord->dy[i] * coord->J[i]); + q = gamma_i * tisheath * nisheath * visheath; + heatflow = q * (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im])); // W - ion_sheath_power_ylow[ip] += power; // Upper Y, so power placed in first guard cell + ion_sheath_power_ylow[ip] += heatflow; // Upper Y, so power placed in first guard cell } } From 5c42cfec203f5c8c876b108f60da6d84f6516c38 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 1 Oct 2023 14:53:57 +0100 Subject: [PATCH 05/19] sheath_boundary_simple: Shifting from field aligned Field energy_flow_ylow should be shifted to non-aligned coordinates. In axisymmetric simulations this doesn't make a difference. Thanks Ben! --- src/sheath_boundary_simple.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index d67395fb..f987d787 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -418,7 +418,7 @@ void SheathBoundarySimple::transform(Options& state) { set(electrons["energy_source"], fromFieldAligned(electron_energy_source)); // Add the total sheath power flux to the tracker of y power flows - add(electrons["energy_flow_ylow"], electron_sheath_power_ylow); + add(electrons["energy_flow_ylow"], fromFieldAligned(electron_sheath_power_ylow)); if (IS_SET_NOBOUNDARY(electrons["velocity"])) { setBoundary(electrons["velocity"], fromFieldAligned(Ve)); @@ -622,6 +622,6 @@ void SheathBoundarySimple::transform(Options& state) { set(species["energy_source"], fromFieldAligned(energy_source)); // Add the total sheath power flux to the tracker of y power flows - add(species["energy_flow_ylow"], ion_sheath_power_ylow); + add(species["energy_flow_ylow"], fromFieldAligned(ion_sheath_power_ylow)); } } From fe6161c2ac2504aa309ee8205cda1c12f83ab3ce Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 1 Oct 2023 17:41:34 +0100 Subject: [PATCH 06/19] Fix fast recycling neutral heat flow calc for target Heat flow now correctly calculated from ion sheath heat flow which is calculated in sheath_boundary_simple. sheath_boundary not yet supported. Improved variable names and comments for clarity --- include/recycling.hxx | 1 + src/recycling.cxx | 87 ++++++++++++++++++++++++++----------------- 2 files changed, 53 insertions(+), 35 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index 3892f4ab..384e990e 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -58,6 +58,7 @@ private: bool diagnose; ///< Save additional post-processing variables? Field3D density_source, energy_source; ///< Recycling particle and energy sources for all locations + Field3D heatflow_ylow, heatflow_xlow; ///< Cell edge fluxes used for calculating fast recycling energy source Field3D target_recycle_density_source, target_recycle_energy_source; ///< Recycling particle and energy sources for target recycling only Field3D wall_recycle_density_source, wall_recycle_energy_source; ///< Recycling particle and energy sources for pfr + sol recycling diff --git a/src/recycling.cxx b/src/recycling.cxx index 545abbf4..96c4a82b 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -156,6 +156,19 @@ void Recycling::transform(Options& state) { ? getNonFinal(species_to["energy_source"]) : 0.0; + // Fast recycling needs to know how much energy the "from" species is losing to the boundary + if (species_from.isSet("energy_flow_ylow")) { + heatflow_ylow = get(species_from["energy_flow_ylow"]); + } else if (channel.target_fast_recycle_fraction > 0) { + throw BoutException("Target fast recycle enabled but no sheath heat flow available, check compatibility with sheath component"); + }; + + // if (species.isSet("energy_flow_xlow")) { + // heatflow_xlow = get(species["energy_flow_xlow"]); + // } else if (sol_fast_recycle_fraction > 0) or (pfr_fast_recycle_fraction > 0) { + // throw BoutException("SOL/PFR fast recycle enabled but no wall heat flow available, check your wall BC choice"); + // }; + // Recycling at the divertor target plates @@ -178,30 +191,32 @@ void Recycling::transform(Options& state) { flux = 0.0; } - // Flow of recycled species inwards + // Flow of recycled neutrals into domain [s-1] BoutReal flow = channel.target_multiplier * flux * (J(r.ind, mesh->ystart) + J(r.ind, mesh->ystart - 1)) - / (sqrt(g_22(r.ind, mesh->ystart)) + sqrt(g_22(r.ind, mesh->ystart - 1))); + / (sqrt(g_22(r.ind, mesh->ystart)) + sqrt(g_22(r.ind, mesh->ystart - 1))); // Omitting *dx*dz because it cancels out when calculating source - // Add to density source - target_recycle_density_source(r.ind, mesh->ystart, jz) += flow + // Calculate sources in the final cell [m^-3 s^-1] + target_recycle_density_source(r.ind, mesh->ystart, jz) += flow // For diagnostic / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - density_source(r.ind, mesh->ystart, jz) += flow + density_source(r.ind, mesh->ystart, jz) += flow // For use in solver / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - // energy of recycled particles - // VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV - // TODO: Is the energy here 3eV * NV or 3/2 * 3eV * NV? It was just 3eV * NV before but it's adding to energy source - BoutReal tisheath = (T(r.ind, mesh->ystart, jz) + T(r.ind, mesh->ystart-1, jz)) * 0.5; // Ion temp at wall - BoutReal neutral_energy = flow * ( - channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction * tisheath // Fast recycling part - + (1 - channel.target_fast_recycle_fraction) * channel.target_energy); // Thermal recycling part + // Energy of recycled particles + BoutReal ion_heatflow = heatflow_ylow(r.ind, mesh->ystart, jz); // Ion heat flux to wall. This is ylow end so take first domain cell - target_recycle_energy_source(r.ind, mesh->ystart, jz) += neutral_energy - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - energy_source(r.ind, mesh->ystart, jz) += neutral_energy - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal neutral_heatflow = + ion_heatflow * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part + + // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] + target_recycle_energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow + / (J(r.ind, mesh->ystart) * dx(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart) * dz(r.ind, mesh->ystart)); + energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow + / (J(r.ind, mesh->ystart) * dx(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart) * dz(r.ind, mesh->ystart)); } } @@ -220,30 +235,32 @@ void Recycling::transform(Options& state) { flux = 0.0; } - // Flow of neutrals inwards + // Flow of recycled neutrals into domain [s-1] BoutReal flow = - channel.target_multiplier * flux * (J(r.ind, mesh->yend) + J(r.ind, mesh->yend + 1)) - / (sqrt(g_22(r.ind, mesh->yend)) + sqrt(g_22(r.ind, mesh->yend + 1))); + channel.target_multiplier * flux + * (J(r.ind, mesh->yend) + J(r.ind, mesh->yend + 1)) + / (sqrt(g_22(r.ind, mesh->yend)) + sqrt(g_22(r.ind, mesh->yend + 1))); // Omitting *dx*dz because it cancels out when calculating source - // Rate of change of neutrals in final cell - // Add to density source - target_recycle_density_source(r.ind, mesh->yend, jz) += flow - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - density_source(r.ind, mesh->yend, jz) += flow + // Calculate sources in the final cell [m^-3 s^-1] + target_recycle_density_source(r.ind, mesh->yend, jz) += flow // For diagnostic / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); + density_source(r.ind, mesh->yend, jz) += flow // For use in solver + / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - // energy of recycled particles - // VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV - // TODO: Is the energy here 3eV * NV or 3/2 * 3eV * NV? It was just 3eV * NV before but it's adding to energy source - BoutReal tisheath = (T(r.ind, mesh->yend, jz) + T(r.ind, mesh->yend+1, jz)) * 0.5; // Ion temp at wall - BoutReal neutral_energy = flow * ( - channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction * tisheath // Fast recycling part - + (1 - channel.target_fast_recycle_fraction) * channel.target_energy); // Thermal recycling part + BoutReal ion_heatflow = heatflow_ylow(r.ind, mesh->yend + 1, jz); // Ion heat flow to wall in [W]. This is yup end so take guard cell - target_recycle_energy_source(r.ind, mesh->yend, jz) += neutral_energy - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - energy_source(r.ind, mesh->yend, jz) += neutral_energy - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal neutral_heatflow = + ion_heatflow * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part + + + // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] + target_recycle_energy_source(r.ind, mesh->yend, jz) += neutral_heatflow + / (J(r.ind, mesh->yend) * dx(r.ind, mesh->yend) * dy(r.ind, mesh->yend) * dz(r.ind, mesh->yend)); + energy_source(r.ind, mesh->yend, jz) += neutral_heatflow + / (J(r.ind, mesh->yend) * dx(r.ind, mesh->yend) * dy(r.ind, mesh->yend) * dz(r.ind, mesh->yend)); } } } From 1d120d597316981a3ece341095f700e5470b2114 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 1 Oct 2023 18:33:47 +0100 Subject: [PATCH 07/19] Extend fast recycling to SOL and PFR recycling Further improvements to code clarity and readability --- include/recycling.hxx | 5 +-- src/recycling.cxx | 76 +++++++++++++++++++++++++++++-------------- 2 files changed, 55 insertions(+), 26 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index 384e990e..748af58e 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -58,12 +58,13 @@ private: bool diagnose; ///< Save additional post-processing variables? Field3D density_source, energy_source; ///< Recycling particle and energy sources for all locations - Field3D heatflow_ylow, heatflow_xlow; ///< Cell edge fluxes used for calculating fast recycling energy source + Field3D energy_flow_ylow, energy_flow_xlow; ///< Cell edge fluxes used for calculating fast recycling energy source + Field3D particle_flow_xlow; ///< Radial wall particle fluxes for recycling calc. No need to get poloidal from here, it's calculated from sheath velocity Field3D target_recycle_density_source, target_recycle_energy_source; ///< Recycling particle and energy sources for target recycling only Field3D wall_recycle_density_source, wall_recycle_energy_source; ///< Recycling particle and energy sources for pfr + sol recycling - Field3D radial_particle_outflow, radial_energy_outflow; ///< Radial fluxes coming from evolve_density and evolve_pressure used in recycling calc + // 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 96c4a82b..2418d874 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -156,24 +156,17 @@ void Recycling::transform(Options& state) { ? getNonFinal(species_to["energy_source"]) : 0.0; - // Fast recycling needs to know how much energy the "from" species is losing to the boundary - if (species_from.isSet("energy_flow_ylow")) { - heatflow_ylow = get(species_from["energy_flow_ylow"]); - } else if (channel.target_fast_recycle_fraction > 0) { - throw BoutException("Target fast recycle enabled but no sheath heat flow available, check compatibility with sheath component"); - }; - - // if (species.isSet("energy_flow_xlow")) { - // heatflow_xlow = get(species["energy_flow_xlow"]); - // } else if (sol_fast_recycle_fraction > 0) or (pfr_fast_recycle_fraction > 0) { - // throw BoutException("SOL/PFR fast recycle enabled but no wall heat flow available, check your wall BC choice"); - // }; - - // Recycling at the divertor target plates if (target_recycle) { + // Fast recycling needs to know how much energy the "from" species is losing to the boundary + if (species_from.isSet("energy_flow_ylow")) { + energy_flow_ylow = get(species_from["energy_flow_ylow"]); + } else if (channel.target_fast_recycle_fraction > 0) { + throw BoutException("Target fast recycle enabled but no sheath heat flow available, check compatibility with sheath component"); + }; + target_recycle_density_source = 0; target_recycle_energy_source = 0; @@ -204,7 +197,7 @@ void Recycling::transform(Options& state) { / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); // Energy of recycled particles - BoutReal ion_heatflow = heatflow_ylow(r.ind, mesh->ystart, jz); // Ion heat flux to wall. This is ylow end so take first domain cell + BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->ystart, jz); // Ion heat flux to wall. This is ylow end so take first domain cell // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] @@ -247,7 +240,7 @@ void Recycling::transform(Options& state) { density_source(r.ind, mesh->yend, jz) += flow // For use in solver / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - BoutReal ion_heatflow = heatflow_ylow(r.ind, mesh->yend + 1, jz); // Ion heat flow to wall in [W]. This is yup end so take guard cell + BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->yend + 1, jz); // Ion heat flow to wall in [W]. This is yup end so take guard cell // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] @@ -268,11 +261,29 @@ void Recycling::transform(Options& state) { wall_recycle_density_source = 0; wall_recycle_energy_source = 0; + // Get edge particle and heat for the species being recycled + if (sol_recycle or pfr_recycle) { + + if (species_from.isSet("energy_flow_xlow")) { + energy_flow_xlow = get(species_from["energy_flow_xlow"]); + } else if ((channel.sol_fast_recycle_fraction > 0) or (channel.pfr_fast_recycle_fraction > 0)) { + throw BoutException("SOL/PFR fast recycle enabled but no cell edge heat flow available, check your wall BC choice"); + }; + + if (species_from.isSet("particle_flow_xlow")) { + particle_flow_xlow = get(species_from["particle_flow_xlow"]); + } else if ((channel.sol_fast_recycle_fraction > 0) or (channel.pfr_fast_recycle_fraction > 0)) { + throw BoutException("SOL/PFR fast recycle enabled but no cell edge particle flow available, check your wall BC choice"); + }; + + } + // Recycling at the SOL edge (2D/3D only) if (sol_recycle) { // Flow out of domain is positive in the positive coordinate direction - radial_particle_outflow = get(species_from["particle_flow_xlow"]); + Field3D radial_particle_outflow = energy_flow_xlow; + Field3D radial_energy_outflow = particle_flow_xlow; if(mesh->lastX()){ // Only do this for the processor which has the edge region for(int iy=0; iy < mesh->LocalNy ; iy++){ @@ -294,9 +305,17 @@ void Recycling::transform(Options& state) { wall_recycle_density_source(mesh->xend, iy, iz) += recycle_particle_flow / volume; density_source(mesh->xend, iy, iz) += recycle_particle_flow / volume; - // For now, this is a fixed temperature - wall_recycle_energy_source(mesh->xend, iy, iz) += channel.sol_energy * recycle_particle_flow / volume; - energy_source(mesh->xend, iy, iz) += channel.sol_energy * recycle_particle_flow / volume; + BoutReal ion_heatflow = radial_energy_outflow(mesh->xend+1, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take guard cell + + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal neutral_heatflow = + ion_heatflow * channel.sol_fast_recycle_energy_factor * channel.sol_fast_recycle_fraction // Fast recycling part + + recycle_particle_flow * (1 - channel.sol_fast_recycle_fraction) * channel.sol_energy; // Thermal recycling part + + // Divide heat flow in [W] by cell volume to get energy source in [m^-3 s^-1] + wall_recycle_energy_source(mesh->xend, iy, iz) += neutral_heatflow / volume; // For diagnostic + energy_source(mesh->xend, iy, iz) += neutral_heatflow / volume; // For use in solver @@ -309,7 +328,8 @@ void Recycling::transform(Options& state) { if (pfr_recycle) { // PFR is flipped compared to edge: x=0 is at the PFR edge. Therefore outflow is in the negative coordinate direction. - radial_particle_outflow = get(species_from["particle_flow_xlow"]) * -1; + Field3D radial_particle_outflow = particle_flow_xlow * -1; + Field3D radial_energy_outflow = energy_flow_xlow * -1; if(mesh->firstX()){ // Only do this for the processor which has the core region if (!mesh->periodicY(mesh->xstart)) { // Only do this for the processor with a periodic Y, i.e. the PFR @@ -333,9 +353,17 @@ void Recycling::transform(Options& state) { wall_recycle_density_source(mesh->xstart, iy, iz) += recycle_particle_flow / volume; density_source(mesh->xstart, iy, iz) += recycle_particle_flow / volume; - // For now, this is a fixed temperature - wall_recycle_energy_source(mesh->xstart, iy, iz) += channel.pfr_energy * recycle_particle_flow / volume; - energy_source(mesh->xstart, iy, iz) += channel.pfr_energy * recycle_particle_flow / volume; + BoutReal ion_heatflow = radial_energy_outflow(mesh->xstart, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take first domain cell + + // Blend fast (ion energy) and thermal (constant energy) recycling + // Calculate returning neutral heat flow in [W] + BoutReal neutral_heatflow = + ion_heatflow * channel.pfr_fast_recycle_energy_factor * channel.pfr_fast_recycle_fraction // Fast recycling part + + recycle_particle_flow * (1 - channel.pfr_fast_recycle_fraction) * channel.pfr_energy; // Thermal recycling part + + // Divide heat flow in [W] by cell volume to get energy source in [m^-3 s^-1] + wall_recycle_energy_source(mesh->xstart, iy, iz) += neutral_heatflow / volume; // For diagnostic + energy_source(mesh->xstart, iy, iz) += neutral_heatflow / volume; // For use in solver } } From b8910c96eb4d39f29af07b6339c395d2bf182657 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 1 Oct 2023 19:04:22 +0100 Subject: [PATCH 08/19] Fix typo in SOL recycling --- include/recycling.hxx | 2 -- src/recycling.cxx | 4 ++-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index 748af58e..dce58fa0 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -64,8 +64,6 @@ private: Field3D target_recycle_density_source, target_recycle_energy_source; ///< Recycling particle and energy sources for target recycling only Field3D wall_recycle_density_source, wall_recycle_energy_source; ///< Recycling particle and energy sources for pfr + sol recycling - // 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 2418d874..b1401c95 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -282,8 +282,8 @@ void Recycling::transform(Options& state) { if (sol_recycle) { // Flow out of domain is positive in the positive coordinate direction - Field3D radial_particle_outflow = energy_flow_xlow; - Field3D radial_energy_outflow = particle_flow_xlow; + Field3D radial_particle_outflow = particle_flow_xlow; + Field3D radial_energy_outflow = energy_flow_xlow; if(mesh->lastX()){ // Only do this for the processor which has the edge region for(int iy=0; iy < mesh->LocalNy ; iy++){ From 7b8e0262cd390fbc7db95a3104719573f185b4f8 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 1 Oct 2023 19:52:50 +0100 Subject: [PATCH 09/19] Fix omission of rec multiplier in fast recycling --- src/recycling.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index b1401c95..c46043fd 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -202,7 +202,7 @@ void Recycling::transform(Options& state) { // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] BoutReal neutral_heatflow = - ion_heatflow * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + ion_heatflow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] @@ -245,7 +245,7 @@ void Recycling::transform(Options& state) { // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] BoutReal neutral_heatflow = - ion_heatflow * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + ion_heatflow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part @@ -310,7 +310,7 @@ void Recycling::transform(Options& state) { // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] BoutReal neutral_heatflow = - ion_heatflow * channel.sol_fast_recycle_energy_factor * channel.sol_fast_recycle_fraction // Fast recycling part + ion_heatflow * channel.sol_multiplier * channel.sol_fast_recycle_energy_factor * channel.sol_fast_recycle_fraction // Fast recycling part + recycle_particle_flow * (1 - channel.sol_fast_recycle_fraction) * channel.sol_energy; // Thermal recycling part // Divide heat flow in [W] by cell volume to get energy source in [m^-3 s^-1] @@ -358,7 +358,7 @@ void Recycling::transform(Options& state) { // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] BoutReal neutral_heatflow = - ion_heatflow * channel.pfr_fast_recycle_energy_factor * channel.pfr_fast_recycle_fraction // Fast recycling part + ion_heatflow * channel.pfr_multiplier * channel.pfr_fast_recycle_energy_factor * channel.pfr_fast_recycle_fraction // Fast recycling part + recycle_particle_flow * (1 - channel.pfr_fast_recycle_fraction) * channel.pfr_energy; // Thermal recycling part // Divide heat flow in [W] by cell volume to get energy source in [m^-3 s^-1] From 69226afee89b4a10553ef5f669cc0cdeb036538a Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 1 Oct 2023 20:47:08 +0100 Subject: [PATCH 10/19] Fix flow calculation to be actual flow and improve code readability --- src/recycling.cxx | 37 +++++++++++++++++-------------------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index c46043fd..ff62e9e4 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -187,14 +187,14 @@ void Recycling::transform(Options& state) { // Flow of recycled neutrals into domain [s-1] BoutReal flow = channel.target_multiplier * flux - * (J(r.ind, mesh->ystart) + J(r.ind, mesh->ystart - 1)) - / (sqrt(g_22(r.ind, mesh->ystart)) + sqrt(g_22(r.ind, mesh->ystart - 1))); // Omitting *dx*dz because it cancels out when calculating source + * (J(r.ind, mesh->ystart) + J(r.ind, mesh->ystart - 1)) / (sqrt(g_22(r.ind, mesh->ystart)) + sqrt(g_22(r.ind, mesh->ystart - 1))) + * 0.5*(dx(r.ind, mesh->ystart) + dx(r.ind, mesh->ystart - 1)) * 0.5*(dz(r.ind, mesh->ystart) + dz(r.ind, mesh->ystart - 1)); + + BoutReal volume = J(r.ind, mesh->ystart) * dx(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart) * dz(r.ind, mesh->ystart); // Calculate sources in the final cell [m^-3 s^-1] - target_recycle_density_source(r.ind, mesh->ystart, jz) += flow // For diagnostic - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); - density_source(r.ind, mesh->ystart, jz) += flow // For use in solver - / (J(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart)); + target_recycle_density_source(r.ind, mesh->ystart, jz) += flow / volume; // For diagnostic + density_source(r.ind, mesh->ystart, jz) += flow / volume; // For use in solver // Energy of recycled particles BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->ystart, jz); // Ion heat flux to wall. This is ylow end so take first domain cell @@ -206,10 +206,8 @@ void Recycling::transform(Options& state) { + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] - target_recycle_energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow - / (J(r.ind, mesh->ystart) * dx(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart) * dz(r.ind, mesh->ystart)); - energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow - / (J(r.ind, mesh->ystart) * dx(r.ind, mesh->ystart) * dy(r.ind, mesh->ystart) * dz(r.ind, mesh->ystart)); + target_recycle_energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow / volume; + energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow / volume; } } @@ -231,15 +229,16 @@ void Recycling::transform(Options& state) { // Flow of recycled neutrals into domain [s-1] BoutReal flow = channel.target_multiplier * flux - * (J(r.ind, mesh->yend) + J(r.ind, mesh->yend + 1)) - / (sqrt(g_22(r.ind, mesh->yend)) + sqrt(g_22(r.ind, mesh->yend + 1))); // Omitting *dx*dz because it cancels out when calculating source + * (J(r.ind, mesh->yend) + J(r.ind, mesh->yend + 1)) / (sqrt(g_22(r.ind, mesh->yend)) + sqrt(g_22(r.ind, mesh->yend + 1))) + * 0.5*(dx(r.ind, mesh->yend) + dx(r.ind, mesh->yend + 1)) * 0.5*(dz(r.ind, mesh->yend) + dz(r.ind, mesh->yend + 1)); + + BoutReal volume = J(r.ind, mesh->yend) * dx(r.ind, mesh->yend) * dy(r.ind, mesh->yend) * dz(r.ind, mesh->yend); // Calculate sources in the final cell [m^-3 s^-1] - target_recycle_density_source(r.ind, mesh->yend, jz) += flow // For diagnostic - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); - density_source(r.ind, mesh->yend, jz) += flow // For use in solver - / (J(r.ind, mesh->yend) * dy(r.ind, mesh->yend)); + target_recycle_density_source(r.ind, mesh->yend, jz) += flow / volume; // For diagnostic + density_source(r.ind, mesh->yend, jz) += flow / volume; // For use in solver + // Energy of recycled particles BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->yend + 1, jz); // Ion heat flow to wall in [W]. This is yup end so take guard cell // Blend fast (ion energy) and thermal (constant energy) recycling @@ -250,10 +249,8 @@ void Recycling::transform(Options& state) { // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] - target_recycle_energy_source(r.ind, mesh->yend, jz) += neutral_heatflow - / (J(r.ind, mesh->yend) * dx(r.ind, mesh->yend) * dy(r.ind, mesh->yend) * dz(r.ind, mesh->yend)); - energy_source(r.ind, mesh->yend, jz) += neutral_heatflow - / (J(r.ind, mesh->yend) * dx(r.ind, mesh->yend) * dy(r.ind, mesh->yend) * dz(r.ind, mesh->yend)); + target_recycle_energy_source(r.ind, mesh->yend, jz) += neutral_heatflow / volume; + energy_source(r.ind, mesh->yend, jz) += neutral_heatflow / volume; } } } From a231c2746651a6c9fd39bdd2667e8adb9690d271 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sun, 1 Oct 2023 21:10:07 +0100 Subject: [PATCH 11/19] Fix ion heat flow sign on ylow recycling --- src/recycling.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index ff62e9e4..17825c3e 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -197,7 +197,7 @@ void Recycling::transform(Options& state) { density_source(r.ind, mesh->ystart, jz) += flow / volume; // For use in solver // Energy of recycled particles - BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->ystart, jz); // Ion heat flux to wall. This is ylow end so take first domain cell + BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->ystart, jz) * -1; // This is ylow end so take first domain cell and flip sign // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] From 4a998a7e8d1b0a5b288439cd348aa13d782c236e Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Oct 2023 11:35:57 +0100 Subject: [PATCH 12/19] Add variables/options/dependencies for pump --- include/recycling.hxx | 8 +++++--- src/recycling.cxx | 25 +++++++++++++++++++++++-- 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index dce58fa0..df94c8a2 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -46,7 +46,7 @@ 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) BoutReal target_fast_recycle_fraction, pfr_fast_recycle_fraction, sol_fast_recycle_fraction; ///< Fraction of ions undergoing fast reflection BoutReal target_fast_recycle_energy_factor, sol_fast_recycle_energy_factor, pfr_fast_recycle_energy_factor; ///< Fraction of energy retained by fast recycled neutrals @@ -54,16 +54,18 @@ private: 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 energy_flow_ylow, energy_flow_xlow; ///< Cell edge fluxes used for calculating fast recycling energy source Field3D particle_flow_xlow; ///< Radial wall particle fluxes for recycling calc. No need to get poloidal from here, it's calculated from sheath velocity + Field2D is_pump; ///< 1 = pump, 0 = no pump. Works only in SOL/PFR. Provided by user in grid file. + Field3D target_recycle_density_source, target_recycle_energy_source; ///< Recycling particle and energy sources for target recycling only Field3D wall_recycle_density_source, wall_recycle_energy_source; ///< Recycling particle and energy sources for pfr + sol recycling - + Field3D pump_recycle_density_source, pump_recycle_energy_source; ///< Recycling particle and energy sources for pump recycling }; diff --git a/src/recycling.cxx b/src/recycling.cxx index 17825c3e..de00311e 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 @@ -21,6 +22,12 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) { .as(), ','); + + // Neutral pump + // 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")); for (const auto& species : species_list) { std::string from = trim(species, " \t\r()"); // The species name in the list @@ -52,6 +59,11 @@ 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) @@ -99,14 +111,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)) { + 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("recycle_fraction 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, target_fast_recycle_fraction, pfr_fast_recycle_fraction, sol_fast_recycle_fraction, target_fast_recycle_energy_factor, sol_fast_recycle_energy_factor, pfr_fast_recycle_energy_factor}); @@ -123,6 +136,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); } } @@ -146,6 +163,10 @@ void Recycling::transform(Options& state) { const Field3D T = get(species_from["temperature"]); // Ion temperature 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 From 016806590dfa31869b4f72cb0672aad47ced77c5 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Oct 2023 14:50:32 +0100 Subject: [PATCH 13/19] Implement pump from master to help merge --- include/recycling.hxx | 2 +- src/recycling.cxx | 207 +++++++++++++++++++++++++++++++++++------- 2 files changed, 174 insertions(+), 35 deletions(-) diff --git a/include/recycling.hxx b/include/recycling.hxx index df94c8a2..8b37d322 100644 --- a/include/recycling.hxx +++ b/include/recycling.hxx @@ -65,7 +65,7 @@ private: Field3D target_recycle_density_source, target_recycle_energy_source; ///< Recycling particle and energy sources for target recycling only Field3D wall_recycle_density_source, wall_recycle_energy_source; ///< Recycling particle and energy sources for pfr + sol recycling - Field3D pump_recycle_density_source, pump_recycle_energy_source; ///< Recycling particle and energy sources for pump recycling + Field3D pump_density_source, pump_energy_source; ///< Recycling particle and energy sources for pump recycling }; diff --git a/src/recycling.cxx b/src/recycling.cxx index de00311e..cd453ef8 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -218,17 +218,17 @@ void Recycling::transform(Options& state) { density_source(r.ind, mesh->ystart, jz) += flow / volume; // For use in solver // Energy of recycled particles - BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->ystart, jz) * -1; // This is ylow end so take first domain cell and flip sign + BoutReal ion_energy_flow = energy_flow_ylow(r.ind, mesh->ystart, jz) * -1; // This is ylow end so take first domain cell and flip sign // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] - BoutReal neutral_heatflow = - ion_heatflow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + BoutReal recycle_energy_flow = + ion_energy_flow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] - target_recycle_energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow / volume; - energy_source(r.ind, mesh->ystart, jz) += neutral_heatflow / volume; + target_recycle_energy_source(r.ind, mesh->ystart, jz) += recycle_energy_flow / volume; + energy_source(r.ind, mesh->ystart, jz) += recycle_energy_flow / volume; } } @@ -260,24 +260,26 @@ void Recycling::transform(Options& state) { density_source(r.ind, mesh->yend, jz) += flow / volume; // For use in solver // Energy of recycled particles - BoutReal ion_heatflow = energy_flow_ylow(r.ind, mesh->yend + 1, jz); // Ion heat flow to wall in [W]. This is yup end so take guard cell + BoutReal ion_energy_flow = energy_flow_ylow(r.ind, mesh->yend + 1, jz); // Ion heat flow to wall in [W]. This is yup end so take guard cell // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] - BoutReal neutral_heatflow = - ion_heatflow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + BoutReal recycle_energy_flow = + ion_energy_flow * channel.target_multiplier * channel.target_fast_recycle_energy_factor * channel.target_fast_recycle_fraction // Fast recycling part + flow * (1 - channel.target_fast_recycle_fraction) * channel.target_energy; // Thermal recycling part // Divide heat flow in [W] by cell volume to get source in [m^-3 s^-1] - target_recycle_energy_source(r.ind, mesh->yend, jz) += neutral_heatflow / volume; - energy_source(r.ind, mesh->yend, jz) += neutral_heatflow / volume; + target_recycle_energy_source(r.ind, mesh->yend, jz) += recycle_energy_flow / volume; + energy_source(r.ind, mesh->yend, jz) += recycle_energy_flow / volume; } } } wall_recycle_density_source = 0; wall_recycle_energy_source = 0; + pump_density_source = 0; + pump_energy_source = 0; // Get edge particle and heat for the species being recycled if (sol_recycle or pfr_recycle) { @@ -311,7 +313,13 @@ void Recycling::transform(Options& state) { BoutReal volume = J(mesh->xend, iy) * dx(mesh->xend, iy) * dy(mesh->xend, iy) * dz(mesh->xend, iy); - // Flow of recycled species back from the edge + // If cell is a pump, overwrite multiplier with pump multiplier + BoutReal multiplier = channel.sol_multiplier; + if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { + multiplier = channel.pump_multiplier; + } + + // Flow of recycled species back from the edge due to recycling // 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 BoutReal recycle_particle_flow = 0; @@ -319,23 +327,72 @@ void Recycling::transform(Options& state) { recycle_particle_flow = channel.sol_multiplier * radial_particle_outflow(mesh->xend+1, iy, iz); } - // Divide by volume to get source - wall_recycle_density_source(mesh->xend, iy, iz) += recycle_particle_flow / volume; - density_source(mesh->xend, iy, iz) += recycle_particle_flow / volume; - - BoutReal ion_heatflow = radial_energy_outflow(mesh->xend+1, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take guard cell + BoutReal ion_energy_flow = radial_energy_outflow(mesh->xend+1, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take guard cell // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] - BoutReal neutral_heatflow = - ion_heatflow * channel.sol_multiplier * channel.sol_fast_recycle_energy_factor * channel.sol_fast_recycle_fraction // Fast recycling part + BoutReal recycle_energy_flow = + ion_energy_flow * channel.sol_multiplier * channel.sol_fast_recycle_energy_factor * channel.sol_fast_recycle_fraction // Fast recycling part + recycle_particle_flow * (1 - channel.sol_fast_recycle_fraction) * channel.sol_energy; // Thermal recycling part - // Divide heat flow in [W] by cell volume to get energy source in [m^-3 s^-1] - wall_recycle_energy_source(mesh->xend, iy, iz) += neutral_heatflow / volume; // For diagnostic - energy_source(mesh->xend, iy, iz) += neutral_heatflow / volume; // For use in solver + // Calculate neutral pump neutral sinks due to neutral impingement + // These are additional to particle sinks due to recycling + BoutReal pump_neutral_energy_sink = 0; + BoutReal pump_neutral_particle_sink = 0; + + if ((is_pump(mesh->xend, iy) == 1.0) and (neutral_pump)) { + + 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 + + // 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 = 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 + // 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 pump_neutral_particle_flow = v_th * nnsheath * dasheath; // [s^-1] + pump_neutral_particle_sink = pump_neutral_particle_flow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + + // Use gamma=2.0 as per Stangeby p.69, total energy of static Maxwellian + BoutReal pump_neutral_energy_flow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] + pump_neutral_energy_sink = pump_neutral_energy_flow / dv * (1 - multiplier); // heatsink [W m^-3] + + // Divide flows by volume to get sources + // Save to pump diagnostic + pump_density_source(mesh->xend, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + pump_energy_source(mesh->xend, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + + + } else { + // Save to wall diagnostic (pump sinks are 0 if not on pump) + wall_recycle_density_source(mesh->xend, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + wall_recycle_energy_source(mesh->xend, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + } + + // Save to solver + density_source(mesh->xend, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + energy_source(mesh->xend, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; - } } @@ -354,12 +411,17 @@ void Recycling::transform(Options& state) { for(int iy=0; iy < mesh->LocalNy ; iy++){ for(int iz=0; iz < mesh->LocalNz; iz++){ - // Volume of cell adjacent to wall which will receive source BoutReal volume = J(mesh->xstart, iy) * dx(mesh->xstart, iy) * dy(mesh->xstart, iy) * dz(mesh->xstart, iy); - // Flow of recycled species back from the edge + // If cell is a pump, overwrite multiplier with pump multiplier + BoutReal multiplier = channel.pfr_multiplier; + if ((is_pump(mesh->xstart, iy) == 1.0) and (neutral_pump)) { + multiplier = channel.pump_multiplier; + } + + // Flow of recycled species back from the edge due to recycling // 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; @@ -367,26 +429,79 @@ void Recycling::transform(Options& state) { recycle_particle_flow = channel.pfr_multiplier * radial_particle_outflow(mesh->xstart, iy, iz); } - // Divide by volume to get source - wall_recycle_density_source(mesh->xstart, iy, iz) += recycle_particle_flow / volume; - density_source(mesh->xstart, iy, iz) += recycle_particle_flow / volume; - - BoutReal ion_heatflow = radial_energy_outflow(mesh->xstart, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take first domain cell + BoutReal ion_energy_flow = radial_energy_outflow(mesh->xstart, iy, iz); // Ion heat flow to wall in [W]. This is on xlow edge so take first domain cell // Blend fast (ion energy) and thermal (constant energy) recycling // Calculate returning neutral heat flow in [W] - BoutReal neutral_heatflow = - ion_heatflow * channel.pfr_multiplier * channel.pfr_fast_recycle_energy_factor * channel.pfr_fast_recycle_fraction // Fast recycling part + BoutReal recycle_energy_flow = + ion_energy_flow * channel.pfr_multiplier * channel.pfr_fast_recycle_energy_factor * channel.pfr_fast_recycle_fraction // Fast recycling part + recycle_particle_flow * (1 - channel.pfr_fast_recycle_fraction) * channel.pfr_energy; // Thermal recycling part - // Divide heat flow in [W] by cell volume to get energy source in [m^-3 s^-1] - wall_recycle_energy_source(mesh->xstart, iy, iz) += neutral_heatflow / volume; // For diagnostic - energy_source(mesh->xstart, iy, iz) += neutral_heatflow / volume; // For use in solver + // Calculate neutral pump neutral sinks due to neutral impingement + // These are additional to particle sinks due to recycling + BoutReal pump_neutral_energy_sink = 0; + BoutReal pump_neutral_particle_sink = 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 + 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 = 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 + // 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 pump_neutral_particle_flow = v_th * nnsheath * dasheath; // [s^-1] + pump_neutral_particle_sink = pump_neutral_particle_flow / dv * (1 - multiplier); // Particle sink [s^-1 m^-3] + + // Use gamma=2.0 as per Stangeby p.69, total energy of static Maxwellian + BoutReal pump_neutral_energy_flow = 2.0 * tnsheath * v_th * nnsheath * dasheath; // [W] + pump_neutral_energy_sink = pump_neutral_energy_flow / dv * (1 - multiplier); // heatsink [W m^-3] + + // Divide flows by volume to get sources + // Save to pump diagnostic + pump_density_source(mesh->xstart, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + pump_energy_source(mesh->xstart, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + + + } else { + // Save to wall diagnostic (pump sinks are 0 if not on pump) + wall_recycle_density_source(mesh->xstart, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + wall_recycle_energy_source(mesh->xstart, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; + + } + + // Save to solver + density_source(mesh->xstart, iy, iz) += recycle_particle_flow/volume - pump_neutral_particle_sink; + energy_source(mesh->xstart, iy, iz) += recycle_energy_flow/volume - pump_neutral_energy_sink; } } } } + } // Put the updated sources back into the state @@ -448,6 +563,30 @@ void Recycling::outputVars(Options& state) { {"long_name", std::string("Wall 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_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_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"}}); + + 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 b649da0392e55c89e306c6baa6e8fc9a36c66ce4 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Oct 2023 15:17:24 +0100 Subject: [PATCH 14/19] Fix typo in neutral pump --- src/recycling.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index 456bbf3b..51db51dd 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -278,8 +278,6 @@ void Recycling::transform(Options& state) { } // Initialise counters of pump recycling fluxes - pump_recycle_density_source = 0; - pump_recycle_energy_source = 0; wall_recycle_density_source = 0; wall_recycle_energy_source = 0; pump_density_source = 0; From 58a8c8435848fdadac8daf1ba51e278ff3f70c0c Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Oct 2023 20:01:44 +0100 Subject: [PATCH 15/19] Fix segmentation fault when using sheath_boundary --- src/recycling.cxx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index 51db51dd..68cb7e48 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -185,9 +185,12 @@ void Recycling::transform(Options& state) { // Fast recycling needs to know how much energy the "from" species is losing to the boundary if (species_from.isSet("energy_flow_ylow")) { energy_flow_ylow = get(species_from["energy_flow_ylow"]); - } else if (channel.target_fast_recycle_fraction > 0) { - throw BoutException("Target fast recycle enabled but no sheath heat flow available, check compatibility with sheath component"); - }; + } else { + energy_flow_ylow = 0; + if (channel.target_fast_recycle_fraction > 0) { + throw BoutException("Target fast recycle enabled but no sheath heat flow available in energy_flow_ylow! \nCurrently only sheath_boundary_simple is supported for fast recycling."); + } + } target_recycle_density_source = 0; target_recycle_energy_source = 0; From 319f62c272182783232d504779bd685a232f3da9 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 3 Oct 2023 09:49:59 +0100 Subject: [PATCH 16/19] Update 1D example recycling flags --- examples/1D-hydrogen-equalT/BOUT.inp | 2 ++ examples/1D-hydrogen/BOUT.inp | 2 ++ examples/1D-neon-source/BOUT.inp | 23 +++++++++++++++++- examples/1D-neon/BOUT.inp | 24 ++++++++++++++++++- examples/1D-recycling/BOUT.inp | 2 ++ .../single-temperature/BOUT.inp | 4 +++- .../two-temperatures/BOUT.inp | 4 +++- 7 files changed, 57 insertions(+), 4 deletions(-) diff --git a/examples/1D-hydrogen-equalT/BOUT.inp b/examples/1D-hydrogen-equalT/BOUT.inp index b164d07f..37eaf69c 100644 --- a/examples/1D-hydrogen-equalT/BOUT.inp +++ b/examples/1D-hydrogen-equalT/BOUT.inp @@ -95,6 +95,8 @@ density_controller_p = 5e2 diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] diff --git a/examples/1D-hydrogen/BOUT.inp b/examples/1D-hydrogen/BOUT.inp index 889602b3..51b0d3d6 100644 --- a/examples/1D-hydrogen/BOUT.inp +++ b/examples/1D-hydrogen/BOUT.inp @@ -98,6 +98,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] diff --git a/examples/1D-neon-source/BOUT.inp b/examples/1D-neon-source/BOUT.inp index ba0b3cc5..0204a8fc 100644 --- a/examples/1D-neon-source/BOUT.inp +++ b/examples/1D-neon-source/BOUT.inp @@ -103,6 +103,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] @@ -169,6 +171,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+] @@ -189,6 +193,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+2] @@ -209,6 +215,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+3] @@ -229,6 +237,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -251,8 +261,9 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction - diagnose = true [Nne+5] @@ -273,6 +284,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -295,6 +308,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -317,6 +332,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -339,6 +356,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -361,6 +380,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true diff --git a/examples/1D-neon/BOUT.inp b/examples/1D-neon/BOUT.inp index 3f4b4d19..2e0b9376 100644 --- a/examples/1D-neon/BOUT.inp +++ b/examples/1D-neon/BOUT.inp @@ -104,7 +104,9 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d -target_recycle_multiplier = 1.0 # Recycling fraction +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy +target_recycle_multiplier = 1 # Recycling fraction [Nd+] @@ -168,6 +170,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction [Nne+] @@ -188,6 +192,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction [Nne+2] @@ -208,6 +214,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction [Nne+3] @@ -228,6 +236,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -250,6 +260,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -272,6 +284,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -294,6 +308,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -316,6 +332,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -338,6 +356,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true @@ -360,6 +380,8 @@ noflow_lower_y = true noflow_upper_y = false # Sheath boundary at upper y recycle_as = ne +target_recycle = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 1 # Recycling fraction diagnose = true diff --git a/examples/1D-recycling/BOUT.inp b/examples/1D-recycling/BOUT.inp index d3ca2049..73666e46 100644 --- a/examples/1D-recycling/BOUT.inp +++ b/examples/1D-recycling/BOUT.inp @@ -101,6 +101,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true recycle_as = d +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1 # Recycling fraction [Nd+] diff --git a/examples/solkit-comparison/single-temperature/BOUT.inp b/examples/solkit-comparison/single-temperature/BOUT.inp index 475b0efb..342e3f4e 100644 --- a/examples/solkit-comparison/single-temperature/BOUT.inp +++ b/examples/solkit-comparison/single-temperature/BOUT.inp @@ -100,7 +100,9 @@ AA = 2 diagnose = true # Output diagnostics for these components? recycle_as = d -target_recycle_multiplier = 1.0 # Recycling fraction +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy +target_recycle_multiplier = 1 # Recycling fraction [Nd+] diff --git a/examples/solkit-comparison/two-temperatures/BOUT.inp b/examples/solkit-comparison/two-temperatures/BOUT.inp index 6f429253..b7c91ee8 100644 --- a/examples/solkit-comparison/two-temperatures/BOUT.inp +++ b/examples/solkit-comparison/two-temperatures/BOUT.inp @@ -100,7 +100,9 @@ AA = 2 diagnose = true # Output diagnostics for these components? recycle_as = d -target_recycle_multiplier = 1.0 # Recycling fraction +target_recycle = true # Target recycling on +target_recycle_energy = 3.5 # Franck-Condon dissociation energy +target_recycle_multiplier = 1 # Recycling fraction [Nd+] From 3874652ab5acb63a1fec576334fd80dd81a81c44 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 3 Oct 2023 09:53:04 +0100 Subject: [PATCH 17/19] Update recycling flags for further examples --- examples/1D-hydrogen/qn/BOUT.inp | 1 + examples/1D-recycling/cvode/BOUT.inp | 2 ++ examples/tokamak/recycling-dthe-drifts/BOUT.inp | 3 +++ examples/tokamak/recycling-dthe/BOUT.inp | 3 +++ examples/tokamak/recycling-dthene/BOUT.inp | 4 ++++ examples/tokamak/recycling/BOUT.inp | 1 + 6 files changed, 14 insertions(+) diff --git a/examples/1D-hydrogen/qn/BOUT.inp b/examples/1D-hydrogen/qn/BOUT.inp index e1b020d1..37d8ca25 100644 --- a/examples/1D-hydrogen/qn/BOUT.inp +++ b/examples/1D-hydrogen/qn/BOUT.inp @@ -100,6 +100,7 @@ thermal_conduction = true # in evolve_pressure diagnose = true # Output diagnostics for these components? recycle_as = d +target_recycle = true target_recycle_multiplier = 1.0 # Recycling fraction [Nd+] diff --git a/examples/1D-recycling/cvode/BOUT.inp b/examples/1D-recycling/cvode/BOUT.inp index f0f11284..03713ba0 100644 --- a/examples/1D-recycling/cvode/BOUT.inp +++ b/examples/1D-recycling/cvode/BOUT.inp @@ -92,6 +92,8 @@ thermal_conduction = true # in evolve_pressure diagnose = true recycle_as = d +target_recycle = true +target_recycle_energy = 3.5 # Franck-Condon dissociation energy target_recycle_multiplier = 1 # Recycling fraction [Nd+] diff --git a/examples/tokamak/recycling-dthe-drifts/BOUT.inp b/examples/tokamak/recycling-dthe-drifts/BOUT.inp index 01cbe8bc..292892ce 100644 --- a/examples/tokamak/recycling-dthe-drifts/BOUT.inp +++ b/examples/tokamak/recycling-dthe-drifts/BOUT.inp @@ -53,6 +53,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nd+] @@ -89,6 +90,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = t +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nt+] @@ -125,6 +127,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = he +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nhe+] diff --git a/examples/tokamak/recycling-dthe/BOUT.inp b/examples/tokamak/recycling-dthe/BOUT.inp index ced75dd8..74fa2eed 100644 --- a/examples/tokamak/recycling-dthe/BOUT.inp +++ b/examples/tokamak/recycling-dthe/BOUT.inp @@ -53,6 +53,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nd+] @@ -89,6 +90,7 @@ anomalous_chi = 1.0 thermal_conduction = true recycle_as = t +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nt+] @@ -125,6 +127,7 @@ anomalous_chi = 1.0 thermal_conduction = true recycle_as = he +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true diff --git a/examples/tokamak/recycling-dthene/BOUT.inp b/examples/tokamak/recycling-dthene/BOUT.inp index 61d87c1f..21c65bb0 100644 --- a/examples/tokamak/recycling-dthene/BOUT.inp +++ b/examples/tokamak/recycling-dthene/BOUT.inp @@ -36,6 +36,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nd+] @@ -76,6 +77,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = t +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nt+] @@ -116,6 +118,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = he +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nhe+] @@ -156,6 +159,7 @@ anomalous_chi = 1 thermal_conduction = true recycle_as = ne +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction [Nne+] diff --git a/examples/tokamak/recycling/BOUT.inp b/examples/tokamak/recycling/BOUT.inp index 57a9b35d..310d0ee0 100644 --- a/examples/tokamak/recycling/BOUT.inp +++ b/examples/tokamak/recycling/BOUT.inp @@ -59,6 +59,7 @@ anomalous_chi = 4 thermal_conduction = true recycle_as = d +target_recycle = true target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true From 62d413829a6e424a2d63642f292477831386d323 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 5 Oct 2023 18:08:41 +0100 Subject: [PATCH 18/19] Add documentation for fast recycling --- docs/sphinx/components.rst | 51 ++++++++++++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 5 deletions(-) diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index 0cb370a7..b7a76614 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -1150,12 +1150,19 @@ Recycling has been implemented at the target, the SOL edge and the PFR edge. Each is off by default and must be activated with a separate flag. Each can be assigned a separate recycle multiplier and recycle energy. -The chosen species must feature an outflow through the boundary - any cells +Configuring thermal recycling +^^^^^^^^^^^^^^^ + +A simple and commonly used way to model recycling is to assume it is fully thermal, +i.e. that every incident ion recombines into a neutral molecule and thermalises with the surface +before becoming re-emitted. Hermes-3 does not yet have a hydrogenic molecule model, and so +the molecules are assumed to instantly dissociate at the Franck-Condon dissociation temperature of 3.5eV. + +In order to set this up, the chosen species must feature an outflow through the boundary - any cells with an inflow have their recycling source set to zero. If a sheath boundary condition is enabled, then this is automatically satisfied at the target through the Bohm condition. If it is not enabled, then the target boundary must be set to `free_o2`, `free_o3` or `decaylength` to -allow an outflow. If SOL or PFR recycling is enabled, a `free_o2`, `free_o3` or `decaylength`on -their respective boundaries is required at all times. +allow an outflow. The recycling component has a `species` option, that is a list of species to recycle. For each of the species in that list, `recycling` will look in @@ -1193,6 +1200,40 @@ Each returning atom has an energy of 3.5eV: pfr_recycle_multiplier = 1 # Recycling fraction pfr_recycle_energy = 3.5 # Energy of recycled particles [eV] +Allowing for fast recycling +^^^^^^^^^^^^^^^ + +In reality, a fraction of incident ions will undergo specular reflection off the surface and +preserve a fraction of their energy. In the popular Monte-Carlo neutral code EIRENE, the +fast recycling fraction and the energy reflection factor are provided by the `TRIM database `_ +as a function of incident angle, surface material and incident particle energy. +Studies found that sheath acceleration can make the ion angle relatively consistent, e.g. 60 degrees; in (`Jae-Sun Park et al 2021 Nucl. Fusion 61 016021 `_). + +The recycled heat flux is: + +.. math:: + + \begin{aligned} + \Gamma_{E_{n}} &= R \times (R_{f} \alpha_{E} \Gamma_{E_{i}}^{sheath} + (1 - R_{f} T_{R} \Gamma_{N_{i}})) \\ + \end{aligned} + +Where :math:`R` is the recycle multiplier, :math:`R_{f}` is the fast reflection fraction, :math:`\alpha_{E}` is the energy reflection factor, +:math:`\Gamma_{E_{i}}^{sheath}` is the incident heat flux from the sheath boundary condition, :math:`T_{R}` is the recycle energy. + +:math:`R_{f}` and :math:`\alpha_{E}` can be set as in the below example. They can also be set to different values for the SOL and PFR by replacing +the word "target" with either "sol" or "pfr". + +.. code-block:: ini + + [d+] + recycle_as = d # Species to recycle as + + target_recycle = true + target_recycle_multiplier = 0.95 # Recycling fraction + target_recycle_energy = 3.5 # Energy of recycled particles [eV] + target_fast_recycle_energy_factor = 0.70 + target_fast_recycle_fraction = 0.80 + Neutral pump ^^^^^^^^^^^^^^^ @@ -1204,8 +1245,8 @@ 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}` -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`. +which is set by the `pump_multiplier` option in the input file. The unrecycled ion flux :math:`\Gamma_{N_{i}}^{unrecycled}` is calculated using the recycling +model and allows for either thermal or fast recycling, but with the difference that the `pump_multiplier` replaces the `recycle_multiplier`. .. math:: From 550ac53f88b116c00561393483bed1f483577bf3 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 5 Oct 2023 18:12:32 +0100 Subject: [PATCH 19/19] Fix minor typo in fast recycling documentation --- docs/sphinx/components.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index b7a76614..64e4a29b 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -1214,11 +1214,11 @@ The recycled heat flux is: .. math:: \begin{aligned} - \Gamma_{E_{n}} &= R \times (R_{f} \alpha_{E} \Gamma_{E_{i}}^{sheath} + (1 - R_{f} T_{R} \Gamma_{N_{i}})) \\ + \Gamma_{E_{n}} &= R \times (R_{f} \alpha_{E} \Gamma_{E_{i}}^{sheath} + (1 - R_{f}) T_{R} \Gamma_{N_{i}})) \\ \end{aligned} Where :math:`R` is the recycle multiplier, :math:`R_{f}` is the fast reflection fraction, :math:`\alpha_{E}` is the energy reflection factor, -:math:`\Gamma_{E_{i}}^{sheath}` is the incident heat flux from the sheath boundary condition, :math:`T_{R}` is the recycle energy. +:math:`\Gamma_{E_{i}}^{sheath}` is the incident heat flux from the sheath boundary condition, :math:`T_{R}` is the recycle energy and :math:`\Gamma_{N_{i}}` is the incident ion flux. :math:`R_{f}` and :math:`\alpha_{E}` can be set as in the below example. They can also be set to different values for the SOL and PFR by replacing the word "target" with either "sol" or "pfr".