Merge branch 'main' into bwibking/reflecting-bc

.devcontainer/devcontainer.json

@@ -0,0 +1,40 @@
+// devcontainer.json
+    {
+      "name": "athenapk-dev",
+      "image": "ghcr.io/parthenon-hpc-lab/cuda11.6-mpi-hdf5-ascent",
+      // disable Dockerfile for now
+      //"build": {
+      //  // Path is relative to the devcontainer.json file.
+      //  "dockerfile": "Dockerfile"
+      //},
+      "hostRequirements": {
+        "cpus": 4
+      },
+      "customizations": {
+        "vscode": {
+          "settings": {},
+          "extensions": [
+            "-ms-vscode.cpptools",
+            "llvm-vs-code-extensions.vscode-clangd",
+            "github.vscode-pull-request-github",
+            "ms-python.python",
+            "ms-toolsai.jupyter",
+            "ms-vscode.live-server",
+            "ms-azuretools.vscode-docker",
+            "swyddfa.esbonio",
+            "tomoki1207.pdf",
+            "ms-vscode.cmake-tools",
+            "ms-vsliveshare.vsliveshare"
+        ]
+        }
+      },
+      "remoteEnv": {
+        "PATH": "${containerEnv:PATH}:/usr/local/hdf5/parallel/bin",
+        "OMPI_MCA_opal_warn_on_missing_libcuda": "0"
+      },
+      //"remoteUser": "ubuntu",
+      // we need to manually checkout the submodules,
+      // but VSCode may try to configure CMake before they are fully checked-out.
+      // workaround TBD
+      "postCreateCommand": "git submodule update --init"
+    }

.github/workflows/ci.yml

@@ -53,6 +53,7 @@ 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/field_loop/field_loop.png
             build/tst/regression/outputs/riemann_hydro/shock_tube.png
             build/tst/regression/outputs/turbulence/parthenon.hst

CHANGELOG.md

@@ -0,0 +1,26 @@
+# Changelog
+
+## Current develop (i.e., `main` branch)
+
+### Added (new features/APIs/variables/...)
+- [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping
+
+### Changed (changing behavior/API/variables/...)
+- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
+
+### Fixed (not changing behavior/API/variables/...)
+
+### Infrastructure
+- [[PR 112]](https://github.com/parthenon-hpc-lab/athenapk/pull/112) Add dev container configuration
+- [[PR 105]](https://github.com/parthenon-hpc-lab/athenapk/pull/105) Bump Parthenon to latest develop (2024-03-13)
+- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Added `CHANGELOG.md`
+
+### Removed (removing behavior/API/varaibles/...)
+
+### Incompatibilities (i.e. breaking changes)
+- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
+  - Updated access to block dimension: `pmb->block_size.nx1` -> `pmb->block_size.nx(X1DIR)` (and similarly x2 and x3)
+  - Update access to mesh size: `pmesh->mesh_size.x1max` -> `pmesh->mesh_size.xmax(X1DIR)` (and similarly x2, x3, and min)
+  - Updated Parthenon `GradMinMod` signature for custom prolongation ops
+  - `GetBlockPointer` returns a raw pointer not a shared one (and updated interfaces to use raw pointers rather than shared ones)
+

README.md

@@ -14,7 +14,8 @@ Current features include
   - HLLE (hydro and MHD), HLLC (hydro), and HLLD (MHD) Riemann solvers
   - adiabatic equation of state
   - MHD based on hyperbolic divergence cleaning following Dedner+ 2002
-  - anisotropic thermal conduction
+  - isotropic and anisotropic thermal conduction
+  - operator-split, second-order RKL2 supertimestepping for diffusive terms
   - optically thin cooling based on tabulated cooling tables with either Townsend 2009 exact integration or operator-split subcycling
 - static and adaptive mesh refinement
 - problem generators for

docs/cluster.md

@@ -76,14 +76,14 @@ $$
     \frac{3 H_0^2}{8 \pi G}.
 $$
 
-Parameters for the HERNQUIST BCG are controlled via:
+Parameters for the `HERNQUIST` BCG are controlled via:
 ```
 <problem/cluster/gravity>
 
 m_bcg_s = 0.001 # in code_mass
 r_bcg_s = 0.004 # in code_length
 ```
-where a HERNQUIST profile adds a gravitational acceleration defined by
+where a `HERNQUIST` profile adds a gravitational acceleration defined by
 
 $$
  g_{BCG}(r) = G \frac{ M_{BCG} }{R_{BCG}^2} \frac{1}{\left( 1 + \frac{r}{R_{BCG}}\right)^2}
@@ -216,7 +216,7 @@ triggering_mode = COLD_GAS # or NONE, BOOSTED_BONDI, BONDI_SCHAYE
 ```
 where `triggering_mode=NONE` will disable AGN triggering. 
 
-With BOOSTED_BONDI accretion, the mass rate of accretion follows
+With `BOOSTED_BONDI` accretion, the mass rate of accretion follows
 
 $$
 \dot{M} = \alpha \frac { 2 \pi G^2 M^2_{SMBH} \hat {\rho} } {
@@ -235,7 +235,7 @@ m_smbh = 1.0e-06 # in code_mass
 accretion_radius = 0.001 # in code_length
 bondi_alpha= 100.0 # unitless
 ```
-With BONDI_SCHAYE accretion, the `$\alpha$` used for BOOSTED_BONDI accretion is modified to depend on the number density following:
+With `BONDI_SCHAYE` accretion, the $\alpha$ used for `BOOSTED_BONDI` accretion is modified to depend on the number density following:
 
 $$
 \alpha =
@@ -252,7 +252,7 @@ bondi_n0= 2.9379989445851786e+72 # in 1/code_length**3
 bondi_beta= 2.0 # unitless
 ```
 
-With both BOOSTED_BONDI and BONDI_SCHAYE accretion, mass is removed from each
+With both `BOOSTED_BONDI` and `BONDI_SCHAYE` accretion, mass is removed from each
 cell within the accretion zone at a mass weighted rate. E.g. the mass in each
 cell within the accretion region changes by
 ```
@@ -264,7 +264,7 @@ unchanged. Thus velocities and temperatures will increase where mass is
 removed.
 
 
-With COLD_GAS accretion, the accretion rate becomes the total mass within the accretion zone equal to or
+With `COLD_GAS` accretion, the accretion rate becomes the total mass within the accretion zone equal to or
 below a defined cold temperature threshold divided by a defined accretion
 timescale. The temperature threshold and accretion timescale are defined by
 ```
@@ -360,13 +360,17 @@ velocity $v_{jet}$ can be set via
 kinetic_jet_temperature = 1e7 # K
 ```
 However, $T_{jet}$ and $v_{jet}$ must be non-negative and fulfill
+
 $$
-v_{jet} = \sqrt{ 2 \left ( \epsilon c^2 - (1 - \epsilon) \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right} \right ) }
+v_{jet} = \sqrt{ 2 \left ( \epsilon c^2 - (1 - \epsilon) \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right) } \right ) }
 $$
+
 to ensure that the sum of rest mass energy, thermal energy, and kinetic energy of the new gas sums to $\dot{M} c^2$. Note that these equations places limits  on $T_{jet}$ and $v_{jet}$, specifically
+
 $$
-v_{jet} \leq c \sqrt{ 2 \epsilon } \qquad \text{and} \qquad \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right} \leq c^2 \frac{ \epsilon}{1 - \epsilon}
+v_{jet} \leq c \sqrt{ 2 \epsilon } \qquad \text{and} \qquad \frac{k_B T_{jet}}{ \mu m_h \left( \gamma - 1 \right) } \leq c^2 \frac{ \epsilon}{1 - \epsilon}
 $$
+
 If the above equations are not satified then an exception will be thrown at
 initialization. If neither $T_{jet}$ nor $v_{jet}$ are specified, then
 $v_{jet}$ will be computed assuming $T_{jet}=0$ and a warning will be given
@@ -429,19 +433,19 @@ where the injected magnetic field follows
 
 $$
 \begin{align}
-\mathcal{B}_r      &=\mathcal{B}_0 2 \frac{h r}{\ell^2} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\\\
-\mathcal{B}_\theta &=\mathcal{B}_0 \alpha \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right ) } \\\\
-\mathcal{B}_h      &=\mathcal{B}_0 2 \left( 1 - \frac{r^2}{\ell^2} \right ) \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\\\
+\mathcal{B}\_r      &= \mathcal{B}\_0 2 \frac{h r}{\ell^2} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\
+\mathcal{B}\_\theta &= \mathcal{B}\_0 \alpha \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right ) } \\
+\mathcal{B}\_h      &= \mathcal{B}\_0 2 \left( 1 - \frac{r^2}{\ell^2} \right ) \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )}
 \end{align}
 $$
 
 which has  the corresponding vector potential field
 
 $$
 \begin{align}
-\mathcal{A}_r &= 0 \\\\
-\mathcal{A}_{\theta} &= \mathcal{B}_0 \ell \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\\\
-\mathcal{A}_h &= \mathcal{B}_0 \ell \frac{\alpha}{2}\exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )}
+\mathcal{A}\_r &= 0 \\
+\mathcal{A}\_{\theta} &= \mathcal{B}\_0 \ell \frac{r}{\ell} \exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )} \\
+\mathcal{A}\_h &= \mathcal{B}\_0 \ell \frac{\alpha}{2}\exp{ \left ( \frac{-r^2 - h^2}{\ell^2} \right )}
 \end{align}
 $$
 
@@ -463,7 +467,7 @@ $$
 where $\dot{\rho}_B$ is set to
 
 $$
-\dot{\rho}_B = \frac{3 \pi}{2} \frac{\dot{M} \left ( 1 - \epsilon \right ) f_{magnetic}}{\ell^3}
+\dot{\rho}_B = \frac{3 \pi}{2} \frac{\dot{M} \left ( 1 - \epsilon \right ) f\_\mathrm{magnetic}}{\ell^3}
 $$
 
 so that the total mass injected matches the accreted mass propotioned to magnetic feedback.
@@ -485,13 +489,17 @@ according to their ratio.
 #### Simple field loop (donut) feedback
 
 Magnetic energy is injected according to the following simple potential
+
 $$
-A_h(r, \theta, h) = B_0 L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+A_h(r, \theta, h) = B_0 L \exp^\left ( -r^2/L^2 \right) \qquad \mathrm{for} \qquad h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
 $$
+
 resultig in a magnetic field configuration of
+
 $$
-B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right)$ for $h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
+B_\theta(r, \theta, h) = 2 B_0 r /L \exp^\left ( -r^2/L^2 \right) \qquad \mathrm{for} \qquad h_\mathrm{offset} \leq |h| \leq h_\mathrm{offset} + h_\mathrm{thickness}
 $$
+
 with all other components being zero.
 
 
@@ -526,7 +534,7 @@ fixed_mass_rate = 1.0
 
 ## SNIA Feedback
 
-Following [Prasad 2020](doi.org/10.1093/mnras/112.2.195), AthenaPK can inject
+Following [Prasad 2020](https://doi.org/10.3847/1538-4357/abc33c), AthenaPK can inject
 mass and energy from type Ia supernovae following the mass profile of the BCG.
 This SNIA feedback can be configured with
 ```
@@ -562,4 +570,4 @@ temperature_threshold = 2e4 # in K
 ```
 
 Note that all parameters need to be specified explicitly for the feedback to work
-(i.e., no hidden default values).
+(i.e., no hidden default values).

docs/development.md

@@ -7,3 +7,9 @@ All pull requests are checked automatically if the code is properly formatted.
 
 The build target `format-athenapk` calls both formatters and automatically format all changes, i.e., before committing changes simply call `make format-athenapk` (or similar).
 Alternatively, leave a comment with `@par-hermes format` in the open pull request to format the code directly on GitHub.
+
+## Dev container
+
+You can open a [GitHub Codespace](https://docs.github.com/en/codespaces) or use VSCode to automatically open local Docker container using the CUDA CI container image.
+
+If you have the [VSCode Dev Container extension](https://marketplace.visualstudio.com/items?itemName=ms-vscode-remote.remote-containers) installed, on opening this repository in VSCode, it will automatically prompt to ask if you want to re-open this repository in a container.

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}`$ which assumes a fully ionized hydrogen plasma [^CM77] with $`\ln \lambda = 40`$ approximating ICM conditions). 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
 

inputs/diffusion.in

@@ -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,11 @@ reconstruction = dc
 gamma = 2.0
 
 <diffusion>
-conduction = thermal_diff
+integrator = unsplit
+conduction = anisotropic
+conduction_coeff = fixed
 thermal_diff_coeff_code = 0.01
+rkl2_max_dt_ratio = 400.0
 
 <parthenon/output0>
 file_type = hdf5

src/CMakeLists.txt

@@ -8,6 +8,7 @@ add_executable(
         eos/adiabatic_glmmhd.cpp
         units.hpp
         eos/adiabatic_hydro.cpp
+        hydro/diffusion/diffusion.cpp
         hydro/diffusion/diffusion.hpp
         hydro/diffusion/conduction.cpp
         hydro/hydro_driver.cpp