diff --git a/examples/tokamak-1D/extra/1D-hydrogen/BOUT.inp b/examples/tokamak-1D/extra/1D-hydrogen/BOUT.inp index 0f2730928..449ee29a2 100644 --- a/examples/tokamak-1D/extra/1D-hydrogen/BOUT.inp +++ b/examples/tokamak-1D/extra/1D-hydrogen/BOUT.inp @@ -18,25 +18,31 @@ timestep = 10000 MXG = 0 # No guard cells in X [mesh] +# 1D simulation, use "y" as the dimension along the fieldline nx = 1 ny = 400 # Resolution along field-line nz = 1 length = 30 # Length of the domain in meters -length_xpt = 10 # Length from midplane to X-point [m] +length_xpt = 10 # Length from midplane to X-point [m] (i.e. this is where the source ends) dymin = 0.1 # Minimum grid spacing near target, as fraction of average. Must be > 0 and < 1 -# Parallel grid spacing +# Parallel grid spacing — grid refinement near the divertor target dy = (length / ny) * (1 + (1-dymin)*(1-y/pi)) -# Calculate where the source ends in grid index +# Calculate where the source ends in grid index (i.e. at the X-point) source = length_xpt / length y_xpt = pi * ( 2 - dymin - sqrt( (2-dymin)^2 - 4*(1-dymin)*source ) ) / (1 - dymin) +# Separatrix index set to -1 for 1D simulations ixseps1 = -1 ixseps2 = -1 +# In 1D, the Jacobian is inversely proportional to volume +# Total flux expansion can be set up with a non-constant J profile +J = 1 + [hermes] # Notes: # - electrons after other species, so the density can be set by quasineutrality @@ -57,14 +63,26 @@ Tnorm = 100 [solver] type = beuler # Backward Euler steady-state solver -snes_type = newtonls +timestep = 1e-2 # Initial timestep + +#### FIRST ORDER IN TIME, FASTER SOLVER BEULER +type = beuler # Backward Euler steady-state solver +snes_type = newtonls # Newton iterations +ksp_type = gmres # GMRES method +pc_type = hypre # Preconditioner type +pc_hypre_type = euclid # Hypre preconditioner type +max_nonlinear_iterations = 10 # Max Newton iterations per timestep +lag_jacobian = 500 # How long to wait before Jacobian recalculation per Newton iteration +#### + +# Solver tolerances +atol = 1e-7 # Absolute tolerance (controls small numbers) +rtol = 1e-5 # Relative tolerance (primary convergence control) -max_nonlinear_iterations = 10 +matrix_free_operator = true diagnose = false # Print diagnostic information? -atol = 1e-7 -rtol = 1e-5 [sheath_boundary] @@ -72,7 +90,7 @@ lower_y = false upper_y = true [neutral_parallel_diffusion] - +diffusion_collisions_mode = multispecies # "afn" or "multispecies" dneut = 10 # (B / Bpol)^2 in neutral diffusion terms #################################### @@ -108,15 +126,18 @@ function = 1 source_shape = H(mesh:y_xpt - y)*1e20 # For upstream_density_feedback [Pd+] + +# Initial condition for ion pressure (in terms of hermes:Nnorm * hermes:Tnorm) function = 1 powerflux = 2.5e7 # Input power flux in W/m^2 +# Divide by length of source region to get heat source in W/m3 source = (powerflux*2/3 / (mesh:length_xpt))*H(mesh:y_xpt - y) # Input power as function of y [NVd+] -function = 0 +function = 0 # Initialise with Vi=0 #################################### @@ -130,38 +151,35 @@ AA = 2 thermal_conduction = true [Nd] - function = 0.001 [Pd] - function = 0.0001 #################################### [e] # Electrons -type = quasineutral, zero_current, evolve_pressure, noflow_boundary +type = quasineutral, evolve_pressure, zero_current, noflow_boundary noflow_upper_y = false charge = -1 AA = 1/1836 - thermal_conduction = true # in evolve_pressure +diagnose = true [Pe] function = `Pd+:function` # Same as ion pressure initially - -source = `Pd+:source` # Same as ion pressure source +source = `Pd+:source` # Same as ion pressure source #################################### [recycling] - species = d+ [reactions] +diagnose = true type = ( d + e -> d+ + 2e, # Deuterium ionisation d+ + e -> d, # Deuterium recombination diff --git a/examples/tokamak-1D/extra/1D-neon-source/BOUT.inp b/examples/tokamak-1D/extra/1D-neon-source/BOUT.inp index 05a4f1dd8..72a019346 100644 --- a/examples/tokamak-1D/extra/1D-neon-source/BOUT.inp +++ b/examples/tokamak-1D/extra/1D-neon-source/BOUT.inp @@ -21,27 +21,30 @@ timestep = 1000 MXG = 0 # No guard cells in X [mesh] +# 1D simulation, use "y" as the dimension along the fieldline nx = 1 ny = 400 # Resolution along field-line nz = 1 length = 30 # Length of the domain in meters -length_xpt = 10 # Length from midplane to X-point [m] +length_xpt = 10 # Length from midplane to X-point [m] (i.e. this is where the source ends) dymin = 0.1 # Minimum grid spacing near target, as fraction of average. Must be > 0 and < 1 -# Parallel grid spacing +# Parallel grid spacing — grid refinement near the divertor target dy = (length / ny) * (1 + (1-dymin)*(1-y/pi)) -# Calculate where the source ends in grid index +# Calculate where the source ends in grid index (i.e. at the X-point) source = length_xpt / length y_xpt = pi * ( 2 - dymin - sqrt( (2-dymin)^2 - 4*(1-dymin)*source ) ) / (1 - dymin) +# Separatrix index set to -1 for 1D simulations ixseps1 = -1 ixseps2 = -1 -[restart_files] -init_missing = true # Initialise any missing variables to zero +# In 1D, the Jacobian is inversely proportional to volume +# Total flux expansion can be set up with a non-constant J profile +J = 1 [hermes] # Notes: @@ -50,10 +53,10 @@ init_missing = true # Initialise any missing variables to zero # - electron_force_balance after thermal force and collisions, so the electric field # includes the force on the electrons components = (d+, d, ne, ne+, ne+2, ne+3, ne+4, ne+5, ne+6, ne+7, ne+8, ne+9, ne+10, e, - sheath_boundary, braginskii_collisions, braginskii_thermal_force, - braginskii_friction, braginskii_heat_exchange, reactions, + sheath_boundary, braginskii_thermal_force, braginskii_collisions, + braginskii_friction, braginskii_heat_exchange, recycling, reactions, electron_force_balance, neutral_parallel_diffusion, - braginskii_ion_viscosity, braginskii_conduction, recycling) + braginskii_ion_viscosity, braginskii_conduction) normalise_metric = true # Normalise the input metric? @@ -62,15 +65,36 @@ Bnorm = 1 Tnorm = 100 [solver] -type = beuler # Backward Euler steady-state solver -max_nonlinear_iterations = 10 - diagnose = true - -atol = 1e-7 -rtol = 1e-5 - -predictor = false +type = snes +snes_type = newtonls +ksp_type = preonly # Direct solver +pc_type = lu +max_nonlinear_iterations = 16 +atol = 4e-8 +rtol = 1e-6 +stol = 0 +maxf = 2000 +maxl = 260 +matrix_free_operator = false +lag_jacobian = 1 +timestep = 0.001 + +pid_controller = true # Change timestep +target_its = 3 # Target number of nonlinear iterations +kP = 0.65 +kI = 0.3 +kD = 0.15 + +[petsc] + +# STRUMPACK as factorization backend for ILU/LU +pc_factor_mat_solver_type = strumpack + +mat_strumpack_reordering = metis # try amd/rcm if needed +mat_strumpack_colperm = true # allow column permutation for nonzero diagonal +pc_factor_shift_type = NONZERO +pc_factor_shift_amount = 1e-12 [sheath_boundary] @@ -78,7 +102,7 @@ lower_y = false upper_y = true [neutral_parallel_diffusion] - +diffusion_collisions_mode = multispecies # "afn" or "multispecies" dneut = 10 # (B / Bpol)^2 in neutral diffusion terms #################################### @@ -105,24 +129,27 @@ 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 +target_recycle_multiplier = 1 # Recycling fraction [Nd+] function = 1 -source_shape = H(mesh:y_xpt - y) * 1e20 +source_shape = H(mesh:y_xpt - y) * 1e20 # for upstream_density_feedback [Pd+] + +# Initial condition for ion pressure (in terms of hermes:Nnorm * hermes:Tnorm) function = 1 powerflux = 2.5e7 # Input power flux in W/m^2 +# Divide by length of source region to get heat source in W/m3 source = (powerflux*2/3 / (mesh:length_xpt))*H(mesh:y_xpt - y) # Input power as function of y [NVd+] -function = 0 +function = 0 # Initialise with Vi=0 #################################### @@ -136,11 +163,9 @@ AA = 2 thermal_conduction = true [Nd] - function = 0.001 [Pd] - function = 0.0001 #################################### @@ -171,8 +196,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+] @@ -193,8 +218,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+2] @@ -215,8 +240,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction [Nne+3] @@ -237,8 +262,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -261,9 +286,10 @@ 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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction + diagnose = true [Nne+5] @@ -284,8 +310,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -308,8 +334,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -332,8 +358,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -356,8 +382,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -380,8 +406,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 = true # Target recycling on +target_recycle_energy = 0 # Simple assumption target_recycle_multiplier = 0.99 # Recycling fraction diagnose = true @@ -395,20 +421,19 @@ function = 1e-7 #################################### [e] # Electrons -type = quasineutral, zero_current, evolve_pressure, noflow_boundary +type = quasineutral, evolve_pressure, zero_current, noflow_boundary noflow_upper_y = false charge = -1 AA = 1/1836 - thermal_conduction = true # in evolve_pressure +diagnose = true [Pe] function = `Pd+:function` # Same as ion pressure initially - -source = `Pd+:source` # Same as ion pressure source +source = `Pd+:source` # Same as ion pressure source #################################### @@ -421,11 +446,11 @@ type = ( d + e -> d+ + 2e, # Deuterium ionisation d+ + e -> d, # Deuterium recombination d + d+ -> d+ + d, # Charge exchange - + ne + e -> ne+ + 2e, # Neon ionisation ne+ + e -> ne, # Neon+ recombination ne+ + d -> ne + d+, # Neon+ charge exchange recombination - + ne+ + e -> ne+2 + 2e, # Neon+ ionisation ne+2 + e -> ne+, # Neon+2 recombination ne+2 + d -> ne+ + d+, # Neon+2 charge exchange recombination diff --git a/examples/tokamak-1D/extra/1D-neon/BOUT.inp b/examples/tokamak-1D/extra/1D-neon/BOUT.inp index a232a6fa6..bae7d4e04 100644 --- a/examples/tokamak-1D/extra/1D-neon/BOUT.inp +++ b/examples/tokamak-1D/extra/1D-neon/BOUT.inp @@ -16,32 +16,35 @@ nout = 100 -timestep = 10000 +timestep = 100 MXG = 0 # No guard cells in X [mesh] +# 1D simulation, use "y" as the dimension along the fieldline nx = 1 ny = 400 # Resolution along field-line nz = 1 length = 30 # Length of the domain in meters -length_xpt = 10 # Length from midplane to X-point [m] +length_xpt = 10 # Length from midplane to X-point [m] (i.e. this is where the source ends) dymin = 0.1 # Minimum grid spacing near target, as fraction of average. Must be > 0 and < 1 -# Parallel grid spacing +# Parallel grid spacing — grid refinement near the divertor target dy = (length / ny) * (1 + (1-dymin)*(1-y/pi)) -# Calculate where the source ends in grid index +# Calculate where the source ends in grid index (i.e. at the X-point) source = length_xpt / length y_xpt = pi * ( 2 - dymin - sqrt( (2-dymin)^2 - 4*(1-dymin)*source ) ) / (1 - dymin) +# Separatrix index set to -1 for 1D simulations ixseps1 = -1 ixseps2 = -1 -[restart_files] -init_missing = true # Initialise any missing variables to zero +# In 1D, the Jacobian is inversely proportional to volume +# Total flux expansion can be set up with a non-constant J profile +J = 1 [hermes] # Notes: @@ -62,16 +65,36 @@ Bnorm = 1 Tnorm = 100 [solver] -type = beuler # Backward Euler steady-state solver -snes_type = newtonls # PETSc nonlinear solver type -max_nonlinear_iterations = 10 - -diagnose = false # Output additional diagnostics? - -atol = 1e-7 -rtol = 1e-5 - -predictor = false +diagnose = true +type = snes +snes_type = newtonls +ksp_type = preonly # Direct solver +pc_type = lu +max_nonlinear_iterations = 16 +atol = 4e-8 +rtol = 1e-6 +stol = 0 +maxf = 2000 +maxl = 260 +matrix_free_operator = false +lag_jacobian = 1 +timestep = 0.001 + +pid_controller = true # Change timestep +target_its = 3 # Target number of nonlinear iterations +kP = 0.65 +kI = 0.3 +kD = 0.15 + +[petsc] + +# STRUMPACK as factorization backend for ILU/LU +pc_factor_mat_solver_type = strumpack + +mat_strumpack_reordering = metis # try amd/rcm if needed +mat_strumpack_colperm = true # allow column permutation for nonzero diagonal +pc_factor_shift_type = NONZERO +pc_factor_shift_amount = 1e-12 [sheath_boundary] @@ -79,7 +102,7 @@ lower_y = false upper_y = true [neutral_parallel_diffusion] - +diffusion_collisions_mode = multispecies # "afn" or "multispecies" dneut = 10 # (B / Bpol)^2 in neutral diffusion terms #################################### @@ -115,15 +138,18 @@ function = 1 source_shape = H(mesh:y_xpt - y) * 1e20 # for upstream_density_feedback [Pd+] + +# Initial condition for ion pressure (in terms of hermes:Nnorm * hermes:Tnorm) function = 1 powerflux = 2.5e7 # Input power flux in W/m^2 +# Divide by length of source region to get heat source in W/m3 source = (powerflux*2/3 / (mesh:length_xpt))*H(mesh:y_xpt - y) # Input power as function of y [NVd+] -function = 0 +function = 0 # Initialise with Vi=0 #################################### @@ -137,11 +163,9 @@ AA = 2 thermal_conduction = true [Nd] - function = 0.001 [Pd] - function = 0.0001 #################################### @@ -395,20 +419,19 @@ function = 1e-7 #################################### [e] # Electrons -type = quasineutral, zero_current, evolve_pressure, noflow_boundary +type = quasineutral, evolve_pressure, zero_current, noflow_boundary noflow_upper_y = false charge = -1 AA = 1/1836 - thermal_conduction = true # in evolve_pressure +diagnose = true [Pe] function = `Pd+:function` # Same as ion pressure initially - -source = `Pd+:source` # Same as ion pressure source +source = `Pd+:source` # Same as ion pressure source #################################### @@ -421,11 +444,11 @@ type = ( d + e -> d+ + 2e, # Deuterium ionisation d+ + e -> d, # Deuterium recombination d + d+ -> d+ + d, # Charge exchange - + ne + e -> ne+ + 2e, # Neon ionisation ne+ + e -> ne, # Neon+ recombination ne+ + d -> ne + d+, # Neon+ charge exchange recombination - + ne+ + e -> ne+2 + 2e, # Neon+ ionisation ne+2 + e -> ne+, # Neon+2 recombination ne+2 + d -> ne+ + d+, # Neon+2 charge exchange recombination diff --git a/examples/tokamak-1D/extra/1D-neon/README.md b/examples/tokamak-1D/extra/1D-neon/README.md index 84d4217f2..996b6f007 100644 --- a/examples/tokamak-1D/extra/1D-neon/README.md +++ b/examples/tokamak-1D/extra/1D-neon/README.md @@ -4,3 +4,11 @@ This is a more complex version of the `1D-recycling` example. It adds neon species with multiple charge states, thermal force, ionisation and recombination between states. + +*Restart from 1D-hydrogen* : Use the restart file from the 1D-hydrogen +example as the initial condition for this simulation. + +*Solver settings* : The solver settings require PETSc compiled with +STRUMPACK . Solving the system of equations with multiple species +is difficult and requires some tuning. See tips and recipes +here: https://github.com/mikekryjak/hermes-perftest/ diff --git a/include/component.hxx b/include/component.hxx index bd3dfe351..4466ed028 100644 --- a/include/component.hxx +++ b/include/component.hxx @@ -23,8 +23,8 @@ #include #include "guarded_options.hxx" -#include "permissions.hxx" #include "hermes_utils.hxx" +#include "permissions.hxx" class Solver; // Time integrator @@ -34,8 +34,9 @@ struct SpeciesInformation { SpeciesInformation(const std::vector& electrons, const std::vector& neutrals, const std::vector& positive_ions, - const std::vector & negative_ions) - : electrons(electrons), neutrals(neutrals), positive_ions(positive_ions), negative_ions(negative_ions), ions(positive_ions) { + const std::vector& negative_ions) + : electrons(electrons), neutrals(neutrals), positive_ions(positive_ions), + negative_ions(negative_ions), ions(positive_ions) { finish_construction(); } @@ -57,19 +58,20 @@ struct SpeciesInformation { } } - std::vector electrons, neutrals, positive_ions, negative_ions, ions, charged, non_electrons, all_species; - - private: - void finish_construction() { - ions = positive_ions; - ions.insert(ions.end(), negative_ions.begin(), negative_ions.end()); - charged = ions; - charged.insert(charged.end(), electrons.begin(), electrons.end()); - non_electrons = ions; - non_electrons.insert(non_electrons.end(), neutrals.begin(), neutrals.end()); - all_species = charged; - all_species.insert(all_species.end(), neutrals.begin(), neutrals.end()); - } + std::vector electrons, neutrals, positive_ions, negative_ions, ions, + charged, non_electrons, all_species; + +private: + void finish_construction() { + ions = positive_ions; + ions.insert(ions.end(), negative_ions.begin(), negative_ions.end()); + charged = ions; + charged.insert(charged.end(), electrons.begin(), electrons.end()); + non_electrons = ions; + non_electrons.insert(non_electrons.end(), neutrals.begin(), neutrals.end()); + all_species = charged; + all_species.insert(all_species.end(), neutrals.begin(), neutrals.end()); + } }; /// Interface for a component of a simulation model @@ -91,31 +93,33 @@ struct Component { /// Modify the given simulation state. This method will wrap the /// state in a GuardedOptions object and pass that to the private /// implementation of transform provided by each component. - void transform(Options &state); - + void transform(Options& state); + /// Use the final simulation state to update internal state /// (e.g. time derivatives) - virtual void finally(const Options &UNUSED(state)) { } + virtual void finally(const Options& UNUSED(state)) {} /// Add extra fields for output, or set attributes e.g docstrings - virtual void outputVars(Options &UNUSED(state)) { } + virtual void outputVars(Options& UNUSED(state)) {} /// Add extra fields to restart files - virtual void restartVars(Options &UNUSED(state)) { } + virtual void restartVars(Options& UNUSED(state)) {} /// Preconditioning - virtual void precon(const Options &UNUSED(state), BoutReal UNUSED(gamma)) { } - + virtual void precon(const Options& UNUSED(state), BoutReal UNUSED(gamma)) {} + /// Create a Component /// /// @param type The name of the component type to create (e.g. "evolve_density") /// @param name The species/name for this instance. /// @param options Component settings: options[name] are specific to this component /// @param solver Time-integration solver - static std::unique_ptr create(const std::string &type, // The type to create - const std::string &name, // The species/name for this instance - Options &options, // Component settings: options[name] are specific to this component - Solver *solver); // Time integration solver + static std::unique_ptr + create(const std::string& type, // The type to create + const std::string& name, // The species/name for this instance + Options& + options, // Component settings: options[name] are specific to this component + Solver* solver); // Time integration solver /// Tell the component the name of all species in the simulation, by type. It /// will use this information to substitute the following placeholders in @@ -140,7 +144,7 @@ struct Component { /// At the end of this function there is a call to /// Permissions::checkNoRemainingSubstitutions. All substitutions /// must be completed or else an exception will be thrown. - void declareAllSpecies(const SpeciesInformation & info); + void declareAllSpecies(const SpeciesInformation& info); protected: /// Set the level of access needed by this component for a particular variable. @@ -168,7 +172,7 @@ private: /// function. It will only allow the reading from/writing to state /// variables with the appropriate permissiosn in /// `state_variable_access`. - virtual void transform_impl(GuardedOptions &state) = 0; + virtual void transform_impl(GuardedOptions& state) = 0; }; /////////////////////////////////////////////////////////////////// @@ -205,33 +209,32 @@ using RegisterComponent = ComponentFactory::RegisterInFactory; /// @tparam T The type the option should be converted to /// /// @param option The Option whose value will be returned -template +template T getNonFinal(const Options& option) { if (!option.isSet()) { throw BoutException("Option {:s} has no value", option.str()); } try { return bout::utils::variantStaticCastOrThrow(option.value); - } catch (const std::bad_cast &e) { + } catch (const std::bad_cast& e) { // Convert to a more useful error message - throw BoutException("Could not convert {:s} to type {:s}", - option.str(), typeid(T).name()); + throw BoutException("Could not convert {:s} to type {:s}", option.str(), + typeid(T).name()); } } -template -T getNonFinal(const GuardedOptions & option) { +template +T getNonFinal(const GuardedOptions& option) { return getNonFinal(option.get()); } #define TOSTRING_(x) #x #define TOSTRING(x) TOSTRING_(x) - namespace hermes { /// Enable a function if and only if `T` is a (subclass of) `GuardedOptions` template using EnableIfGuardedOption = std::enable_if_t>; -} +} // namespace hermes /// Faster non-printing getter for Options /// If this fails, it will throw BoutException @@ -243,7 +246,7 @@ using EnableIfGuardedOption = std::enable_if_t +template T get(const Options& option, [[maybe_unused]] const std::string& location = "") { #if CHECKLEVEL >= 1 // Mark option as final, both inside the domain and the boundary @@ -252,8 +255,8 @@ T get(const Options& option, [[maybe_unused]] const std::string& location = "") #endif return getNonFinal(option); } -template -T get(const GuardedOptions & option, const std::string& location = "") { +template +T get(const GuardedOptions& option, const std::string& location = "") { return get(option.get(), location); } @@ -261,25 +264,23 @@ T get(const GuardedOptions & option, const std::string& location = "") { /// Sets the final flag so setting the value /// afterwards will lead to an error bool isSetFinal(const Options& option, const std::string& location = ""); -bool isSetFinal(const GuardedOptions & option, const std::string& location = ""); +bool isSetFinal(const GuardedOptions& option, const std::string& location = ""); #if CHECKLEVEL >= 1 /// A wrapper around isSetFinal() which captures debugging information /// /// Usage: /// if (IS_SET(option["value"])); -#define IS_SET(option) \ - isSetFinal(option, __FILE__ ":" TOSTRING(__LINE__)) +#define IS_SET(option) isSetFinal(option, __FILE__ ":" TOSTRING(__LINE__)) #else -#define IS_SET(option) \ - isSetFinal(option) +#define IS_SET(option) isSetFinal(option) #endif /// Check if an option can be fetched /// Sets the final flag so setting the value in the domain /// afterwards will lead to an error bool isSetFinalNoBoundary(const Options& option, const std::string& location = ""); -bool isSetFinalNoBoundary(const GuardedOptions & option, const std::string& location = ""); +bool isSetFinalNoBoundary(const GuardedOptions& option, const std::string& location = ""); #if CHECKLEVEL >= 1 /// A wrapper around isSetFinalNoBoundary() which captures debugging information @@ -289,8 +290,7 @@ bool isSetFinalNoBoundary(const GuardedOptions & option, const std::string& loca #define IS_SET_NOBOUNDARY(option) \ isSetFinalNoBoundary(option, __FILE__ ":" TOSTRING(__LINE__)) #else -#define IS_SET_NOBOUNDARY(option) \ - isSetFinalNoBoundary(option) +#define IS_SET_NOBOUNDARY(option) isSetFinalNoBoundary(option) #endif #if CHECKLEVEL >= 1 @@ -298,11 +298,9 @@ bool isSetFinalNoBoundary(const GuardedOptions & option, const std::string& loca /// /// Usage: /// auto var = GET_VALUE(Field3D, option["value"]); -#define GET_VALUE(Type, option) \ - get(option, __FILE__ ":" TOSTRING(__LINE__)) +#define GET_VALUE(Type, option) get(option, __FILE__ ":" TOSTRING(__LINE__)) #else -#define GET_VALUE(Type, option) \ - get(option) +#define GET_VALUE(Type, option) get(option) #endif /// Faster non-printing getter for Options @@ -317,8 +315,9 @@ bool isSetFinalNoBoundary(const GuardedOptions & option, const std::string& loca /// /// @param option The Option whose value will be returned /// @param location An optional string to indicate where this value is used -template -T getNoBoundary(const Options& option, [[maybe_unused]] const std::string& location = "") { +template +T getNoBoundary(const Options& option, + [[maybe_unused]] const std::string& location = "") { #if CHECKLEVEL >= 1 // Mark option as final inside the domain const_cast(option).attributes["final-domain"] = location; @@ -326,7 +325,7 @@ T getNoBoundary(const Options& option, [[maybe_unused]] const std::string& locat return getNonFinal(option); } -template> +template > T getNoBoundary(GO&& option, const std::string& location = "") { return getNoBoundary(std::forward(option).get(Regions::Interior), location); } @@ -339,8 +338,7 @@ T getNoBoundary(GO&& option, const std::string& location = "") { #define GET_NOBOUNDARY(Type, option) \ getNoBoundary(option, __FILE__ ":" TOSTRING(__LINE__)) #else -#define GET_NOBOUNDARY(Type, option) \ - getNoBoundary(option) +#define GET_NOBOUNDARY(Type, option) getNoBoundary(option) #endif /// Check whether value is valid, returning true @@ -352,9 +350,9 @@ bool hermesDataInvalid([[maybe_unused]] const T& value) { /// Check Field3D values. /// Doesn't check boundary cells -template<> +template <> inline bool hermesDataInvalid(const Field3D& value) { - for (auto& i : value.getRegion("RGN_NOBNDRY")) { + for (const auto& i : value.getRegion("RGN_NOBNDRY")) { if (!std::isfinite(value[i])) { return true; } @@ -364,7 +362,7 @@ inline bool hermesDataInvalid(const Field3D& value) { /// Check Field2D values. /// Doesn't check boundary cells -template<> +template <> inline bool hermesDataInvalid(const Field2D& value) { for (auto& i : value.getRegion("RGN_NOBNDRY")) { if (!std::isfinite(value[i])) { @@ -381,7 +379,7 @@ inline bool hermesDataInvalid(const Field2D& value) { /// This is to prevent values being modified after use. /// /// @tparam T The type of the value to set. Usually this is inferred -template +template Options& set(Options& option, T value) { // Check that the value has not already been used #if CHECKLEVEL >= 1 @@ -404,7 +402,7 @@ Options& set(Options& option, T value) { return option; } -template> +template > decltype(auto) set(GO&& option, T value) { set(std::forward(option).getWritable(), value); return std::forward(option); @@ -418,7 +416,7 @@ decltype(auto) set(GO&& option, T value) { /// or getNonFinal. /// /// @tparam T The type of the value to set. Usually this is inferred -template +template Options& setBoundary(Options& option, T value) { // Check that the value has not already been used #if CHECKLEVEL >= 1 @@ -431,7 +429,7 @@ Options& setBoundary(Options& option, T value) { return option; } -template> +template > decltype(auto) setBoundary(GO&& option, T value) { setBoundary(std::forward(option).getWritable(Regions::Boundaries), value); return std::forward(option); @@ -445,22 +443,25 @@ decltype(auto) setBoundary(GO&& option, T value) { /// /// @param option The value to modify (or set if not already set) /// @param value The quantity to add. -template +template Options& add(Options& option, T value) { if (!option.isSet()) { return set(option, value); } else { try { - return set(option, value + bout::utils::variantStaticCastOrThrow(option.value)); - } catch (const std::bad_cast &e) { + return set(option, + value + + bout::utils::variantStaticCastOrThrow( + option.value)); + } catch (const std::bad_cast& e) { // Convert to a more useful error message - throw BoutException("Could not convert {:s} to type {:s}", - option.str(), typeid(T).name()); + throw BoutException("Could not convert {:s} to type {:s}", option.str(), + typeid(T).name()); } } } -template> +template > decltype(auto) add(GO&& option, T value) { add(std::forward(option).getWritable(), value); return std::forward(option); @@ -471,50 +472,61 @@ decltype(auto) add(GO&& option, T value) { /// /// @param option The value to modify (or set if not already set) /// @param value The quantity to add. -template +template Options& subtract(Options& option, T value) { if (!option.isSet()) { return set(option, -value); } else { try { - return set(option, bout::utils::variantStaticCastOrThrow(option.value) - value); - } catch (const std::bad_cast &e) { + return set(option, bout::utils::variantStaticCastOrThrow( + option.value) + - value); + } catch (const std::bad_cast& e) { // Convert to a more useful error message - throw BoutException("Could not convert {:s} to type {:s}", - option.str(), typeid(T).name()); + throw BoutException("Could not convert {:s} to type {:s}", option.str(), + typeid(T).name()); } } } -template> +template > decltype(auto) subtract(GO&& option, T value) { subtract(std::forward(option).getWritable(), value); return std::forward(option); } -template -void set_with_attrs(Options& option, T value, std::initializer_list> attrs) { +template +void set_with_attrs( + Options& option, T value, + std::initializer_list> attrs) { option.force(value); option.setAttributes(attrs); } -template> -void set_with_attrs(GO&& option, T value, std::initializer_list> attrs) { +template > +void set_with_attrs( + GO&& option, T value, + std::initializer_list> attrs) { set_with_attrs(std::forward(option).getWritable(), value, attrs); } #if CHECKLEVEL >= 1 -template<> -inline void set_with_attrs(Options& option, Field3D value, std::initializer_list> attrs) { +template <> +inline void set_with_attrs( + Options& option, Field3D value, + std::initializer_list> attrs) { if (!value.isAllocated()) { - throw BoutException("set_with_attrs: Field3D assigned to {:s} is not allocated", option.str()); + throw BoutException("set_with_attrs: Field3D assigned to {:s} is not allocated", + option.str()); } option.force(value); option.setAttributes(attrs); } -template> -inline void set_with_attrs(GO&& option, Field3D value, std::initializer_list> attrs) { +template > +inline void set_with_attrs( + GO&& option, Field3D value, + std::initializer_list> attrs) { set_with_attrs(std::forward(option).getWritable(), std::move(value), attrs); } #endif diff --git a/include/neutral_parallel_diffusion.hxx b/include/neutral_parallel_diffusion.hxx index 9b93f85eb..5e3702613 100644 --- a/include/neutral_parallel_diffusion.hxx +++ b/include/neutral_parallel_diffusion.hxx @@ -3,6 +3,15 @@ #define NEUTRAL_PARALLEL_DIFFUSION_H #include "component.hxx" +#include "guarded_options.hxx" +#include "permissions.hxx" +#include +#include +#include + +#include +#include +#include /// Add effective diffusion of neutrals in a 1D system, /// by projecting cross-field diffusion onto parallel distance. @@ -34,24 +43,27 @@ struct NeutralParallelDiffusion : public Component { .as(); diagnose = options["diagnose"] - .doc("Output additional diagnostics?") - .withDefault(false); + .doc("Output additional diagnostics?") + .withDefault(false); + + diffusion_collisions_mode = + options["diffusion_collisions_mode"] + .doc("Can be afn: CX and IZ, or multispecies: CX and nn, ni, ne (if enabled)") + .withDefault("afn"); - diffusion_collisions_mode = options["diffusion_collisions_mode"] - .doc("Can be afn: CX and IZ, or multispecies: CX and nn, ni, ne (if enabled)") - .withDefault("afn"); - equation_fix = options["equation_fix"] - .doc("Fix correcting pressure advection and conductivity factors?") - .withDefault(true); + .doc("Fix correcting pressure advection and conductivity factors?") + .withDefault(true); - perpendicular_conduction = options["perpendicular_conduction"] - .doc("Enable parallel projection of perpendicular conduction?") - .withDefault(true); + perpendicular_conduction = + options["perpendicular_conduction"] + .doc("Enable parallel projection of perpendicular conduction?") + .withDefault(true); - perpendicular_viscosity = options["perpendicular_viscosity"] - .doc("Enable parallel projection of perpendicular viscosity?") - .withDefault(true); + perpendicular_viscosity = + options["perpendicular_viscosity"] + .doc("Enable parallel projection of perpendicular viscosity?") + .withDefault(true); // FIXME: strictly speaking, momentum is not optional if velocity has been set substitutePermissions("optional_inputs", {"pressure", "velocity", "momentum"}); @@ -63,24 +75,27 @@ struct NeutralParallelDiffusion : public Component { } /// Save variables to the output - void outputVars(Options &state) override; + void outputVars(Options& state) override; + private: BoutReal dneut; ///< cross-field diffusion projection (B / Bpol)^2 bool diagnose; ///< Output diagnostics? - std::map> collision_names; ///< Collisions used for collisionality - std::string diffusion_collisions_mode; ///< Collision selection, either afn or multispecies - Field3D nu; ///< Collision frequency for conduction - bool equation_fix; ///< Fix incorrect 3/2 factor in pressure advection? + std::map> + collision_names; ///< Collisions used for collisionality + std::string + diffusion_collisions_mode; ///< Collision selection, either afn or multispecies + Field3D nu; ///< Collision frequency for conduction + bool equation_fix; ///< Fix incorrect 3/2 factor in pressure advection? bool perpendicular_conduction; ///< Enable conduction? - bool perpendicular_viscosity; ///< Enable viscosity? + bool perpendicular_viscosity; ///< Enable viscosity? /// Per-species diagnostics struct Diagnostics { Field3D Dn; ///< Diffusion coefficient - Field3D S; ///< Particle source - Field3D E; ///< Energy source - Field3D F; ///< Momentum source + Field3D S; ///< Particle source + Field3D E; ///< Energy source + Field3D F; ///< Momentum source }; /// Store diagnostics for each species diff --git a/src/braginskii_collisions.cxx b/src/braginskii_collisions.cxx index d66716e69..edf87ce00 100644 --- a/src/braginskii_collisions.cxx +++ b/src/braginskii_collisions.cxx @@ -18,7 +18,9 @@ #include "../include/braginskii_collisions.hxx" #include "../include/component.hxx" +#include "../include/guarded_options.hxx" #include "../include/hermes_utils.hxx" +#include "../include/permissions.hxx" BraginskiiCollisions::BraginskiiCollisions(const std::string& name, Options& alloptions, Solver*) @@ -148,7 +150,7 @@ void BraginskiiCollisions::collide(GuardedOptions& species1, GuardedOptions& spe }); add(species2["collision_frequency"], nu); // Total collision frequency - std::string const coll_name = + const std::string coll_name = species2.name() + std::string("_") + species1.name() + std::string("_coll"); set(species2["collision_frequencies"][coll_name], nu); // Collision frequency for individual reaction @@ -359,7 +361,7 @@ void BraginskiiCollisions::transform_impl(GuardedOptions& state) { const BoutReal charge2 = Z2 * SI::qe; // in Coulombs // Ion-ion collisions - Field3D const nu_12 = filledFrom(density1, [&](auto& i) { + const Field3D nu_12 = filledFrom(density1, [&](auto& i) { const BoutReal Tlim1 = softFloor(temperature1[i], 0.1); const BoutReal Tlim2 = softFloor(temperature2[i], 0.1); @@ -367,10 +369,10 @@ void BraginskiiCollisions::transform_impl(GuardedOptions& state) { const BoutReal Nlim2 = softFloor(density2[i], 1e10); // Coulomb logarithm - BoutReal const coulomb_log = + const BoutReal coulomb_log = 29.91 - log((Z1 * Z2 * (AA1 + AA2)) / (AA1 * Tlim2 + AA2 * Tlim1) - * sqrt(Nlim1 * SQ(Z1) / Tlim1 + Nlim2 * SQ(Z2) / Tlim2)); + * sqrt((Nlim1 * SQ(Z1) / Tlim1) + (Nlim2 * SQ(Z2) / Tlim2))); // Calculate v_a^2, v_b^2 const BoutReal v1sq = 2 * Tlim1 * SI::qe / mass1; @@ -390,9 +392,13 @@ void BraginskiiCollisions::transform_impl(GuardedOptions& state) { } else { // species1 charged, species2 neutral + if (!ion_neutral) { + continue; + } + // Scattering of charged species 1 // Neutral density - Field3D const Nn = GET_NOBOUNDARY(Field3D, species2["density"]); + const Field3D Nn = GET_NOBOUNDARY(Field3D, species2["density"]); BoutReal a0 = 5e-19; // Cross-section [m^2] @@ -509,8 +515,8 @@ void BraginskiiCollisions::outputVars(Options& state) { const std::map& level2 = section.getChildren(); for (auto s2 = std::begin(level2); s2 != std::end(level2); ++s2) { std::string A = s1->first; - std::string const B = s2->first; - std::string const AB = A + B; + const std::string B = s2->first; + const std::string AB = A + B; // Collision frequencies set_with_attrs(state[std::string("K") + AB + std::string("_coll")], diff --git a/src/neutral_parallel_diffusion.cxx b/src/neutral_parallel_diffusion.cxx index 8f9bf83c5..0eee96d90 100644 --- a/src/neutral_parallel_diffusion.cxx +++ b/src/neutral_parallel_diffusion.cxx @@ -1,8 +1,18 @@ #include "../include/neutral_parallel_diffusion.hxx" +#include "../include/component.hxx" +#include "../include/guarded_options.hxx" #include "../include/hermes_utils.hxx" +#include +#include #include +#include #include +#include +#include +#include + +#include using bout::globals::mesh; @@ -20,68 +30,71 @@ void NeutralParallelDiffusion::transform_impl(GuardedOptions& state) { } // Initialise collision_names for each species if not already done if (collision_names.find(species_name) == collision_names.end()) { - collision_names[species_name] = {}; // Initialise with empty vector + collision_names[species_name] = {}; // Initialise with empty vector } const BoutReal AA = GET_VALUE(BoutReal, species["AA"]); // Atomic mass const Field3D Nn = GET_VALUE(Field3D, species["density"]); const Field3D Tn = GET_VALUE(Field3D, species["temperature"]); - const Field3D Pn = IS_SET(species["pressure"]) ? - GET_VALUE(Field3D, species["pressure"]) : Nn * Tn; + const Field3D Pn = + IS_SET(species["pressure"]) ? GET_VALUE(Field3D, species["pressure"]) : Nn * Tn; // Collisionality // Multispecies mode: in, en, nn, cx // New mode: cx, iz (in line with SOLPS AFN, Horsten 2017) - if (collision_names[species_name].empty()) { // Calculate only once - at the beginning + if (collision_names[species_name].empty()) { // Calculate only once - at the beginning if (diffusion_collisions_mode == "afn") { for (const auto& collision : species["collision_frequencies"].getChildren()) { - std::string collision_name = collision.first; + const std::string collision_name = collision.first; - if (// Charge exchange - (collisionSpeciesMatch( - collision_name, species_name, "+", "cx", "partial")) or + if ( // Charge exchange + (collisionSpeciesMatch(collision_name, species_name, "+", "cx", "partial")) + or // Ionisation - (collisionSpeciesMatch( - collision_name, species_name, "+", "iz", "partial"))) { - - collision_names[species_name].push_back(collision_name); - } + (collisionSpeciesMatch(collision_name, species_name, "+", "iz", + "partial"))) { + + collision_names[species_name].push_back(collision_name); + } } - // Multispecies mode: all collisions and CX are included + // Multispecies mode: all collisions and CX are included } else if (diffusion_collisions_mode == "multispecies") { for (const auto& collision : species["collision_frequencies"].getChildren()) { - std::string collision_name = collision.first; + const std::string collision_name = collision.first; - if (// Charge exchange - (collisionSpeciesMatch( - collision_name, species_name, "", "cx", "partial")) or + if ( // Charge exchange + (collisionSpeciesMatch(collision_name, species_name, "", "cx", "partial")) + or // Any collision (ne, ni, nn) - (collisionSpeciesMatch( - collision_name, species_name, "", "coll", "partial"))) { - - collision_names[species_name].push_back(collision_name); - } + (collisionSpeciesMatch(collision_name, species_name, "", "coll", + "partial"))) { + + collision_names[species_name].push_back(collision_name); + } } } else { - throw BoutException("\tdiffusion_collisions_mode for {:s} must be either multispecies or afn", species_name); + throw BoutException( + "\tdiffusion_collisions_mode for {:s} must be either multispecies or afn", + species_name); } if (collision_names[species_name].empty()) { - throw BoutException("\tNo collisions found for {:s} in neutral_parallel_diffusion for selected collisions mode", species_name); + throw BoutException("\tNo collisions found for {:s} in " + "neutral_parallel_diffusion for selected collisions mode", + species_name); } // Write chosen collisions to log file output_info.write("\t{:s} neutral diffusion collisionality mode: '{:s}' using ", - species_name, diffusion_collisions_mode); - for (const auto& collision : collision_names[species_name]) { + species_name, diffusion_collisions_mode); + for (const auto& collision : collision_names[species_name]) { output_info.write("{:s} ", collision); } output_info.write("\n"); - } // Collect the collisionalities based on list of names @@ -89,13 +102,15 @@ void NeutralParallelDiffusion::transform_impl(GuardedOptions& state) { for (const auto& collision_name : collision_names[species_name]) { nu += GET_VALUE(Field3D, species["collision_frequencies"][collision_name]); } + // If nu is zero then Dn will be NaN + nu = floor(nu, 1e-5); // Diffusion coefficient BoutReal advection_factor = 0; BoutReal kappa_factor = 0; if (equation_fix) { - advection_factor = (5. / 2); // This is equivalent to 5/3 if on pressure basis + advection_factor = (5. / 2); // This is equivalent to 5/3 if on pressure basis kappa_factor = (5. / 2); } else { advection_factor = (3. / 2); @@ -108,22 +123,22 @@ void NeutralParallelDiffusion::transform_impl(GuardedOptions& state) { mesh->communicate(Dn); // Cross-field diffusion calculated from pressure gradient - // This is the pressure-diffusion approximation + // This is the pressure-diffusion approximation Field3D logPn = log(floor(Pn, 1e-7)); logPn.applyBoundary("neumann"); // Particle advection - Field3D S = FV::Div_par_K_Grad_par(Dn * Nn, logPn); + const Field3D S = FV::Div_par_K_Grad_par(Dn * Nn, logPn); add(species["density_source"], S); Field3D kappa_n = kappa_factor * Nn * Dn; kappa_n.applyBoundary("neumann"); // Heat transfer - Field3D E = + FV::Div_par_K_Grad_par( - Dn * advection_factor * Pn, logPn); // Pressure advection + Field3D E = + +FV::Div_par_K_Grad_par(Dn * advection_factor * Pn, logPn); // Pressure advection if (perpendicular_conduction) { - E += FV::Div_par_K_Grad_par(kappa_n, Tn); // Conduction + E += FV::Div_par_K_Grad_par(kappa_n, Tn); // Conduction } add(species["energy_source"], E); @@ -135,10 +150,10 @@ void NeutralParallelDiffusion::transform_impl(GuardedOptions& state) { // Transport Processes in Gases", 1972 // - Field3D Vn = GET_VALUE(Field3D, species["velocity"]); - Field3D NVn = GET_VALUE(Field3D, species["momentum"]); + const Field3D Vn = GET_VALUE(Field3D, species["velocity"]); + const Field3D NVn = GET_VALUE(Field3D, species["momentum"]); - Field3D eta_n = (2. / 5) * kappa_n; + const Field3D eta_n = (2. / 5) * kappa_n; // Momentum diffusion F = FV::Div_par_K_Grad_par(NVn * Dn, logPn) + FV::Div_par_K_Grad_par(eta_n, Vn); @@ -150,7 +165,7 @@ void NeutralParallelDiffusion::transform_impl(GuardedOptions& state) { auto search = diagnostics.find(species_name); if (search == diagnostics.end()) { // First time, create diagnostic - diagnostics.emplace(species_name, Diagnostics {Dn, S, E, F}); + diagnostics.emplace(species_name, Diagnostics{Dn, S, E, F}); } else { // Update diagnostic values auto& d = search->second; @@ -177,40 +192,40 @@ void NeutralParallelDiffusion::outputVars(Options& state) { const std::string& species_name = it.first; const auto& d = it.second; set_with_attrs(state[std::string("D") + species_name + std::string("_Dpar")], d.Dn, - {{"time_dimension", "t"}, - {"units", "m^2/s"}, - {"conversion", rho_s0 * Cs0}, - {"standard_name", "diffusion coefficient"}, - {"long_name", species_name + " particle diffusion coefficient"}, - {"species", species_name}, - {"source", "neutral_parallel_diffusion"}}); + {{"time_dimension", "t"}, + {"units", "m^2/s"}, + {"conversion", rho_s0 * Cs0}, + {"standard_name", "diffusion coefficient"}, + {"long_name", species_name + " particle diffusion coefficient"}, + {"species", species_name}, + {"source", "neutral_parallel_diffusion"}}); set_with_attrs(state[std::string("S") + species_name + std::string("_Dpar")], d.S, - {{"time_dimension", "t"}, - {"units", "s^-1"}, - {"conversion", Nnorm * Omega_ci}, - {"standard_name", "particle diffusion"}, - {"long_name", species_name + " particle source due to diffusion"}, - {"species", species_name}, - {"source", "neutral_parallel_diffusion"}}); + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "particle diffusion"}, + {"long_name", species_name + " particle source due to diffusion"}, + {"species", species_name}, + {"source", "neutral_parallel_diffusion"}}); set_with_attrs(state[std::string("E") + species_name + std::string("_Dpar")], d.E, - {{"time_dimension", "t"}, - {"units", "W / m^3"}, - {"conversion", SI::qe * Tnorm * Nnorm * Omega_ci}, - {"standard_name", "energy diffusion"}, - {"long_name", species_name + " energy source due to diffusion"}, - {"species", species_name}, - {"source", "neutral_parallel_diffusion"}}); + {{"time_dimension", "t"}, + {"units", "W / m^3"}, + {"conversion", SI::qe * Tnorm * Nnorm * Omega_ci}, + {"standard_name", "energy diffusion"}, + {"long_name", species_name + " energy source due to diffusion"}, + {"species", species_name}, + {"source", "neutral_parallel_diffusion"}}); set_with_attrs(state[std::string("F") + species_name + std::string("_Dpar")], d.F, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"standard_name", "momentum diffusion"}, - {"long_name", species_name + " momentum source due to diffusion"}, - {"species", species_name}, - {"source", "neutral_parallel_diffusion"}}); + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum diffusion"}, + {"long_name", species_name + " momentum source due to diffusion"}, + {"species", species_name}, + {"source", "neutral_parallel_diffusion"}}); } } } diff --git a/tests/unit/test_braginskii_collisions.cxx b/tests/unit/test_braginskii_collisions.cxx index 1ad9b0a99..de2e4818b 100644 --- a/tests/unit/test_braginskii_collisions.cxx +++ b/tests/unit/test_braginskii_collisions.cxx @@ -29,7 +29,7 @@ TEST_F(BraginskiiCollisionsTest, CreateComponent) { options["units"]["meters"] = 1.0; options["units"]["seconds"] = 1.0; options["units"]["inv_meters_cubed"] = 1e19; - BraginskiiCollisions const component("test", options, nullptr); + const BraginskiiCollisions component("test", options, nullptr); } TEST_F(BraginskiiCollisionsTest, OnlyElectrons) { @@ -167,3 +167,69 @@ TEST_F(BraginskiiCollisionsTest, TnormDependence) { ASSERT_DOUBLE_EQ(get(state["species"]["d+"]["collision_frequency"])(0, 0, 0), get(state2["species"]["d+"]["collision_frequency"])(0, 0, 0)); } + +TEST_F(BraginskiiCollisionsTest, IonNeutral) { + Options options{ + {"units", + {{"eV", 1.0}, {"meters", 1.0}, {"seconds", 1.0}, {"inv_meters_cubed", 1.0}}}, + {"test", + {{"electron_neutral", false}, + {"ion_neutral", true}, + {"electron_electron", false}, + {"ion_ion", false}, + {"neutral_neutral", false}}}}; + + BraginskiiCollisions component("test", options, nullptr); + + // Put ion in between neutral species to test both orders (ion-neutral and neutral-ion) + Options state{ + {"species", + {{"d", {{"density", 1e18}, {"temperature", 3}, {"AA", 2}}}, + {"d+", {{"density", 2e19}, {"temperature", 20}, {"charge", 1}, {"AA", 2}}}, + {"d2", {{"density", 1e17}, {"temperature", 1}, {"AA", 4}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d", "d2"}, {"d+"}, {})); + component.transform(state); + + ASSERT_TRUE(state["species"]["d"].isSet("collision_frequency")); + ASSERT_TRUE(state["species"]["d"]["collision_frequencies"].isSet("d_d+_coll")); + + ASSERT_TRUE(state["species"]["d2"].isSet("collision_frequency")); + ASSERT_TRUE(state["species"]["d2"]["collision_frequencies"].isSet("d2_d+_coll")); + + ASSERT_TRUE(state["species"]["d+"].isSet("collision_frequency")); + ASSERT_TRUE(state["species"]["d+"]["collision_frequencies"].isSet("d+_d_coll")); + ASSERT_TRUE(state["species"]["d+"]["collision_frequencies"].isSet("d+_d2_coll")); + + Options options2{ + {"units", + {{"eV", 1.0}, {"meters", 1.0}, {"seconds", 1.0}, {"inv_meters_cubed", 1.0}}}, + {"test", + {{"electron_neutral", false}, + {"ion_neutral", false}, // Turn off ion_neutral collisions + {"electron_electron", false}, + {"ion_ion", false}, + {"neutral_neutral", false}}}}; + + BraginskiiCollisions component2("test", options2, nullptr); + + Options state2{ + {"species", + {{"d", {{"density", 1e18}, {"temperature", 3}, {"AA", 2}}}, + {"d+", {{"density", 2e19}, {"temperature", 20}, {"charge", 1}, {"AA", 2}}}, + {"d2", {{"density", 1e17}, {"temperature", 1}, {"AA", 4}}}}}}; + + component2.declareAllSpecies(SpeciesInformation({}, {"d", "d2"}, {"d+"}, {})); + component2.transform(state2); + + // All ion-neutral collision frequencies should be missing + ASSERT_FALSE(state2["species"]["d"].isSet("collision_frequency")); + ASSERT_FALSE(state2["species"]["d"]["collision_frequencies"].isSet("d_d+_coll")); + + ASSERT_FALSE(state2["species"]["d2"].isSet("collision_frequency")); + ASSERT_FALSE(state2["species"]["d2"]["collision_frequencies"].isSet("d2_d+_coll")); + + ASSERT_FALSE(state2["species"]["d+"].isSet("collision_frequency")); + ASSERT_FALSE(state2["species"]["d+"]["collision_frequencies"].isSet("d+_d_coll")); + ASSERT_FALSE(state2["species"]["d+"]["collision_frequencies"].isSet("d+_d2_coll")); +} diff --git a/tests/unit/test_neutral_parallel_diffusion.cxx b/tests/unit/test_neutral_parallel_diffusion.cxx new file mode 100644 index 000000000..850885895 --- /dev/null +++ b/tests/unit/test_neutral_parallel_diffusion.cxx @@ -0,0 +1,218 @@ +#include + +#include +#include +#include +#include +#include +#include "gtest/gtest.h" + +#include "../../include/component.hxx" +#include "../../include/neutral_parallel_diffusion.hxx" +#include "fake_mesh_fixture.hxx" + +#include + +/// Global mesh +namespace bout { +namespace globals { +extern Mesh* mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; + +// Reuse the "standard" fixture for FakeMesh +using NeutralParallelDiffusionTest = FakeMeshFixture; + +TEST_F(NeutralParallelDiffusionTest, ComponentRequiresDneut) { + Options options; // No dneut setting + EXPECT_THROW(const NeutralParallelDiffusion component("test", options, nullptr), + BoutException); +} + +TEST_F(NeutralParallelDiffusionTest, CreateComponent) { + Options options{{"test", {{"dneut", 1.0}}}}; + const NeutralParallelDiffusion component("test", options, nullptr); +} + +TEST_F(NeutralParallelDiffusionTest, NoNeutrals) { + Options options{{"test", {{"dneut", 1.0}}}}; + NeutralParallelDiffusion component("test", options, nullptr); + + // State without any neutrals + Options state{{"species", + {{"e", {{"AA", 1. / 1836}, {"charge", -1.0}}}, + {"d+", {{"AA", 2.0}, {"charge", 1.0}}}, + {"d2-", {{"AA", 2.0}, {"charge", -1.0}}}}}}; + + component.declareAllSpecies(SpeciesInformation({"e"}, {}, {"d+"}, {"d2-"})); + component.transform(state); + + ASSERT_FALSE(state["species"]["e"].isSet("density_source")); + ASSERT_FALSE(state["species"]["d+"].isSet("density_source")); + ASSERT_FALSE(state["species"]["d2-"].isSet("density_source")); +} + +TEST_F(NeutralParallelDiffusionTest, NeutralNoCollisions) { + Options options{{"test", {{"dneut", 1.0}}}}; + NeutralParallelDiffusion component("test", options, nullptr); + + // State without any collision frequencies + Options state{ + {"species", {{"d", {{"AA", 2.0}, {"density", 1e19}, {"temperature", 1.0}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d"}, {}, {})); + EXPECT_THROW(component.transform(state), BoutException); +} + +TEST_F(NeutralParallelDiffusionTest, NeutralNoVelocity) { + Options options{{"test", {{"dneut", 1.0}}}}; + NeutralParallelDiffusion component("test", options, nullptr); + + // State without any collision frequencies + Options state{{"species", + {{"d", + {{"AA", 2.0}, + {"density", 1e19}, + {"temperature", 1.0}, + {"collision_frequencies", {{"d_d+_iz", 1.0}}}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d"}, {}, {})); + component.transform(state); + + ASSERT_TRUE(state["species"]["d"].isSet("density_source")); + ASSERT_TRUE(state["species"]["d"].isSet("energy_source")); + ASSERT_FALSE(state["species"]["d"].isSet("momentum_source")); +} + +TEST_F(NeutralParallelDiffusionTest, NeutralWithVelocity) { + Options options{{"test", {{"dneut", 1.0}}}}; + NeutralParallelDiffusion component("test", options, nullptr); + + // State without any collision frequencies + Options state{{"species", + {{"d", + {{"AA", 2.0}, + {"density", 1e19}, + {"temperature", 1.0}, + {"velocity", 0.0}, + {"momentum", 0.0}, + {"collision_frequencies", {{"d_d+_iz", 1.0}}}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d"}, {}, {})); + component.transform(state); + + ASSERT_TRUE(state["species"]["d"].isSet("density_source")); + ASSERT_TRUE(state["species"]["d"].isSet("energy_source")); + ASSERT_TRUE(state["species"]["d"].isSet("momentum_source")); +} + +TEST_F(NeutralParallelDiffusionTest, NeutralZeroCollisionFrequency) { + Options options{{"test", {{"dneut", 1.0}}}}; + NeutralParallelDiffusion component("test", options, nullptr); + + // State without any collision frequencies + Options state{{"species", + {{"d", + {{"AA", 2.0}, + {"density", 1e19}, + {"temperature", 1.0}, + {"velocity", 0.0}, + {"momentum", 0.0}, + {"collision_frequencies", {{"d_d+_iz", 0.0}, {"d_d+_cx", 0.0}}}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d"}, {}, {})); + component.transform(state); + + ASSERT_TRUE(state["species"]["d"].isSet("density_source")); + ASSERT_TRUE(state["species"]["d"].isSet("energy_source")); + ASSERT_TRUE(state["species"]["d"].isSet("momentum_source")); + + // Check that density source is finite + const Field3D density_source = get(state["species"]["d"]["density_source"]); + const auto& region = density_source.getRegion("RGN_NOBNDRY"); + ASSERT_TRUE(std::all_of(region.begin(), region.end(), [&](const auto& i) { + return std::isfinite(density_source[i]); + })); +} + +TEST_F(NeutralParallelDiffusionTest, InvalidCollisionsMode) { + Options options{ + {"test", + {{"dneut", 1.0}, {"diffusion_collisions_mode", "multiecies"}}}}; // Mis-spelled + + NeutralParallelDiffusion component("test", options, nullptr); + // State without any collision frequencies + Options state{{"species", + {{"d", + {{"AA", 2.0}, + {"density", 1e19}, + {"temperature", 1.0}, + {"velocity", 0.0}, + {"momentum", 0.0}, + {"collision_frequencies", {{"d_d+_iz", 0.0}, {"d_d+_cx", 0.0}}}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d"}, {}, {})); + + EXPECT_THROW(component.transform(state);, BoutException); +} + +TEST_F(NeutralParallelDiffusionTest, Multispecies) { + Options options{ + {"test", {{"dneut", 1.0}, {"diffusion_collisions_mode", "multispecies"}}}}; + NeutralParallelDiffusion component("test", options, nullptr); + + // State without any collision frequencies + Options state{{"species", + {{"d", + {{"AA", 2.0}, + {"density", 1e19}, + {"temperature", 1.0}, + {"velocity", 0.0}, + {"momentum", 0.0}, + {"collision_frequencies", {{"d_d+_iz", 0.0}, {"d_d+_cx", 0.0}}}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d"}, {}, {})); + component.transform(state); + + ASSERT_TRUE(state["species"]["d"].isSet("density_source")); + ASSERT_TRUE(state["species"]["d"].isSet("energy_source")); + ASSERT_TRUE(state["species"]["d"].isSet("momentum_source")); + + // Check that density source is finite + const Field3D density_source = get(state["species"]["d"]["density_source"]); + const auto& region = density_source.getRegion("RGN_NOBNDRY"); + ASSERT_TRUE(std::all_of(region.begin(), region.end(), [&](const auto& i) { + return std::isfinite(density_source[i]); + })); +} + +TEST_F(NeutralParallelDiffusionTest, Diagnose) { + Options options{{"test", {{"dneut", 1.0}, {"diagnose", true}}}}; + + NeutralParallelDiffusion component("test", options, nullptr); + + // State without any collision frequencies + Options state{{"species", + {{"d", + {{"AA", 2.0}, + {"density", 1e19}, + {"temperature", 1.0}, + {"velocity", 0.0}, + {"momentum", 0.0}, + {"collision_frequencies", {{"d_d+_iz", 0.0}, {"d_d+_cx", 0.0}}}}}}}}; + + component.declareAllSpecies(SpeciesInformation({}, {"d"}, {}, {})); + component.transform(state); + + Options outputs = { + {"Tnorm", 1.0}, {"Nnorm", 1.0}, {"Omega_ci", 1.0}, {"Cs0", 1.0}, {"rho_s0", 1.0}}; + + component.outputVars(outputs); + + ASSERT_TRUE(outputs.isSet("Dd_Dpar")); + ASSERT_TRUE(outputs.isSet("Sd_Dpar")); + ASSERT_TRUE(outputs.isSet("Ed_Dpar")); +}