Merge branch 'main' into exacb

.devcontainer/devcontainer.json

@@ -1,7 +1,7 @@
 // devcontainer.json
     {
       "name": "athenapk-dev",
-      "image": "ghcr.io/parthenon-hpc-lab/cuda11.6-mpi-hdf5-ascent",
+      "image": "ghcr.io/parthenon-hpc-lab/cuda11.6-noascent",
       // disable Dockerfile for now
       //"build": {
       //  // Path is relative to the devcontainer.json file.

.github/workflows/ci.yml

@@ -42,7 +42,7 @@ jobs:
           cd build
           # Pick GPU with most available memory
           export CUDA_VISIBLE_DEVICES=$(nvidia-smi --query-gpu=memory.free,index --format=csv,nounits,noheader | sort -nr | head -1 | awk '{ print $NF }')
-          ctest -L ${{ matrix.parallel }}
+          ctest -L ${{ matrix.parallel }} --timeout 3600
       - uses: actions/upload-artifact@v3
         if: ${{ always() }}
         with:
@@ -57,6 +57,8 @@ jobs:
             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.out1.hst

CHANGELOG.md

@@ -3,15 +3,22 @@
 ## Current develop (i.e., `main` branch)
 
 ### Added (new features/APIs/variables/...)
+- [[PR 89]](https://github.com/parthenon-hpc-lab/athenapk/pull/89) Add viscosity and resistivity
 - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 122]](https://github.com/parthenon-hpc-lab/athenapk/pull/122) Fixed sqrt(4pi) factor in CGS Gauss unit and add unit doc
+- [[PR 119]](https://github.com/parthenon-hpc-lab/athenapk/pull/119) Fixed Athena++ paper test case for KHI pgen. Added turbulence pgen doc.
 - [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97) Fixed Schure cooling curve. Removed SD one. Added description of cooling function conventions.
 - [[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/...)
+- [[PR 128]](https://github.com/parthenon-hpc-lab/athenapk/pull/128) Fixed `dt_diff` in RKL2
 
 ### Infrastructure
+- [[PR 128]](https://github.com/parthenon-hpc-lab/athenapk/pull/128) Bump Parthenon to support `dn` based outputs
+- [[PR 124]](https://github.com/parthenon-hpc-lab/athenapk/pull/124) Bump Kokkos 4.4.1 (and Parthenon to include view-of-view fix)
+- [[PR 117]](https://github.com/parthenon-hpc-lab/athenapk/pull/117) Update devcontainer.json to latest CI container
 - [[PR 114]](https://github.com/parthenon-hpc-lab/athenapk/pull/114) Bump Parthenon 24.08 and Kokkos to 4.4.00
 - [[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)
@@ -20,6 +27,8 @@
 ### Removed (removing behavior/API/varaibles/...)
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 124]](https://github.com/parthenon-hpc-lab/athenapk/pull/124) Enrolling custom boundary conditions changed
+  - Boundary conditions can now be enrolled using a string that can be subsequently be used in the input file (see, e.g., cloud problem generator)
 - [[PR 114]](https://github.com/parthenon-hpc-lab/athenapk/pull/114) Bump Parthenon 24.08 and Kokkos to 4.4.00
   - Changed signature of `UserWorkBeforeOutput` to include `SimTime` as last paramter
   - Fixes bitwise idential restarts for AMR simulations (the derefinement counter is now included)

README.md

@@ -14,8 +14,13 @@ 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
-  - isotropic and anisotropic thermal conduction
-  - operator-split, second-order RKL2 supertimestepping for diffusive terms
+  - diffusion processes
+    - isotropic and anisotropic thermal conduction
+    - viscosity
+    - resistivity
+  - diffusion integrator
+    - unsplit
+    - operator-split, second-order RKL2 supertimestepping
   - 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
@@ -126,21 +131,8 @@ the `file_type = hdf5` format, see
 [VisIt](https://wci.llnl.gov/simulation/computer-codes/visit/).
 In ParaView, select the "XDMF Reader" when prompted.
 
-2. With [yt](https://yt-project.org/) -- though currently through a custom frontend
-that is not yet part of the main yt branch and, thus, has to be installed manually, e.g.,
-as follows:
-```bash
-cd ~/src # or any other folder of choice
-git clone https://github.com/forrestglines/yt.git
-cd yt
-git checkout parthenon-frontend
-
-# If you're using conda or virtualenv
-pip install -e .
-# OR alternatively, if you using the plain Python environment
-pip install --user -e .
-```
-Afterwards, `*.phdf` files can be read as usual with `yt.load()`.
+2. With [yt](https://yt-project.org/)
+As of versions >=4.4 `*.phdf` files can be read as usual with `yt.load()`.
 
 3. Using [Ascent](https://github.com/Alpine-DAV/ascent) (for in situ visualization and analysis).
 This requires Ascent to be installed/available at compile time of AthenaPK.

docs/README.md

@@ -0,0 +1,37 @@
+# AthenaPK documentation
+
+Note that we're aware the rendering of equations in markdown on GitHub for the documentation
+is currently not great.
+Eventually, the docs will transition to Read the Docs (or similar).
+
+## Overview
+
+The documentation currently includes
+
+- [Configuring solvers in the input file](input.md)
+  - [An extended overview of various conventions for optically thin cooling implementations](cooling_notes.md)
+  - [Notebooks to calculate cooling tables from literature](cooling)
+- [Brief notes on developing code for AthenaPK](development.md)
+- [How to add a custom/user problem generator](pgen.md)
+- [Units](units.md)
+- Detailed descriptions of more complex problem generators
+  - [Galaxy Cluster and Cluster-like Problem Setup](cluster.md)
+  - [Driven turbulence](turbulence.md)
+
+## Tutorial
+
+An AthenaPK tutorial was given as part of the **Towards exascale-ready astrophysics**
+workshop https://indico3-jsc.fz-juelich.de/event/169/ taking place 25-27 Sep 2024 online.
+
+The material is currently located at https://github.com/pgrete/athenapk_tutorial
+
+While the instructions for building the code are specific to the workshop environment
+the tutorial itself should translate directly to other environments/systems.
+
+## Parthenon documenation
+
+Many paramters/options are directly controlled through the Parthenon framework
+(both with regard to building and in the input file).
+
+While the [Parthenon documenation](https://parthenon-hpc-lab.github.io/parthenon) is
+more geared towards developers it also contains useful information for users.

docs/input.md

@@ -67,8 +67,65 @@ To control the floors, following parameters can be set in the `<hydro>` block:
 *Note* the pressure floor will take precedence over the temperature floor in the
 conserved to primitive conversion if both are defined.
 
+#### Units
+
+See(here)[units.md].
+
 #### Diffusive processes
 
+Diffusive processes in AthenaPK can be configured in the `<diffusion>` block of the input file.
+```
+<diffusion>
+integrator = unsplit       # alternatively: rkl2 (for rkl2 integrator (operator split integrator)
+#rkl2_max_dt_ratio = 100.0 # limits the ratio between the parabolic and hyperbolic timesteps (only used for RKL2 operator split integrator)
+#cfl = 1.0                 # Additional safety factor applied in the caluclation of the diffusive timestep (used in both unsplit and RKL2 integration schemes). Defaults to hyperbolic cfl.
+
+conduction = anisotropic               # none (disabled), or isotropic, or anisotropic
+conduction_coeff = fixed               # alternative: spitzer
+thermal_diff_coeff_code = 0.01         # fixed coefficent in code units (code_length^2/code_time)
+#spitzer_cond_in_erg_by_s_K_cm = 4.6e7 # spitzer coefficient in cgs units (requires definition of a unit system)
+#conduction_sat_phi = 0.3              # fudge factor to account for uncertainties in saturated fluxes
+
+
+viscosity = none            # none (disabled) or isotropic
+viscosity_coeff = fixed
+mom_diff_coeff_code = 0.25  # fixed coefficent of the kinetmatic viscosity in code units (code_length^2/code_time)
+
+resistivity = none          # none (disabled) or ohmic
+resistivity_coeff = fixed
+ohm_diff_coeff_code = 0.25  # fixed coefficent of the magnetic (ohmic) diffusivity code units (code_length^2/code_time)
+```
+(An)isotropic thermal conduction (with fixed or Spitzer coefficient), and isotropic viscosity and
+resistivity with fixed coefficient are currently implemented.
+They can be integrated in an unsplit manner or operator split using a second-order accurate RKL2
+supertimestepping algorithm.
+More details are described in the following.
+
+#### Integrators
+
+Diffusive processes can be integrated in either an unsplit
+fashion (`diffusion/integrator=unsplit`) or operator split using a second-order accurate
+RKL2 super timestepping algorithm (`diffusion/integrator=rkl2`) following [^M+14].
+
+In the unsplit case, the diffusive processes are included at the end of every stage in
+the main integration loop and the global timestep is limited accordingly.
+A separate CFL can be set for the diffusive processes via `diffusion/cfl=...`, which
+defaults to the hyperbolic value if not set.
+
+In the RKL2 case, the global timestep is not limited by the diffusive processes by default.
+However, as reported by [^V+17] a large number of stages
+($`s \approx \sqrt(\Delta t_{hyp}/\Delta t_{par}) \geq 20`$) in the supertimestepping
+(in combination with anisotropic, limited diffusion) may lead to a loss in accuracy, which
+is why the difference between hyperbolic and parabolic timesteps can be limited by
+`diffusion/rkl2_max_dt_ratio=...` and a warning is shown if the ratio is above 400.
+
+[^M+14]:
+    C. D. Meyer, D. S. Balsara, and T. D. Aslam, “A stabilized Runge–Kutta–Legendre method for explicit super-time-stepping of parabolic and mixed equations,” Journal of Computational Physics, vol. 257, pp. 594–626, 2014, doi: https://doi.org/10.1016/j.jcp.2013.08.021.
+
+[^V+17]:
+    B. Vaidya, D. Prasad, A. Mignone, P. Sharma, and L. Rickler, “Scalable explicit implementation of anisotropic diffusion with Runge–Kutta–Legendre super-time stepping,” Monthly Notices of the Royal Astronomical Society, vol. 472, no. 3, pp. 3147–3160, 2017, doi: 10.1093/mnras/stx2176.
+
+
 ##### 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.
@@ -140,6 +197,29 @@ Default value corresponds to the typical value used in literature and goes back
 [^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
 
+#### Viscosity/Momentum diffusion
+
+Only isotropic viscosity with a (spatially and temporally) fixed coefficient in code units
+(`code_length`^2/`code_time`) is currently implemented.
+To enable set (in the `<diffusion>` block)
+```
+viscosity = isotropic
+viscosity_coeff = fixed
+mom_diff_coeff_code = 0.25  # fixed coefficent of the kinetmatic viscosity in code units (code_length^2/code_time)
+```
+
+#### Resistivity/Ohmic diffusion
+
+Only resistivity with a (spatially and temporally) fixed coefficient in code units
+(`code_length`^2/`code_time`)is currently implemented.
+To enable set (in the `<diffusion>` block)
+```
+resistivity = ohmic
+resistivity_coeff = fixed
+ohm_diff_coeff_code = 0.25  # fixed coefficent of the magnetic (ohmic) diffusivity code units (code_length^2/code_time)
+```
+
+
 ### Additional MHD options in `<hydro>` block
 
 Parameter: `glmmhd_source` (string)

docs/pgen.md

@@ -1,4 +1,4 @@
-# Problem generators
+# Custom problem generators
 
 Different simulation setups in AthenaPK are controlled via the so-called problem generators.
 New problem generators can easily be added and we are happy to accept and merge contibuted problem

docs/turbulence.md

@@ -0,0 +1,107 @@
+# Driven turbulence simulations
+
+The turbulence problem generator uses explicit inverse Fourier transformations (iFTs)
+on each meshblock in order to reduce communication during the iFT.
+Thus, it is only efficient if comparatively few modes are used (say < 100).
+
+Quite generally, driven turbulence simulations start from uniform initial conditions
+(uniform density and pressure, some initial magnetic field configuration in case of an
+MHD setup, and the fluid at rest) and reach a state of stationary, isotropic (or anisotropic
+depending on the strength of the background magnetic field) turbulence after one to few
+large eddy turnover times (again depending on the background magnetic field strength).
+The large eddy turnover time is usually defined as `T = L/U` with `L` being the scale
+of the largest eddies and `U` the root mean square Mach number in the stationary regime.
+
+The current implementation uses the following forcing spectrum
+`(k/k_peak)^2 * (2 - (k/k_peak)^2)`.
+Here, `k_peak` is the peak wavenumber of the forcing spectrum. It is related the scales of the largest eddies as
+`L = 1/k_f` given that a box size of 1 is currently assumed/hardcoded.
+
+## Problem setup
+
+An example parameter file can be found in `inputs/turbulence.in`.
+
+A typical setup contains the following blocks in the input file:
+
+```
+<job>
+problem_id = turbulence
+
+<problem/turbulence>
+rho0         = 1.0      # initial mean density
+p0           = 1.0      # initial mean pressure
+b0           = 0.01     # initial magnetic field strength
+b_config     = 0        # 0 - net flux; 1 - no net flux uniform B; 2 - non net flux sin B; 4 - field loop
+kpeak        = 2.0      # characteristic wavenumber
+corr_time    = 1.0      # autocorrelation time of the OU forcing process
+rseed        = 20190729 # random seed of the OU forcing process
+sol_weight   = 1.0      # solenoidal weight of the acceleration field
+accel_rms    = 0.5      # root mean square value of the acceleration field
+num_modes    = 30       # number of wavemodes
+
+<modes>
+k_1_0	= +2
+k_1_1	= -1
+k_1_2	= +0
+k_2_0	= +1
+...
+```
+
+The following parameters can be changed to control both the initial state:
+
+- `rho0` initial mean density
+- `p0` initial mean thermal pressure
+- `b0` initial mean magnetic field strength
+- `b_config`
+  - `0`: net flux case (uniform B_x)
+  - `1`: no net flux case (uniform B_x with changing sign in each half of the box)
+  - `2`: no net flux with initial sinosoidal B_x field
+  - `3`: deprecated
+  - `4`: closed field loop/cylinder in the box (in x-y plane) located at
+    - `x0=0.5` (default)
+    - `y0=0.5` (default)
+    - `z0=0.5` (default)
+    - and radius `loop_rad=0.25`
+
+as well as the driving field:
+
+- `kpeak` peak wavenumber of the forcing spectrum. Make sure to update the wavemodes to match `kpeak`, see below.
+- `corr_time` autocorrelation time of the acceleration field (in code units).
+Using delta-in-time forcing, i.e., a very low value, is discouraged, see [Grete et al. 2018 ApJL](https://iopscience.iop.org/article/10.3847/2041-8213/aac0f5).
+- `rseed` random seed for the OU process. Only change for new simulation, but keep unchanged for restarting simulations.
+- `sol_weight` solenoidal weight of the acceleration field. `1.0` is purely solenoidal/rotational and `0.0` is purely dilatational/compressive. Any value between `0.0` and `1.0` is possible. The parameter is related to the resulting rotational power in the 3D acceleration field as
+`1. - ((1-sol_weight)^2/(1-2*sol_weight+3*sol_weight^2))`, see eq (9) in [Federrath et al. 2010 A&A](
+https://doi.org/10.1051/0004-6361/200912437).
+- `accel_rms` root mean square value of the acceleration (controls the "strength" of the forcing field)
+- `num_modes` number of wavemodes that are specified in the `<modes>` section of the parameter file.
+The modes are specified manually as an explicit inverse FT is performed and only modes set are included (all others are assumed to be 0).
+This is done to make the global inverse FT possible without any
+expensive communication between blocks but this becomes excessively
+expensiv for large number of modes.
+Typically using a few tens of modes is a good choice in practice.
+In order to generate a set of modes run the `inputs/generate_fmturb_modes.py` script and replace
+the corresponding parts of the parameter file with the output of the script.
+Within the script, the top three variables (`k_peak`, `k_high`, and `k_low`) need to be adjusted in
+order to generate a complete set (i.e., all) of wavemodes.
+Important, the `k_peak` in the script should match the `k_peak` set
+in the input file.
+Alternatively, wavemodes can be chosen/defined manually, e.g., if not all wavemodes are desired or
+only individual modes should be forced.
+
+## Typical results
+
+The results shown here are obtained from running simulations with the parameters given in the next section.
+
+### High level temporal evolution
+![image](img/turb_evol.png)
+
+### Power spectra
+![image](img/turb_spec.png)
+
+### Consistency of acceleration field
+As each meshblock does a full iFT of all modes the following slices from a run with 8 meshblocks 
+illustrate that there is no discontinuities at the meshblock boundary.
+
+Plot shows x-, y-, and z-acceleration (in rows top to bottom) slices in the x-, y-, and z-direction (in columns from left to right).
+
+![image](img/turb_acc.png)