WIP Integration test

.github/workflows/ci.yml

@@ -53,6 +53,9 @@ jobs:
             build/tst/regression/outputs/cluster_hse/analytic_comparison.png
             build/tst/regression/outputs/cluster_tabular_cooling/convergence.png
             build/tst/regression/outputs/aniso_therm_cond_ring_conv/ring_convergence.png
+            build/tst/regression/outputs/aniso_therm_cond_gauss_conv/cond.png
+            build/tst/regression/outputs/diffusion/ohm.png
+            build/tst/regression/outputs/diffusion/visc.png
             build/tst/regression/outputs/field_loop/field_loop.png
             build/tst/regression/outputs/riemann_hydro/shock_tube.png
             build/tst/regression/outputs/turbulence/parthenon.hst

docs/input.md

@@ -69,15 +69,22 @@ conserved to primitive conversion if both are defined.
 
 #### Diffusive processes
 
-##### Anisotropic thermal conduction (required MHD)
+##### Isotropic (hydro and MHD )and anisotropic thermal conduction (only MHD)
 In the presence of magnetic fields thermal conduction is becoming anisotropic with the flux along
 the local magnetic field direction typically being much stronger than the flux perpendicular to the magnetic field.
 
 From a theoretical point of view, thermal conduction is included in the system of MHD equations by an additional
 term in the total energy equation:
 ```math
-\delta_t E + \nabla \cdot (... + \mathbf{F}) \quad \mathrm{with}\\
-\mathbf{F} = - \kappa \mathbf{\hat b} (\mathbf{\hat b \cdot \nabla T})
+\delta_t E + \nabla \cdot (... + \mathbf{F}_\mathrm{c})
+```
+where the full thermal conduction flux $`\mathbf{F}_\mathrm{c}`$ contains both the classic thermal conduction
+```math
+\mathbf{F}_\mathrm{classic} = - \kappa \mathbf{\hat b} (\mathbf{\hat b \cdot \nabla T})
+```
+as well as the saturated flux (as introduced by [^CM77])
+```math
+\mathbf{F}_\mathrm{sat} = - 5 \phi \rho^{-1/2} p^{3/2} \mathrm{sgn}(\mathbf{\hat b \cdot \nabla T}) \mathbf{\hat b}
 ```
 
 From an implementation point of view, two options implemented and can be configured within a `<diffusion>` block in the input file.
@@ -86,23 +93,52 @@ the integration step (before flux correction in case of AMR, and calculating the
 Moreover, they are implemented explicitly, i.e., they add a (potentially very restrictive) constraint to the timestep due to the scaling with $`\propto \Delta_x^2`$.
 Finally, we employ limiters for calculating the temperature gradients following Sharma & Hammett (2007)[^SH07].
 This prevents unphysical conduction against the gradient, which may be introduced because the off-axis gradients are not centered on the interfaces.
+Similarly, to account for the different nature of classic and saturated fluxes (parabolic and hyperbolic, respectively),
+we follow [^M+12] and use a smooth transition
+```math
+\mathbf{F}_\mathrm{c} = \frac{q}{q + F_\mathrm{classic}} \mathbf{F}_\mathrm{classic} \quad \mathrm{with} \quad q = 5 \phi \rho^{-1/2} p^{3/2}
+```
+and upwinding of the hyperbolic, saturated fluxes.
 
-To enable conduction, set
+To enable thermal conduction, set
 
 Parameter: `conduction` (string)
 - `none` : No thermal conduction
+- `isotropic` : Isotropic thermal conduction
+- `anisotropic` : Anisotropic thermal conduction
+
+In addition the coefficient (or diffusivity) needs to be set
+
+Parameter: `conduction_coeff` (string)
 - `spitzer` : Anisotropic thermal conduction with a temperature dependent classic Spitzer thermal conductivity
   $`\kappa (T) = c_\kappa T^{5/2} \mathrm{erg/s/K/cm}`$ and
-  $`c_\kappa`$ being constant prefactor (set via `diffusion/spitzer_cond_in_erg_by_s_K_cm` with a default value of $`4.6\times10^{-7}`$). Note, as indicated by the units in the input parameter name, this kind of thermal conductivity requires a full set of units
+  $`c_\kappa`$ being constant prefactor (set via the additional `diffusion/spitzer_cond_in_erg_by_s_K_cm` parameter with a default value of $`4.6\times10^{-7}`$). Note, as indicated by the units in the input parameter name, this kind of thermal conductivity requires a full set of units
   to be defined for the simulation.
-- `thermal_diff` : Contrary to a temperature dependent conductivity, a simple thermal diffusivity can be used instead for which
+- `fixed` : Contrary to a temperature dependent conductivity, a simple thermal diffusivity can be used instead for which
 the conduction flux is $`\mathbf{F} = - \chi \rho \mathbf{\hat b} (\mathbf{\hat b \cdot \nabla \frac{p_\mathrm{th}}{\rho}})`$
-Here, the strength, $`\chi`$, is controlled via the `thermal_diff_coeff_code` parameter in code units.
+Here, the strength, $`\chi`$, is controlled via the additional `thermal_diff_coeff_code` parameter in code units.
 Given the dimensions of $`L^2/T`$ it is referred to a thermal diffusivity rather than thermal conductivity.
 
+Parameter: `conduction_sat_phi` (float)
+- Default value 0.3\
+Factor to account for the uncertainty in the estimated of saturated fluxes, see [^CM77].
+Default value corresponds to the typical value used in literature and goes back to [^MMM80] and [^BM82].
+
+
 [^SH07]:
     P. Sharma and G. W. Hammett, "Preserving monotonicity in anisotropic diffusion," Journal of Computational Physics, vol. 227, no. 1, Art. no. 1, 2007, doi: https://doi.org/10.1016/j.jcp.2007.07.026.
 
+[^M+12]:
+    A. Mignone, C. Zanni, P. Tzeferacos, B. van Straalen, P. Colella, and G. Bodo, “THE PLUTO CODE FOR ADAPTIVE MESH COMPUTATIONS IN ASTROPHYSICAL FLUID DYNAMICS,” The Astrophysical Journal Supplement Series, vol. 198, Art. no. 1, Dec. 2011, doi: https://doi.org/10.1088/0067-0049/198/1/7
+
+[^CM77]:
+    L. Cowie and C. F. McKee, “The evaporation of spherical clouds in a hot gas. I. Classical and saturated mass loss rates.,” , vol. 211, pp. 135–146, Jan. 1977, doi: https://doi.org/10.1086/154911
+
+[^MMM80]:
+    C. E. Max, C. F. McKee, and W. C. Mead, “A model for laser driven ablative implosions,” The Physics of Fluids, vol. 23, Art. no. 8, 1980, doi: https://doi.org/10.1063/1.863183
+
+[^BM82]:
+    S. A. Balbus and C. F. McKee, “The evaporation of spherical clouds in a hot gas. III - Suprathermal evaporation,” , vol. 252, pp. 529–552, Jan. 1982, doi: https://doi.org/10.1086/159581
 
 ### Additional MHD options in `<hydro>` block
 

docs/turbulence.md

@@ -0,0 +1,68 @@
+# Turbulence problem generator
+
+The problem generator has been refactored from its implementation in K-Athena.
+Until the documentation is fully adapted, please consult https://gitlab.com/pgrete/kathena/-/wikis/turbulence for general information.
+
+A sample input file is provided in [turbulence.in](../inputs/turbulence.in)
+
+## Blob injection
+
+In order to study the evolution of (cold) clouds in a turbulent environment,
+blobs can be injected to the simulation.
+The injection is a one off mechanism controlled by the following parameters
+
+```
+<problem/turbulence>
+# only one of the following three conditions can be set at a given time
+inject_once_at_time = -1.0
+inject_once_at_cycle = -1
+inject_once_on_restart = false # should not be set in the input file, but only via the command line upon restart
+
+inject_n_blobs = -1 # number of blob to inject
+
+# then for the given number of blobs follow parameters need to be given (starting to count with 0)
+inject_blob_radius_0 = ... # float, in code length units, no default value
+inject_blob_loc_0 = ...,...,... # location vector of three comma-separated floats, in code length units, no default value
+inject_blob_chi_0 = ... # float, dimensionless, no default value, density ratio to existing value
+
+inject_blob_radius_1 = ...
+...
+```
+
+In practice, this will result in blobs being seeded at a given time, cycle, or upon restart
+by adjusting the density within the blob's radius by a factor of $\chi$ and
+at the same time adjusting the temperature by a factor of $1/\chi$ so that the
+blob remain in pressure equilibrium.
+
+While this is an action that is performed once, it can be repeated upon restart (or a later
+time) by resetting the variables.
+
+A current restriction is that the blobs cannot be seeded across a domain boundary (i.e.,
+the periodicity of the box is not taken into account).
+
+## Rescaling **not recommended*
+
+*The rescaling described in the following is generally not recommended, as it result in a
+state that is not naturally reached.
+Moreover, given the artificial nature of a hard reset, some time after the rescaling is
+required for the system to readjust.
+
+For non-isothermal simulations, the plasma will eventually heat up over time due to dissipation.
+One possibility to remove that extra heat (or add heat), is to rescale the temperature in the simulation.
+This can be done via the following parameters:
+
+```
+<problem/turbulence>
+# only one of the following three conditions can be set at a given time
+rescale_once_at_time = -1.0
+rescale_once_at_cycle = -1
+rescale_once_on_restart = false # should not be set in the input file, but only via the command line upon restart
+
+rescale_to_rms_Ms = -1.0
+```
+
+As the parameters suggest, rescaling is a one off action (though it can be repeated when
+the parameters are set again for a subsequent restart).
+The density and velocity field are not changed, only the specific internal energy is
+adjusted so that the volume-weighted root mean squared sonic Mach number matches
+the target value.

inputs/diffusion.in

@@ -1,9 +1,9 @@
 # AthenaPK - a performance portable block structured AMR MHD code
-# Copyright (c) 2021, Athena Parthenon Collaboration. All rights reserved.
+# Copyright (c) 2021-2023, Athena Parthenon Collaboration. All rights reserved.
 # Licensed under the BSD 3-Clause License (the "LICENSE");
 
 <comment>
-problem   = Thermal diffusion setup
+problem   = Diffusion setup (for thermal, momentum and Ohmic diffusion tests)
 
 <job>
 problem_id = diffusion
@@ -14,7 +14,8 @@ Bx = 1.0     # Bx for x1 step function (permutated for iprobs in other direction
 By = 0.0     # By for x1 step function (permutated for iprobs in other directions)
 
 #iprob = 10   # Diffusion of Gaussian profile in x1 direction
-sigma = 0.1  # standard deviation of Gaussian for iprob=10
+t0 = 0.5      # Temporal offset for initial Gaussian profile
+amp = 1e-6    # Amplitude of Gaussian profile
 
 iprob = 20   # ring diffusion in x1-x2 plane; 21 for x2-x3 plane; 22 for x3-x1 plane
 
@@ -59,8 +60,17 @@ reconstruction = dc
 gamma = 2.0
 
 <diffusion>
-conduction = thermal_diff
+integrator = rkl2
+conduction = anisotropic
+conduction_coeff = fixed
 thermal_diff_coeff_code = 0.01
+viscosity = none            # none (disabled), isotropic, or anisotropic
+viscosity_coeff = fixed
+mom_diff_coeff_code = 0.25
+resistivity = none          # none (disabled) or isotropic
+resistivity_coeff = fixed
+ohm_diff_coeff_code = 0.25
+rkl2_max_dt_ratio = 100.0
 
 <parthenon/output0>
 file_type = hdf5

src/CMakeLists.txt

@@ -7,18 +7,27 @@ add_executable(
         eos/adiabatic_glmmhd.cpp
         units.hpp
         eos/adiabatic_hydro.cpp
-        hydro/diffusion/diffusion.hpp
         hydro/diffusion/conduction.cpp
+        hydro/diffusion/diffusion.cpp
+        hydro/diffusion/diffusion.hpp
+        hydro/diffusion/resistivity.cpp
+        hydro/diffusion/viscosity.cpp
         hydro/hydro_driver.cpp
         hydro/hydro.cpp
         hydro/glmmhd/dedner_source.cpp
         hydro/prolongation/custom_ops.hpp
         hydro/srcterms/gravitational_field.hpp
         hydro/srcterms/tabular_cooling.hpp
         hydro/srcterms/tabular_cooling.cpp
+        outputs/per_block.cpp
         refinement/gradient.cpp
         refinement/other.cpp
+        tracers/tracers.cpp
+        tracers/tracers.hpp
         utils/few_modes_ft.cpp
+        utils/few_modes_ft.hpp
+        utils/interpolation.hpp
+        utils/robust.hpp
 )
 
 add_subdirectory(pgen)