Add isotropic viscosity and resistivity (#89)

* Add simple anisotropic step function test

* Separate FillDerived and EstimateTimestep in driver in prep for STS list

* Add diffflux parameter

* Add RKL2 STS task list

* Add calc of RKL2 stages

* Remove unncessary register for rkl2

* Adopt STS RKL2 variable naming

* Move calc of dt_diff into PreStep

* Make tlim an argument for diff step test

* Adjust RKL2 conv test to gaussian profile

* Add conv panel to conv plot

* auto-format

* rename diffusion integrator parameter

* Add isotropic thermal conduction

* Add isotropic cond to conv test

* Add RKL2 conv test

* Add new dt max ratio for rkl2 param

* Add prolongation and fluxcorrect to RKL2 task list

* Use base container as active STS container (workaround some AMR bug for using prolong/restric with non-base containers)

* Add isotropic Spitzer thermal conduction timestep

* Calc isotropic, non-const thermal diff

* Fix calc of saturated heat flux

* Add LimO3 limiter

* Add limo3 convergence

* Fix LimO3 recon

* Fix saturated conduction prefactor

* Remove calc of saturated conduction from cond coeff

* Add upwinded saturated conduction in x-dir

* Add saturated conduction prefactor

* Add x2 and x3 sat cond fluxes

* Increase default rkl2 ratio to 400 and allow flux correction for all integrators

* Remove parabolic timestep constraint for saturated conduction limit regime

* Add perturb to cloud pgen

* Add perturb to B (knowing this is not great...)

* Revert "Add perturb to B (knowing this is not great...)"

This reverts commit 1d0cb19.

* Revert "Add perturb to cloud pgen"

This reverts commit 6df018f.

* Limit cooling to upper bound of TFloor and cooling table cutoff

* Add oblique B field

* Update coords and driver

* Fix test cases and add success check

* Add isotropic shear viscosity

* Add viscosity test problem

* Add viscosity convergence test

* Add Ohmic resistivity

* Remove visc pgen and move to diffusion pgen

* Add resis. conv test to diffusion one

* Add linwave3d decay diffusion test and fix parabolic dt

* Fix test thresholds

* Add changelog and readme

* Cleanup leftover code

* Add doc

* Fix interface

* Fix more interface

* Fix even more interfaces

* And fixing a typo

* Increase ctest timeout

* Address BWO review comments

.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,6 +3,7 @@
 ## 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/...)

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

docs/input.md

@@ -69,6 +69,62 @@ conserved to primitive conversion if both are defined.
 
 #### 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.
+Note that if this limit is enforced the `dt=` shown on the terminal might not be the actual
+`dt` taken in the code as the limit is currently enforced only after the output
+has been printed.
+
+[^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 +196,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)

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-2024, 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
@@ -60,11 +60,17 @@ reconstruction = dc
 gamma = 2.0
 
 <diffusion>
-integrator = unsplit
+integrator = rkl2
 conduction = anisotropic
 conduction_coeff = fixed
 thermal_diff_coeff_code = 0.01
-rkl2_max_dt_ratio = 400.0
+viscosity = none            # none (disabled), isotropic, or anisotropic
+viscosity_coeff = fixed
+mom_diff_coeff_code = 0.25
+resistivity = none          # none (disabled) or ohmic
+resistivity_coeff = fixed
+ohm_diff_coeff_code = 0.25
+rkl2_max_dt_ratio = 100.0
 
 <parthenon/output0>
 file_type = hdf5

src/CMakeLists.txt

@@ -7,9 +7,11 @@ add_executable(
         eos/adiabatic_glmmhd.cpp
         units.hpp
         eos/adiabatic_hydro.cpp
+        hydro/diffusion/conduction.cpp
         hydro/diffusion/diffusion.cpp
         hydro/diffusion/diffusion.hpp
-        hydro/diffusion/conduction.cpp
+        hydro/diffusion/resistivity.cpp
+        hydro/diffusion/viscosity.cpp
         hydro/hydro_driver.cpp
         hydro/hydro.cpp
         hydro/glmmhd/dedner_source.cpp

src/hydro/diffusion/conduction.cpp

@@ -179,8 +179,8 @@ Real EstimateConductionTimestep(MeshData<Real> *md) {
         },
         Kokkos::Min<Real>(min_dt_cond));
   }
-
-  return fac * min_dt_cond;
+  const auto &cfl_diff = hydro_pkg->Param<Real>("cfl_diff");
+  return cfl_diff * fac * min_dt_cond;
 }
 
 //---------------------------------------------------------------------------------------

src/hydro/diffusion/diffusion.cpp

@@ -1,6 +1,6 @@
 //========================================================================================
 // AthenaPK - a performance portable block structured AMR astrophysical MHD code.
-// Copyright (c) 2021, Athena-Parthenon Collaboration. All rights reserved.
+// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved.
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
 //! \file diffusion.cpp
@@ -27,5 +27,27 @@ TaskStatus CalcDiffFluxes(StateDescriptor *hydro_pkg, MeshData<Real> *md) {
       ThermalFluxGeneral(md);
     }
   }
+  const auto &viscosity = hydro_pkg->Param<Viscosity>("viscosity");
+  if (viscosity != Viscosity::none) {
+    const auto &mom_diff = hydro_pkg->Param<MomentumDiffusivity>("mom_diff");
+
+    if (viscosity == Viscosity::isotropic &&
+        mom_diff.GetCoeffType() == ViscosityCoeff::fixed) {
+      MomentumDiffFluxIsoFixed(md);
+    } else {
+      MomentumDiffFluxGeneral(md);
+    }
+  }
+  const auto &resistivity = hydro_pkg->Param<Resistivity>("resistivity");
+  if (resistivity != Resistivity::none) {
+    const auto &ohm_diff = hydro_pkg->Param<OhmicDiffusivity>("ohm_diff");
+
+    if (resistivity == Resistivity::ohmic &&
+        ohm_diff.GetCoeffType() == ResistivityCoeff::fixed) {
+      OhmicDiffFluxIsoFixed(md);
+    } else {
+      OhmicDiffFluxGeneral(md);
+    }
+  }
   return TaskStatus::complete;
 }

src/hydro/diffusion/diffusion.hpp

@@ -99,6 +99,71 @@ void ThermalFluxIsoFixed(MeshData<Real> *md);
 //! Calculate thermal conduction (general case incl. anisotropic and saturated)
 void ThermalFluxGeneral(MeshData<Real> *md);
 
+struct MomentumDiffusivity {
+ private:
+  Real mbar_, me_, kb_;
+  Viscosity viscosity_;
+  ViscosityCoeff viscosity_coeff_type_;
+  // "free" coefficient/prefactor. Value depends on viscosity set in the constructor.
+  Real coeff_;
+
+ public:
+  KOKKOS_INLINE_FUNCTION
+  MomentumDiffusivity(Viscosity viscosity, ViscosityCoeff viscosity_coeff_type,
+                      Real coeff, Real mbar, Real me, Real kb)
+      : viscosity_(viscosity), viscosity_coeff_type_(viscosity_coeff_type), coeff_(coeff),
+        mbar_(mbar), me_(me), kb_(kb) {}
+
+  KOKKOS_INLINE_FUNCTION
+  Real Get(const Real pres, const Real rho) const;
+
+  KOKKOS_INLINE_FUNCTION
+  Viscosity GetType() const { return viscosity_; }
+
+  KOKKOS_INLINE_FUNCTION
+  ViscosityCoeff GetCoeffType() const { return viscosity_coeff_type_; }
+};
+
+Real EstimateViscosityTimestep(MeshData<Real> *md);
+
+//! Calculate isotropic viscosity with fixed coefficient
+void MomentumDiffFluxIsoFixed(MeshData<Real> *md);
+//! Calculate viscosity (general case incl. anisotropic)
+void MomentumDiffFluxGeneral(MeshData<Real> *md);
+
+struct OhmicDiffusivity {
+ private:
+  Real mbar_, me_, kb_;
+  Resistivity resistivity_;
+  ResistivityCoeff resistivity_coeff_type_;
+  // "free" coefficient/prefactor. Value depends on resistivity set in the constructor.
+  Real coeff_;
+
+ public:
+  KOKKOS_INLINE_FUNCTION
+  OhmicDiffusivity(Resistivity resistivity, ResistivityCoeff resistivity_coeff_type,
+                   Real coeff, Real mbar, Real me, Real kb)
+      : resistivity_(resistivity), resistivity_coeff_type_(resistivity_coeff_type),
+        coeff_(coeff), mbar_(mbar), me_(me), kb_(kb) {}
+
+  KOKKOS_INLINE_FUNCTION
+  Real Get(const Real pres, const Real rho) const;
+
+  KOKKOS_INLINE_FUNCTION
+  Resistivity GetType() const { return resistivity_; }
+
+  KOKKOS_INLINE_FUNCTION
+  ResistivityCoeff GetCoeffType() const { return resistivity_coeff_type_; }
+};
+
+Real EstimateResistivityTimestep(MeshData<Real> *md);
+
+//! Calculate isotropic resistivity with fixed coefficient
+void OhmicDiffFluxIsoFixed(MeshData<Real> *md);
+
+//! Calculate resistivity (general case incl. Spitzer)
+void OhmicDiffFluxGeneral(MeshData<Real> *md);
+
 // Calculate all diffusion fluxes, i.e., update the .flux views in md
 TaskStatus CalcDiffFluxes(StateDescriptor *hydro_pkg, MeshData<Real> *md);