Merge branch 'master' into fast-recycling

docs/sphinx/components.rst

@@ -1189,9 +1189,57 @@ 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 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`.
+
+.. 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:

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))

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
+