diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 1d8fbe03..8fcfcba2 100644
--- a/.github/workflows/ci.yml
+++ b/.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
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 369e8579..65d55b73 100644
--- a/CHANGELOG.md
+++ b/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/...)
diff --git a/README.md b/README.md
index 13d7280a..fd66e824 100644
--- a/README.md
+++ b/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
diff --git a/docs/input.md b/docs/input.md
index c51aa8e7..d5e48d4f 100644
--- a/docs/input.md
+++ b/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)
diff --git a/inputs/diffusion.in b/inputs/diffusion.in
index 44d02bcb..915dd088 100644
--- a/inputs/diffusion.in
+++ b/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
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ddcf09b7..ae2a521a 100644
--- a/src/CMakeLists.txt
+++ b/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
diff --git a/src/hydro/diffusion/conduction.cpp b/src/hydro/diffusion/conduction.cpp
index d6e57d30..40256024 100644
--- a/src/hydro/diffusion/conduction.cpp
+++ b/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;
 }
 
 //---------------------------------------------------------------------------------------
diff --git a/src/hydro/diffusion/diffusion.cpp b/src/hydro/diffusion/diffusion.cpp
index 20617ca2..d8b1c59d 100644
--- a/src/hydro/diffusion/diffusion.cpp
+++ b/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;
 }
diff --git a/src/hydro/diffusion/diffusion.hpp b/src/hydro/diffusion/diffusion.hpp
index 7d522902..9e03f25a 100644
--- a/src/hydro/diffusion/diffusion.hpp
+++ b/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);
 
diff --git a/src/hydro/diffusion/resistivity.cpp b/src/hydro/diffusion/resistivity.cpp
new file mode 100644
index 00000000..da8bf9ad
--- /dev/null
+++ b/src/hydro/diffusion/resistivity.cpp
@@ -0,0 +1,240 @@
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR astrophysical MHD code.
+// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved.
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+//! \file resistivity.cpp
+//! \brief
+
+// Parthenon headers
+#include <cmath>
+#include <parthenon/package.hpp>
+
+// AthenaPK headers
+#include "../../main.hpp"
+#include "config.hpp"
+#include "diffusion.hpp"
+#include "kokkos_abstraction.hpp"
+#include "utils/error_checking.hpp"
+
+using namespace parthenon::package::prelude;
+
+KOKKOS_INLINE_FUNCTION
+Real OhmicDiffusivity::Get(const Real pres, const Real rho) const {
+  if (resistivity_coeff_type_ == ResistivityCoeff::fixed) {
+    return coeff_;
+  } else if (resistivity_coeff_type_ == ResistivityCoeff::spitzer) {
+    PARTHENON_FAIL("needs impl");
+  } else {
+    PARTHENON_FAIL("Unknown Resistivity coeff");
+  }
+}
+
+Real EstimateResistivityTimestep(MeshData<Real> *md) {
+  // get to package via first block in Meshdata (which exists by construction)
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+  const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+
+  IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
+  IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
+
+  Real min_dt_resist = std::numeric_limits<Real>::max();
+  const auto ndim = prim_pack.GetNdim();
+
+  Real fac = 0.5;
+  if (ndim == 2) {
+    fac = 0.25;
+  } else if (ndim == 3) {
+    fac = 1.0 / 6.0;
+  }
+
+  const auto &ohm_diff = hydro_pkg->Param<OhmicDiffusivity>("ohm_diff");
+
+  if (ohm_diff.GetType() == Resistivity::ohmic &&
+      ohm_diff.GetCoeffType() == ResistivityCoeff::fixed) {
+    // TODO(pgrete): once mindx is properly calculated before this loop, we can get rid of
+    // it entirely.
+    // Using 0.0 as parameters rho and p as they're not used anyway for a fixed coeff.
+    const auto ohm_diff_coeff = ohm_diff.Get(0.0, 0.0);
+    Kokkos::parallel_reduce(
+        "EstimateResistivityTimestep (ohmic fixed)",
+        Kokkos::MDRangePolicy<Kokkos::Rank<4>>(
+            DevExecSpace(), {0, kb.s, jb.s, ib.s},
+            {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1},
+            {1, 1, 1, ib.e + 1 - ib.s}),
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &min_dt) {
+          const auto &coords = prim_pack.GetCoords(b);
+          min_dt =
+              fmin(min_dt, SQR(coords.Dxc<1>(k, j, i)) / (ohm_diff_coeff + TINY_NUMBER));
+          if (ndim >= 2) {
+            min_dt = fmin(min_dt,
+                          SQR(coords.Dxc<2>(k, j, i)) / (ohm_diff_coeff + TINY_NUMBER));
+          }
+          if (ndim >= 3) {
+            min_dt = fmin(min_dt,
+                          SQR(coords.Dxc<3>(k, j, i)) / (ohm_diff_coeff + TINY_NUMBER));
+          }
+        },
+        Kokkos::Min<Real>(min_dt_resist));
+  } else {
+    PARTHENON_THROW("Needs impl.");
+  }
+
+  const auto &cfl_diff = hydro_pkg->Param<Real>("cfl_diff");
+  return cfl_diff * fac * min_dt_resist;
+}
+
+//---------------------------------------------------------------------------------------
+//! Calculate isotropic resistivity with fixed coefficient
+
+void OhmicDiffFluxIsoFixed(MeshData<Real> *md) {
+  auto pmb = md->GetBlockData(0)->GetBlockPointer();
+  IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  std::vector<parthenon::MetadataFlag> flags_ind({Metadata::Independent});
+  auto cons_pack = md->PackVariablesAndFluxes(flags_ind);
+  auto hydro_pkg = pmb->packages.Get("Hydro");
+
+  auto const &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+
+  const int ndim = pmb->pmy_mesh->ndim;
+
+  const auto &ohm_diff = hydro_pkg->Param<OhmicDiffusivity>("ohm_diff");
+  // Using fixed and uniform coefficient so it's safe to get it outside the kernel.
+  // Using 0.0 as parameters rho and p as they're not used anyway for a fixed coeff.
+  const auto eta = ohm_diff.Get(0.0, 0.0);
+
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "Resist. X1 fluxes (ohmic)", DevExecSpace(), 0,
+      cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e + 1,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        const auto &coords = prim_pack.GetCoords(b);
+        auto &cons = cons_pack(b);
+        const auto &prim = prim_pack(b);
+
+        // Face centered current densities
+        // j2 = d3B1 - d1B3
+        const auto d3B1 =
+            ndim > 2 ? (0.5 * (prim(IB1, k + 1, j, i - 1) + prim(IB1, k + 1, j, i)) -
+                        0.5 * (prim(IB1, k - 1, j, i - 1) + prim(IB1, k - 1, j, i))) /
+                           (2.0 * coords.Dxf<3, 1>(k, j, i))
+                     : 0.0;
+
+        const auto d1B3 =
+            (prim(IB3, k, j, i) - prim(IB3, k, j, i - 1)) / coords.Dxc<1>(k, j, i);
+
+        const auto j2 = d3B1 - d1B3;
+
+        // j3 = d1B2 - d2B1
+        const auto d1B2 =
+            (prim(IB2, k, j, i) - prim(IB2, k, j, i - 1)) / coords.Dxc<1>(k, j, i);
+
+        const auto d2B1 =
+            ndim > 1 ? (0.5 * (prim(IB1, k, j + 1, i - 1) + prim(IB1, k, j + 1, i)) -
+                        0.5 * (prim(IB1, k, j - 1, i - 1) + prim(IB1, k, j - 1, i))) /
+                           (2.0 * coords.Dxf<2, 1>(k, j, i))
+                     : 0.0;
+
+        const auto j3 = d1B2 - d2B1;
+
+        cons.flux(X1DIR, IB2, k, j, i) += -eta * j3;
+        cons.flux(X1DIR, IB3, k, j, i) += eta * j2;
+        cons.flux(X1DIR, IEN, k, j, i) +=
+            0.5 * eta *
+            ((prim(IB3, k, j, i - 1) + prim(IB3, k, j, i)) * j2 -
+             (prim(IB2, k, j, i - 1) + prim(IB2, k, j, i)) * j3);
+      });
+
+  if (ndim < 2) {
+    return;
+  }
+
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "Resist. X2 fluxes (ohmic)", parthenon::DevExecSpace(), 0,
+      cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e + 1, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        const auto &coords = prim_pack.GetCoords(b);
+        auto &cons = cons_pack(b);
+        const auto &prim = prim_pack(b);
+
+        // Face centered current densities
+        // j3 = d1B2 - d2B1
+        const auto d1B2 = (0.5 * (prim(IB2, k, j - 1, i + 1) + prim(IB2, k, j, i + 1)) -
+                           0.5 * (prim(IB2, k, j - 1, i - 1) + prim(IB2, k, j, i - 1))) /
+                          (2.0 * coords.Dxf<1, 2>(k, j, i));
+
+        const auto d2B1 =
+            (prim(IB1, k, j, i) - prim(IB1, k, j - 1, i)) / coords.Dxc<2>(k, j, i);
+
+        const auto j3 = d1B2 - d2B1;
+
+        // j1 = d2B3 - d3B2
+        const auto d2B3 =
+            (prim(IB3, k, j, i) - prim(IB3, k, j - 1, i)) / coords.Dxc<2>(k, j, i);
+
+        const auto d3B2 =
+            ndim > 2 ? (0.5 * (prim(IB2, k + 1, j - 1, i) + prim(IB2, k + 1, j, i)) -
+                        0.5 * (prim(IB2, k - 1, j - 1, i) + prim(IB2, k - 1, j, i))) /
+                           (2.0 * coords.Dxf<3, 2>(k, j, i))
+                     : 0.0;
+
+        const auto j1 = d2B3 - d3B2;
+
+        cons.flux(X2DIR, IB1, k, j, i) += eta * j3;
+        cons.flux(X2DIR, IB3, k, j, i) += -eta * j1;
+        cons.flux(X2DIR, IEN, k, j, i) +=
+            0.5 * eta *
+            ((prim(IB1, k, j - 1, i) + prim(IB1, k, j, i)) * j3 -
+             (prim(IB3, k, j - 1, i) + prim(IB3, k, j, i)) * j1);
+      });
+
+  if (ndim < 3) {
+    return;
+  }
+
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "Resist. X3 fluxes (ohmic)", parthenon::DevExecSpace(), 0,
+      cons_pack.GetDim(5) - 1, kb.s, kb.e + 1, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        const auto &coords = prim_pack.GetCoords(b);
+        auto &cons = cons_pack(b);
+        const auto &prim = prim_pack(b);
+
+        // Face centered current densities
+        // j1 = d2B3 - d3B2
+        const auto d2B3 = (0.5 * (prim(IB3, k - 1, j + 1, i) + prim(IB3, k, j + 1, i)) -
+                           0.5 * (prim(IB3, k - 1, j - 1, i) + prim(IB3, k, j - 1, i))) /
+                          (2.0 * coords.Dxf<2, 3>(k, j, i));
+
+        const auto d3B2 =
+            (prim(IB2, k, j, i) - prim(IB2, k - 1, j, i)) / coords.Dxc<3>(k, j, i);
+
+        const auto j1 = d2B3 - d3B2;
+
+        // j2 = d3B1 - d1B3
+        const auto d3B1 =
+            (prim(IB1, k, j, i) - prim(IB1, k - 1, j, i)) / coords.Dxc<3>(k, j, i);
+
+        const auto d1B3 = (0.5 * (prim(IB3, k - 1, j, i + 1) + prim(IB3, k, j, i + 1)) -
+                           0.5 * (prim(IB3, k - 1, j, i - 1) + prim(IB3, k, j, i - 1))) /
+                          (2.0 * coords.Dxf<1, 3>(k, j, i));
+
+        const auto j2 = d3B1 - d1B3;
+
+        cons.flux(X3DIR, IB1, k, j, i) += -eta * j2;
+        cons.flux(X3DIR, IB2, k, j, i) += eta * j1;
+        cons.flux(X3DIR, IEN, k, j, i) +=
+            0.5 * eta *
+            ((prim(IB2, k - 1, j, i) + prim(IB2, k, j, i)) * j1 -
+             (prim(IB1, k - 1, j, i) + prim(IB1, k, j, i)) * j2);
+      });
+}
+
+//---------------------------------------------------------------------------------------
+//! TODO(pgrete) Calculate Ohmic diffusion, general case, e.g., with varying (Spitzer)
+//! coefficient
+
+void OhmicDiffFluxGeneral(MeshData<Real> *md) { PARTHENON_THROW("Needs impl."); }
\ No newline at end of file
diff --git a/src/hydro/diffusion/viscosity.cpp b/src/hydro/diffusion/viscosity.cpp
new file mode 100644
index 00000000..f222c869
--- /dev/null
+++ b/src/hydro/diffusion/viscosity.cpp
@@ -0,0 +1,295 @@
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR astrophysical MHD code.
+// Copyright (c) 2021-2024, Athena-Parthenon Collaboration. All rights reserved.
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+// AthenaXXX astrophysical plasma code
+// Copyright(C) 2020 James M. Stone <jmstone@ias.edu> and the Athena code team
+// Licensed under the 3-clause BSD License (the "LICENSE")
+//========================================================================================
+//! \file viscosity.cpp
+//! \brief
+
+// Parthenon headers
+#include <cmath>
+#include <parthenon/package.hpp>
+
+// AthenaPK headers
+#include "../../main.hpp"
+#include "config.hpp"
+#include "diffusion.hpp"
+#include "kokkos_abstraction.hpp"
+#include "utils/error_checking.hpp"
+
+using namespace parthenon::package::prelude;
+
+KOKKOS_INLINE_FUNCTION
+Real MomentumDiffusivity::Get(const Real pres, const Real rho) const {
+  if (viscosity_coeff_type_ == ViscosityCoeff::fixed) {
+    return coeff_;
+  } else {
+    PARTHENON_FAIL("Unknown viscosity coeff");
+  }
+}
+
+Real EstimateViscosityTimestep(MeshData<Real> *md) {
+  // get to package via first block in Meshdata (which exists by construction)
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+  const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+
+  IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
+  IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
+
+  Real min_dt_visc = std::numeric_limits<Real>::max();
+  const auto ndim = prim_pack.GetNdim();
+
+  Real fac = 0.5;
+  if (ndim == 2) {
+    fac = 0.25;
+  } else if (ndim == 3) {
+    fac = 1.0 / 6.0;
+  }
+
+  const auto gm1 = hydro_pkg->Param<Real>("AdiabaticIndex");
+  const auto &mom_diff = hydro_pkg->Param<MomentumDiffusivity>("mom_diff");
+
+  if (mom_diff.GetType() == Viscosity::isotropic &&
+      mom_diff.GetCoeffType() == ViscosityCoeff::fixed) {
+    // TODO(pgrete): once mindx is properly calculated before this loop, we can get rid of
+    // it entirely.
+    // Using 0.0 as parameters rho and p as they're not used anyway for a fixed coeff.
+    const auto mom_diff_coeff = mom_diff.Get(0.0, 0.0);
+    Kokkos::parallel_reduce(
+        "EstimateViscosityTimestep (iso fixed)",
+        Kokkos::MDRangePolicy<Kokkos::Rank<4>>(
+            DevExecSpace(), {0, kb.s, jb.s, ib.s},
+            {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1},
+            {1, 1, 1, ib.e + 1 - ib.s}),
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &min_dt) {
+          const auto &coords = prim_pack.GetCoords(b);
+          min_dt =
+              fmin(min_dt, SQR(coords.Dxc<1>(k, j, i)) / (mom_diff_coeff + TINY_NUMBER));
+          if (ndim >= 2) {
+            min_dt = fmin(min_dt,
+                          SQR(coords.Dxc<2>(k, j, i)) / (mom_diff_coeff + TINY_NUMBER));
+          }
+          if (ndim >= 3) {
+            min_dt = fmin(min_dt,
+                          SQR(coords.Dxc<3>(k, j, i)) / (mom_diff_coeff + TINY_NUMBER));
+          }
+        },
+        Kokkos::Min<Real>(min_dt_visc));
+  } else {
+    PARTHENON_THROW("Needs impl.");
+  }
+
+  const auto &cfl_diff = hydro_pkg->Param<Real>("cfl_diff");
+  return cfl_diff * fac * min_dt_visc;
+}
+
+//---------------------------------------------------------------------------------------
+//! Calculate isotropic viscosity with fixed coefficient
+
+void MomentumDiffFluxIsoFixed(MeshData<Real> *md) {
+  auto pmb = md->GetBlockData(0)->GetBlockPointer();
+  IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  std::vector<parthenon::MetadataFlag> flags_ind({Metadata::Independent});
+  auto cons_pack = md->PackVariablesAndFluxes(flags_ind);
+  auto hydro_pkg = pmb->packages.Get("Hydro");
+
+  auto const &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+
+  const int ndim = pmb->pmy_mesh->ndim;
+
+  const auto &mom_diff = hydro_pkg->Param<MomentumDiffusivity>("mom_diff");
+  // Using fixed and uniform coefficient so it's safe to get it outside the kernel.
+  // Using 0.0 as parameters rho and p as they're not used anyway for a fixed coeff.
+  const auto nu = mom_diff.Get(0.0, 0.0);
+
+  const int scratch_level =
+      hydro_pkg->Param<int>("scratch_level"); // 0 is actual scratch (tiny); 1 is HBM
+  const int nx1 = pmb->cellbounds.ncellsi(IndexDomain::entire);
+
+  size_t scratch_size_in_bytes = parthenon::ScratchPad1D<Real>::shmem_size(nx1) * 3;
+
+  parthenon::par_for_outer(
+      DEFAULT_OUTER_LOOP_PATTERN, "Visc. X1 fluxes (iso)", DevExecSpace(),
+      scratch_size_in_bytes, scratch_level, 0, cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s,
+      jb.e,
+      KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k, const int j) {
+        const auto &coords = prim_pack.GetCoords(b);
+        auto &cons = cons_pack(b);
+        const auto &prim = prim_pack(b);
+        parthenon::ScratchPad1D<Real> fvx(member.team_scratch(scratch_level), nx1);
+        parthenon::ScratchPad1D<Real> fvy(member.team_scratch(scratch_level), nx1);
+        parthenon::ScratchPad1D<Real> fvz(member.team_scratch(scratch_level), nx1);
+
+        // Add [2(dVx/dx)-(2/3)dVx/dx, dVy/dx, dVz/dx]
+        par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
+          fvx(i) = 4.0 * (prim(IV1, k, j, i) - prim(IV1, k, j, i - 1)) /
+                   (3.0 * coords.Dxc<1>(i));
+          fvy(i) = (prim(IV2, k, j, i) - prim(IV2, k, j, i - 1)) / coords.Dxc<1>(i);
+          fvz(i) = (prim(IV3, k, j, i) - prim(IV3, k, j, i - 1)) / coords.Dxc<1>(i);
+        });
+        member.team_barrier();
+
+        // In 2D/3D Add [(-2/3)dVy/dy, dVx/dy, 0]
+        if (ndim > 1) {
+          par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
+            fvx(i) -= ((prim(IV2, k, j + 1, i) + prim(IV2, k, j + 1, i - 1)) -
+                       (prim(IV2, k, j - 1, i) + prim(IV2, k, j - 1, i - 1))) /
+                      (6.0 * coords.Dxc<2>(j));
+            fvy(i) += ((prim(IV1, k, j + 1, i) + prim(IV1, k, j + 1, i - 1)) -
+                       (prim(IV1, k, j - 1, i) + prim(IV1, k, j - 1, i - 1))) /
+                      (4.0 * coords.Dxc<2>(j));
+          });
+          member.team_barrier();
+        }
+
+        // In 3D Add [(-2/3)dVz/dz, 0,  dVx/dz]
+        if (ndim > 2) {
+          par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
+            fvx(i) -= ((prim(IV3, k + 1, j, i) + prim(IV3, k + 1, j, i - 1)) -
+                       (prim(IV3, k - 1, j, i) + prim(IV3, k - 1, j, i - 1))) /
+                      (6.0 * coords.Dxc<3>(k));
+            fvz(i) += ((prim(IV1, k + 1, j, i) + prim(IV1, k + 1, j, i - 1)) -
+                       (prim(IV1, k - 1, j, i) + prim(IV1, k - 1, j, i - 1))) /
+                      (4.0 * coords.Dxc<3>(k));
+          });
+          member.team_barrier();
+        }
+
+        // Sum viscous fluxes into fluxes of conserved variables; including energy fluxes
+        par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
+          Real nud = 0.5 * nu * (prim(IDN, k, j, i) + prim(IDN, k, j, i - 1));
+          cons.flux(X1DIR, IV1, k, j, i) -= nud * fvx(i);
+          cons.flux(X1DIR, IV2, k, j, i) -= nud * fvy(i);
+          cons.flux(X1DIR, IV3, k, j, i) -= nud * fvz(i);
+          cons.flux(X1DIR, IEN, k, j, i) -=
+              0.5 * nud *
+              ((prim(IV1, k, j, i - 1) + prim(IV1, k, j, i)) * fvx(i) +
+               (prim(IV2, k, j, i - 1) + prim(IV2, k, j, i)) * fvy(i) +
+               (prim(IV3, k, j, i - 1) + prim(IV3, k, j, i)) * fvz(i));
+        });
+      });
+
+  if (ndim < 2) {
+    return;
+  }
+  /* Compute viscous fluxes in 2-direction  --------------------------------------*/
+  parthenon::par_for_outer(
+      DEFAULT_OUTER_LOOP_PATTERN, "Visc. X2 fluxes (iso)", parthenon::DevExecSpace(),
+      scratch_size_in_bytes, scratch_level, 0, cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s,
+      jb.e + 1,
+      KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k, const int j) {
+        const auto &coords = prim_pack.GetCoords(b);
+        auto &cons = cons_pack(b);
+        const auto &prim = prim_pack(b);
+        parthenon::ScratchPad1D<Real> fvx(member.team_scratch(scratch_level), nx1);
+        parthenon::ScratchPad1D<Real> fvy(member.team_scratch(scratch_level), nx1);
+        parthenon::ScratchPad1D<Real> fvz(member.team_scratch(scratch_level), nx1);
+
+        // Add [(dVx/dy+dVy/dx), 2(dVy/dy)-(2/3)(dVx/dx+dVy/dy), dVz/dy]
+        par_for_inner(member, ib.s, ib.e, [&](const int i) {
+          fvx(i) = (prim(IV1, k, j, i) - prim(IV1, k, j - 1, i)) / coords.Dxc<2>(j) +
+                   ((prim(IV2, k, j, i + 1) + prim(IV2, k, j - 1, i + 1)) -
+                    (prim(IV2, k, j, i - 1) + prim(IV2, k, j - 1, i - 1))) /
+                       (4.0 * coords.Dxc<1>(i));
+          fvy(i) = (prim(IV2, k, j, i) - prim(IV2, k, j - 1, i)) * 4.0 /
+                       (3.0 * coords.Dxc<2>(j)) -
+                   ((prim(IV1, k, j, i + 1) + prim(IV1, k, j - 1, i + 1)) -
+                    (prim(IV1, k, j, i - 1) + prim(IV1, k, j - 1, i - 1))) /
+                       (6.0 * coords.Dxc<1>(i));
+          fvz(i) = (prim(IV3, k, j, i) - prim(IV3, k, j - 1, i)) / coords.Dxc<2>(j);
+        });
+        member.team_barrier();
+
+        // In 3D Add [0, (-2/3)dVz/dz, dVy/dz]
+        if (ndim > 2) {
+          par_for_inner(member, ib.s, ib.e, [&](const int i) {
+            fvy(i) -= ((prim(IV3, k + 1, j, i) + prim(IV3, k + 1, j - 1, i)) -
+                       (prim(IV3, k - 1, j, i) + prim(IV3, k - 1, j - 1, i))) /
+                      (6.0 * coords.Dxc<3>(k));
+            fvz(i) += ((prim(IV2, k + 1, j, i) + prim(IV2, k + 1, j - 1, i)) -
+                       (prim(IV2, k - 1, j, i) + prim(IV2, k - 1, j - 1, i))) /
+                      (4.0 * coords.Dxc<3>(k));
+          });
+          member.team_barrier();
+        }
+
+        // Sum viscous fluxes into fluxes of conserved variables; including energy fluxes
+        par_for_inner(member, ib.s, ib.e, [&](const int i) {
+          Real nud = 0.5 * nu * (prim(IDN, k, j, i) + prim(IDN, k, j - 1, i));
+          cons.flux(X2DIR, IV1, k, j, i) -= nud * fvx(i);
+          cons.flux(X2DIR, IV2, k, j, i) -= nud * fvy(i);
+          cons.flux(X2DIR, IV3, k, j, i) -= nud * fvz(i);
+          cons.flux(X2DIR, IEN, k, j, i) -=
+              0.5 * nud *
+              ((prim(IV1, k, j - 1, i) + prim(IV1, k, j, i)) * fvx(i) +
+               (prim(IV2, k, j - 1, i) + prim(IV2, k, j, i)) * fvy(i) +
+               (prim(IV3, k, j - 1, i) + prim(IV3, k, j, i)) * fvz(i));
+        });
+      });
+  /* Compute heat fluxes in 3-direction, 3D problem ONLY  ---------------------*/
+  if (ndim < 3) {
+    return;
+  }
+
+  parthenon::par_for_outer(
+      DEFAULT_OUTER_LOOP_PATTERN, "Visc. X3 fluxes (iso)", parthenon::DevExecSpace(),
+      scratch_size_in_bytes, scratch_level, 0, cons_pack.GetDim(5) - 1, kb.s, kb.e + 1,
+      jb.s, jb.e,
+      KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k, const int j) {
+        const auto &coords = prim_pack.GetCoords(b);
+        auto &cons = cons_pack(b);
+        const auto &prim = prim_pack(b);
+
+        parthenon::ScratchPad1D<Real> fvx(member.team_scratch(scratch_level), nx1);
+        parthenon::ScratchPad1D<Real> fvy(member.team_scratch(scratch_level), nx1);
+        parthenon::ScratchPad1D<Real> fvz(member.team_scratch(scratch_level), nx1);
+
+        // Add [(dVx/dz+dVz/dx), (dVy/dz+dVz/dy), 2(dVz/dz)-(2/3)(dVx/dx+dVy/dy+dVz/dz)]
+        par_for_inner(member, ib.s, ib.e, [&](const int i) {
+          fvx(i) = (prim(IV1, k, j, i) - prim(IV1, k - 1, j, i)) / coords.Dxc<3>(k) +
+                   ((prim(IV3, k, j, i + 1) + prim(IV3, k - 1, j, i + 1)) -
+                    (prim(IV3, k, j, i - 1) + prim(IV3, k - 1, j, i - 1))) /
+                       (4.0 * coords.Dxc<1>(i));
+          fvy(i) = (prim(IV2, k, j, i) - prim(IV2, k - 1, j, i)) / coords.Dxc<3>(k) +
+                   ((prim(IV3, k, j + 1, i) + prim(IV3, k - 1, j + 1, i)) -
+                    (prim(IV3, k, j - 1, i) + prim(IV3, k - 1, j - 1, i))) /
+                       (4.0 * coords.Dxc<2>(j));
+          fvz(i) = (prim(IV3, k, j, i) - prim(IV3, k - 1, j, i)) * 4.0 /
+                       (3.0 * coords.Dxc<3>(k)) -
+                   ((prim(IV1, k, j, i + 1) + prim(IV1, k - 1, j, i + 1)) -
+                    (prim(IV1, k, j, i - 1) + prim(IV1, k - 1, j, i - 1))) /
+                       (6.0 * coords.Dxc<1>(i)) -
+                   ((prim(IV2, k, j + 1, i) + prim(IV2, k - 1, j + 1, i)) -
+                    (prim(IV2, k, j - 1, i) + prim(IV2, k - 1, j - 1, i))) /
+                       (6.0 * coords.Dxc<2>(j));
+        });
+        member.team_barrier();
+
+        // Sum viscous fluxes into fluxes of conserved variables; including energy fluxes
+        par_for_inner(member, ib.s, ib.e, [&](const int i) {
+          Real nud = 0.5 * nu * (prim(IDN, k, j, i) + prim(IDN, k - 1, j, i));
+          cons.flux(X3DIR, IV1, k, j, i) -= nud * fvx(i);
+          cons.flux(X3DIR, IV2, k, j, i) -= nud * fvy(i);
+          cons.flux(X3DIR, IV3, k, j, i) -= nud * fvz(i);
+          cons.flux(X3DIR, IEN, k, j, i) -=
+              0.5 * nud *
+              ((prim(IV1, k - 1, j, i) + prim(IV1, k, j, i)) * fvx(i) +
+               (prim(IV2, k - 1, j, i) + prim(IV2, k, j, i)) * fvy(i) +
+               (prim(IV3, k - 1, j, i) + prim(IV3, k, j, i)) * fvz(i));
+        });
+      });
+}
+
+//---------------------------------------------------------------------------------------
+//! TODO(pgrete) Calculate momentum diffusion, general case, i.e., anisotropic and/or with
+//! varying coefficient
+
+void MomentumDiffFluxGeneral(MeshData<Real> *md) { PARTHENON_THROW("Needs impl."); }
diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp
index 54a97b7f..7bc95be7 100644
--- a/src/hydro/hydro.cpp
+++ b/src/hydro/hydro.cpp
@@ -64,13 +64,28 @@ void PreStepMeshUserWorkInLoop(Mesh *pmesh, ParameterInput *pin, SimTime &tm) {
   auto hydro_pkg = pmesh->block_list[0]->packages.Get("Hydro");
   const auto num_partitions = pmesh->DefaultNumPartitions();
 
-  if ((hydro_pkg->Param<DiffInt>("diffint") == DiffInt::rkl2) &&
-      (hydro_pkg->Param<Conduction>("conduction") != Conduction::none)) {
+  if (hydro_pkg->Param<DiffInt>("diffint") == DiffInt::rkl2) {
     auto dt_diff = std::numeric_limits<Real>::max();
-    for (auto i = 0; i < num_partitions; i++) {
-      auto &md = pmesh->mesh_data.GetOrAdd("base", i);
+    if (hydro_pkg->Param<Conduction>("conduction") != Conduction::none) {
+      for (auto i = 0; i < num_partitions; i++) {
+        auto &md = pmesh->mesh_data.GetOrAdd("base", i);
 
-      dt_diff = std::min(dt_diff, EstimateConductionTimestep(md.get()));
+        dt_diff = std::min(dt_diff, EstimateConductionTimestep(md.get()));
+      }
+    }
+    if (hydro_pkg->Param<Viscosity>("viscosity") != Viscosity::none) {
+      for (auto i = 0; i < num_partitions; i++) {
+        auto &md = pmesh->mesh_data.GetOrAdd("base", i);
+
+        dt_diff = std::min(dt_diff, EstimateViscosityTimestep(md.get()));
+      }
+    }
+    if (hydro_pkg->Param<Resistivity>("resistivity") != Resistivity::none) {
+      for (auto i = 0; i < num_partitions; i++) {
+        auto &md = pmesh->mesh_data.GetOrAdd("base", i);
+
+        dt_diff = std::min(dt_diff, EstimateResistivityTimestep(md.get()));
+      }
     }
 #ifdef MPI_PARALLEL
     PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &dt_diff, 1, MPI_PARTHENON_REAL,
@@ -544,6 +559,67 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     }
     pkg->AddParam<>("conduction", conduction);
 
+    auto viscosity = Viscosity::none;
+    auto viscosity_str = pin->GetOrAddString("diffusion", "viscosity", "none");
+    if (viscosity_str == "isotropic") {
+      viscosity = Viscosity::isotropic;
+    } else if (viscosity_str != "none") {
+      PARTHENON_FAIL("Unknown viscosity method. Options are: none, isotropic");
+    }
+    // If viscosity is enabled, process supported coefficients
+    if (viscosity != Viscosity::none) {
+      auto viscosity_coeff_str =
+          pin->GetOrAddString("diffusion", "viscosity_coeff", "none");
+      auto viscosity_coeff = ViscosityCoeff::none;
+
+      if (viscosity_coeff_str == "fixed") {
+        viscosity_coeff = ViscosityCoeff::fixed;
+        Real mom_diff_coeff_code = pin->GetReal("diffusion", "mom_diff_coeff_code");
+        auto mom_diff = MomentumDiffusivity(viscosity, viscosity_coeff,
+                                            mom_diff_coeff_code, 0.0, 0.0, 0.0);
+        pkg->AddParam<>("mom_diff", mom_diff);
+
+      } else {
+        PARTHENON_FAIL("Viscosity is enabled but no coefficient is set. Please "
+                       "set diffusion/viscosity_coeff to 'fixed' and "
+                       "diffusion/mom_diff_coeff_code to the desired value.");
+      }
+    }
+    pkg->AddParam<>("viscosity", viscosity);
+
+    auto resistivity = Resistivity::none;
+    auto resistivity_str = pin->GetOrAddString("diffusion", "resistivity", "none");
+    if (resistivity_str == "ohmic") {
+      resistivity = Resistivity::ohmic;
+    } else if (resistivity_str != "none") {
+      PARTHENON_FAIL("Unknown resistivity method. Options are: none, ohmic");
+    }
+    // If resistivity is enabled, process supported coefficients
+    if (resistivity != Resistivity::none) {
+      auto resistivity_coeff_str =
+          pin->GetOrAddString("diffusion", "resistivity_coeff", "none");
+      auto resistivity_coeff = ResistivityCoeff::none;
+
+      if (resistivity_coeff_str == "spitzer") {
+        // If this is implemented, check how the Spitzer coeff for thermal conduction is
+        // handled.
+        PARTHENON_FAIL("needs impl");
+
+      } else if (resistivity_coeff_str == "fixed") {
+        resistivity_coeff = ResistivityCoeff::fixed;
+        Real ohm_diff_coeff_code = pin->GetReal("diffusion", "ohm_diff_coeff_code");
+        auto ohm_diff = OhmicDiffusivity(resistivity, resistivity_coeff,
+                                         ohm_diff_coeff_code, 0.0, 0.0, 0.0);
+        pkg->AddParam<>("ohm_diff", ohm_diff);
+
+      } else {
+        PARTHENON_FAIL("Resistivity is enabled but no coefficient is set. Please "
+                       "set diffusion/resistivity_coeff to 'fixed' and "
+                       "diffusion/ohm_diff_coeff_code to the desired value.");
+      }
+    }
+    pkg->AddParam<>("resistivity", resistivity);
+
     auto diffint_str = pin->GetOrAddString("diffusion", "integrator", "none");
     auto diffint = DiffInt::none;
     if (diffint_str == "unsplit") {
@@ -558,6 +634,10 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     }
     if (diffint != DiffInt::none) {
       pkg->AddParam<Real>("dt_diff", 0.0, true); // diffusive timestep constraint
+      // As in Athena++ a cfl safety factor is also applied to the theoretical limit.
+      // By default it is equal to the hyperbolic cfl.
+      auto cfl_diff = pin->GetOrAddReal("diffusion", "cfl", pkg->Param<Real>("cfl"));
+      pkg->AddParam<>("cfl_diff", cfl_diff);
     }
     pkg->AddParam<>("diffint", diffint);
 
@@ -790,9 +870,16 @@ Real EstimateTimestep(MeshData<Real> *md) {
   }
 
   // For RKL2 STS, the diffusive timestep is calculated separately in the driver
-  if ((hydro_pkg->Param<DiffInt>("diffint") == DiffInt::unsplit) &&
-      (hydro_pkg->Param<Conduction>("conduction") != Conduction::none)) {
-    min_dt = std::min(min_dt, EstimateConductionTimestep(md));
+  if (hydro_pkg->Param<DiffInt>("diffint") == DiffInt::unsplit) {
+    if (hydro_pkg->Param<Conduction>("conduction") != Conduction::none) {
+      min_dt = std::min(min_dt, EstimateConductionTimestep(md));
+    }
+    if (hydro_pkg->Param<Viscosity>("viscosity") != Viscosity::none) {
+      min_dt = std::min(min_dt, EstimateViscosityTimestep(md));
+    }
+    if (hydro_pkg->Param<Resistivity>("resistivity") != Resistivity::none) {
+      min_dt = std::min(min_dt, EstimateResistivityTimestep(md));
+    }
   }
 
   if (ProblemEstimateTimestep != nullptr) {
diff --git a/src/main.cpp b/src/main.cpp
index 93f984bb..fa0c4dec 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -55,6 +55,7 @@ int main(int argc, char *argv[]) {
     pman.app_input->InitUserMeshData = linear_wave_mhd::InitUserMeshData;
     pman.app_input->ProblemGenerator = linear_wave_mhd::ProblemGenerator;
     pman.app_input->UserWorkAfterLoop = linear_wave_mhd::UserWorkAfterLoop;
+    Hydro::ProblemInitPackageData = linear_wave_mhd::ProblemInitPackageData;
   } else if (problem == "cpaw") {
     pman.app_input->InitUserMeshData = cpaw::InitUserMeshData;
     pman.app_input->ProblemGenerator = cpaw::ProblemGenerator;
diff --git a/src/main.hpp b/src/main.hpp
index 033118ad..4c81c307 100644
--- a/src/main.hpp
+++ b/src/main.hpp
@@ -1,5 +1,5 @@
 // AthenaPK - a performance portable block structured AMR MHD code
-// Copyright (c) 2020-2021, Athena Parthenon Collaboration. All rights reserved.
+// Copyright (c) 2020-2024, Athena Parthenon Collaboration. All rights reserved.
 // Licensed under the 3-Clause License (the "LICENSE")
 
 #ifndef MAIN_HPP_
@@ -37,6 +37,10 @@ enum class Fluid { undefined, euler, glmmhd };
 enum class Cooling { none, tabular };
 enum class Conduction { none, isotropic, anisotropic };
 enum class ConductionCoeff { none, fixed, spitzer };
+enum class Viscosity { none, isotropic };
+enum class ViscosityCoeff { none, fixed };
+enum class Resistivity { none, ohmic };
+enum class ResistivityCoeff { none, fixed, spitzer };
 enum class DiffInt { none, unsplit, rkl2 };
 
 enum class Hst { idx, ekin, emag, divb };
diff --git a/src/pgen/diffusion.cpp b/src/pgen/diffusion.cpp
index 24ae60f2..269926f8 100644
--- a/src/pgen/diffusion.cpp
+++ b/src/pgen/diffusion.cpp
@@ -1,6 +1,6 @@
 
 // AthenaPK - a performance portable block structured AMR MHD code
-// Copyright (c) 2021-2023, Athena Parthenon Collaboration. All rights reserved.
+// Copyright (c) 2021-2024, Athena Parthenon Collaboration. All rights reserved.
 // Licensed under the 3-Clause License (the "LICENSE")
 
 // Parthenon headers
@@ -11,11 +11,13 @@
 
 // AthenaPK headers
 #include "../main.hpp"
+#include "utils/error_checking.hpp"
 
 namespace diffusion {
 using namespace parthenon::driver::prelude;
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
+  auto hydro_pkg = pmb->packages.Get("Hydro");
   IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
   IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
   IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
@@ -23,24 +25,39 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto &mbd = pmb->meshblock_data.Get();
   auto &u = mbd->Get("cons").data;
 
+  const auto gamma = pin->GetReal("hydro", "gamma");
+  const bool mhd_enabled = hydro_pkg->Param<Fluid>("fluid") == Fluid::glmmhd;
+
   const auto Bx = pin->GetOrAddReal("problem/diffusion", "Bx", 0.0);
   const auto By = pin->GetOrAddReal("problem/diffusion", "By", 0.0);
 
   const auto iprob = pin->GetInteger("problem/diffusion", "iprob");
+  PARTHENON_REQUIRE_THROWS(mhd_enabled || !(iprob == 0 || iprob == 1 || iprob == 2 ||
+                                            iprob == 10 || iprob == 20 || iprob == 40),
+                           "Selected iprob for diffusion pgen requires MHD enabled.")
   Real t0 = 0.5;
   Real diff_coeff = 0.0;
   Real amp = 1e-6;
-  // Get parameters for Gaussian profile
-  if (iprob == 10) {
-    diff_coeff = pin->GetReal("diffusion", "thermal_diff_coeff_code");
+  // Get common parameters for Gaussian profile
+  if ((iprob == 10) || (iprob == 30) || (iprob == 40)) {
     t0 = pin->GetOrAddReal("problem/diffusion", "t0", t0);
     amp = pin->GetOrAddReal("problem/diffusion", "amp", amp);
   }
+  // Heat diffusion of 1D Gaussian
+  if (iprob == 10) {
+    diff_coeff = pin->GetReal("diffusion", "thermal_diff_coeff_code");
+    // Viscous diffusion of 1D Gaussian
+  } else if (iprob == 30) {
+    diff_coeff = pin->GetReal("diffusion", "mom_diff_coeff_code");
+    // Ohmic diffusion of 1D Gaussian
+  } else if (iprob == 40) {
+    diff_coeff = pin->GetReal("diffusion", "ohm_diff_coeff_code");
+  }
 
   auto &coords = pmb->coords;
 
   pmb->par_for(
-      "ProblemGenerator: Diffusion Step", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      "ProblemGenerator: Diffusion", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int k, const int j, const int i) {
         u(IDN, k, j, i) = 1.0;
 
@@ -48,11 +65,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
         u(IM2, k, j, i) = 0.0;
         u(IM3, k, j, i) = 0.0;
 
-        u(IB1, k, j, i) = 0.0;
-        u(IB2, k, j, i) = 0.0;
-        u(IB3, k, j, i) = 0.0;
+        if (mhd_enabled) {
+          u(IB1, k, j, i) = 0.0;
+          u(IB2, k, j, i) = 0.0;
+          u(IB3, k, j, i) = 0.0;
+        }
 
-        Real eint;
+        Real eint = -1.0;
         // step function x1
         if (iprob == 0) {
           u(IB1, k, j, i) = Bx;
@@ -111,12 +130,31 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
           u(IB3, k, j, i) = y / r;
           u(IB1, k, j, i) = -x / r;
           eint = std::abs(r - 0.6) < 0.1 && std::abs(phi) < M_PI / 12.0 ? 12.0 : 10.0;
+          // Viscous diffusion of 1D Gaussian
+        } else if (iprob == 30) {
+          u(IM2, k, j, i) =
+              u(IDN, k, j, i) * amp /
+              std::pow(std::sqrt(4. * M_PI * diff_coeff * t0), 1.0) *
+              std::exp(-(std::pow(coords.Xc<1>(i), 2.)) / (4. * diff_coeff * t0));
+          eint = 1.0 / (gamma * (gamma - 1.0)); // c_s = 1 everywhere
+          // Ohmic diffusion of 1D Gaussian
+        } else if (iprob == 40) {
+          u(IB2, k, j, i) =
+              amp / std::pow(std::sqrt(4. * M_PI * diff_coeff * t0), 1.0) *
+              std::exp(-(std::pow(coords.Xc<1>(i), 2.)) / (4. * diff_coeff * t0));
+          eint = 1.0 / (gamma * (gamma - 1.0)); // c_s = 1 everywhere
         }
+
+        PARTHENON_REQUIRE(eint > 0.0, "Missing init of eint");
         u(IEN, k, j, i) =
             u(IDN, k, j, i) * eint +
-            0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + SQR(u(IB3, k, j, i)) +
-                   (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) /
-                       u(IDN, k, j, i));
+            0.5 * ((SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) /
+                   u(IDN, k, j, i));
+
+        if (mhd_enabled) {
+          u(IEN, k, j, i) +=
+              0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + SQR(u(IB3, k, j, i)));
+        }
       });
 }
 } // namespace diffusion
diff --git a/src/pgen/linear_wave_mhd.cpp b/src/pgen/linear_wave_mhd.cpp
index 67affbe7..6e2c105c 100644
--- a/src/pgen/linear_wave_mhd.cpp
+++ b/src/pgen/linear_wave_mhd.cpp
@@ -709,4 +709,38 @@ void Eigensystem(const Real d, const Real v1, const Real v2, const Real v3, cons
   left_eigenmatrix[6][6] = left_eigenmatrix[0][6];
 }
 
+// For decaying wave with diffusive processes test problem, dump max V_2
+Real HstMaxV2(MeshData<Real> *md) {
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+
+  const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+
+  IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
+  IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
+
+  Real max_v2 = 0.0;
+
+  Kokkos::parallel_reduce(
+      "HstMaxV2",
+      Kokkos::MDRangePolicy<Kokkos::Rank<4>>(
+          parthenon::DevExecSpace(), {0, kb.s, jb.s, ib.s},
+          {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1},
+          {1, 1, 1, ib.e + 1 - ib.s}),
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lmax) {
+        lmax = Kokkos::fmax(lmax, Kokkos::fabs(prim_pack(b, IV2, k, j, i)));
+      },
+      Kokkos::Max<Real>(max_v2));
+
+  return max_v2;
+}
+
+void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg) {
+  if (pin->GetOrAddBoolean("problem/linear_wave", "dump_max_v2", false)) {
+    auto hst_vars = pkg->Param<parthenon::HstVar_list>(parthenon::hist_param_key);
+    hst_vars.emplace_back(parthenon::HistoryOutputVar(
+        parthenon::UserHistoryOperation::max, HstMaxV2, "MaxAbsV2"));
+    pkg->UpdateParam(parthenon::hist_param_key, hst_vars);
+  }
+}
 } // namespace linear_wave_mhd
diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp
index b3e55aad..ddea80f3 100644
--- a/src/pgen/pgen.hpp
+++ b/src/pgen/pgen.hpp
@@ -25,6 +25,7 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin);
 void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
 void UserWorkAfterLoop(Mesh *mesh, parthenon::ParameterInput *pin,
                        parthenon::SimTime &tm);
+void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg);
 } // namespace linear_wave_mhd
 
 namespace cpaw {
diff --git a/tst/regression/CMakeLists.txt b/tst/regression/CMakeLists.txt
index c43b4de5..f4fa246e 100644
--- a/tst/regression/CMakeLists.txt
+++ b/tst/regression/CMakeLists.txt
@@ -49,6 +49,12 @@ setup_test_both("aniso_therm_cond_ring_multid" "--driver ${PROJECT_BINARY_DIR}/b
 setup_test_both("aniso_therm_cond_gauss_conv" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \
   --driver_input ${PROJECT_SOURCE_DIR}/inputs/diffusion.in --num_steps 24" "convergence")
 
+setup_test_both("diffusion" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \
+  --driver_input ${PROJECT_SOURCE_DIR}/inputs/diffusion.in --num_steps 12" "convergence")
+
+  setup_test_both("diffusion_linwave3d" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \
+  --driver_input ${PROJECT_SOURCE_DIR}/inputs/linear_wave3d.in --num_steps 2" "convergence")
+
 setup_test_both("field_loop" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \
   --driver_input ${PROJECT_SOURCE_DIR}/inputs/field_loop.in --num_steps 12" "convergence")
 
diff --git a/tst/regression/test_suites/aniso_therm_cond_ring_conv/aniso_therm_cond_ring_conv.py b/tst/regression/test_suites/aniso_therm_cond_ring_conv/aniso_therm_cond_ring_conv.py
index ba13129e..41a2dcc6 100644
--- a/tst/regression/test_suites/aniso_therm_cond_ring_conv/aniso_therm_cond_ring_conv.py
+++ b/tst/regression/test_suites/aniso_therm_cond_ring_conv/aniso_therm_cond_ring_conv.py
@@ -53,6 +53,10 @@ def Prepare(self, parameters, step):
             "parthenon/meshblock/nx3=1",
             "problem/diffusion/iprob=20",
             "parthenon/time/tlim=200.0",
+            # Work around for RKL2 integrator (that, by default, does not limit the
+            # timestep, which in newer versions of Parthenon results in triggering
+            # a fail-safe given the default init value of numeric_limits max.
+            "parthenon/time/dt_ceil=200.0",
             "parthenon/output0/dt=200.0",
             f"parthenon/output0/id={res}",
         ]
diff --git a/tst/regression/test_suites/aniso_therm_cond_ring_multid/aniso_therm_cond_ring_multid.py b/tst/regression/test_suites/aniso_therm_cond_ring_multid/aniso_therm_cond_ring_multid.py
index 3f933507..74083a63 100644
--- a/tst/regression/test_suites/aniso_therm_cond_ring_multid/aniso_therm_cond_ring_multid.py
+++ b/tst/regression/test_suites/aniso_therm_cond_ring_multid/aniso_therm_cond_ring_multid.py
@@ -87,6 +87,7 @@ def Prepare(self, parameters, step):
             "parthenon/time/tlim=200.0",
             "parthenon/output0/dt=200.0",
             f"parthenon/output0/id={step}",
+            "diffusion/integrator=unsplit",
         ]
 
         return parameters
diff --git a/tst/regression/test_suites/diffusion/__init__.py b/tst/regression/test_suites/diffusion/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tst/regression/test_suites/diffusion/diffusion.py b/tst/regression/test_suites/diffusion/diffusion.py
new file mode 100644
index 00000000..47a0a8b3
--- /dev/null
+++ b/tst/regression/test_suites/diffusion/diffusion.py
@@ -0,0 +1,217 @@
+# ========================================================================================
+# AthenaPK - a performance portable block structured AMR MHD code
+# Copyright (c) 2023-2024, Athena Parthenon Collaboration. All rights reserved.
+# Licensed under the 3-clause BSD License, see LICENSE file for details
+# ========================================================================================
+
+# Modules
+import math
+import numpy as np
+import matplotlib
+
+matplotlib.use("agg")
+import matplotlib.pylab as plt
+import sys
+import os
+import itertools
+import utils.test_case
+from scipy.optimize import curve_fit
+
+# To prevent littering up imported folders with .pyc files or __pycache_ folder
+sys.dont_write_bytecode = True
+
+diff_cfgs = ["visc", "ohm"]
+int_cfgs = ["unsplit", "rkl2"]
+res_cfgs = [256, 512, 1024]
+tlim = 2.0
+diff_coeff = 0.25
+
+all_cfgs = list(itertools.product(diff_cfgs, res_cfgs, int_cfgs))
+
+
+def get_outname(all_cfg):
+    diff, res, int_cfg = all_cfg
+    return f"{diff}_{res}_{int_cfg}"
+
+
+class TestCase(utils.test_case.TestCaseAbs):
+    def Prepare(self, parameters, step):
+        assert parameters.num_ranks <= 4, "Use <= 4 ranks for diffusion test."
+
+        diff, res, int_cfg = all_cfgs[step - 1]
+
+        outname = get_outname(all_cfgs[step - 1])
+
+        if diff == "visc":
+            fluid_ = "euler"
+            iprob_ = 30
+            viscosity_ = "isotropic"
+            resistivity_ = "none"
+        elif diff == "ohm":
+            fluid_ = "glmmhd"
+            iprob_ = 40
+            viscosity_ = "none"
+            resistivity_ = "ohmic"
+
+        parameters.driver_cmd_line_args = [
+            "parthenon/mesh/nx1=%d" % res,
+            "parthenon/meshblock/nx1=64",
+            "parthenon/mesh/x1min=-6.0",
+            "parthenon/mesh/x1max=6.0",
+            "parthenon/mesh/ix1_bc=outflow",
+            "parthenon/mesh/ox1_bc=outflow",
+            "parthenon/mesh/nx2=1",
+            "parthenon/meshblock/nx2=1",
+            "parthenon/mesh/x2min=-1.0",
+            "parthenon/mesh/x2max=1.0",
+            "parthenon/mesh/nx3=1",
+            "parthenon/meshblock/nx3=1",
+            f"parthenon/output0/id={outname}",
+            f"parthenon/time/tlim={tlim}",
+            # Work around for RKL2 integrator (that, by default, does not limit the
+            # timestep, which in newer versions of Parthenon results in triggering
+            # a fail-safe given the default init value of numeric_limits max.
+            "parthenon/time/dt_ceil=%f" % tlim,
+            f"hydro/fluid={fluid_}",
+            "hydro/gamma=1.4",
+            "hydro/cfl=0.8",
+            "hydro/integrator=rk2",
+            f"problem/diffusion/iprob={iprob_}",
+            f"problem/diffusion/Bx=0.0",
+            f"problem/diffusion/By=0.0",
+            "diffusion/conduction=none",
+            f"diffusion/viscosity={viscosity_}",
+            f"diffusion/resistivity={resistivity_}",
+            # we can set both as their activity is controlled separately
+            f"diffusion/mom_diff_coeff_code={diff_coeff}",
+            f"diffusion/ohm_diff_coeff_code={diff_coeff}",
+            f"diffusion/integrator={int_cfg}",
+            f"diffusion/rkl2_max_dt_ratio=200",
+        ]
+
+        return parameters
+
+    def Analyse(self, parameters):
+        sys.path.insert(
+            1,
+            parameters.parthenon_path
+            + "/scripts/python/packages/parthenon_tools/parthenon_tools",
+        )
+
+        try:
+            import phdf
+        except ModuleNotFoundError:
+            print("Couldn't find module to read Parthenon hdf5 files.")
+            return False
+
+        tests_passed = True
+
+        def get_ref(x):
+            return (
+                1e-6
+                / np.sqrt(4.0 * np.pi * diff_coeff * (0.5 + tlim))
+                * np.exp(-(x**2.0) / (4.0 * diff_coeff * (0.5 + tlim)))
+            )
+
+        num_diff = len(diff_cfgs)
+        for idx_diff in range(num_diff):
+            num_rows = len(res_cfgs)
+            num_cols = len(int_cfgs)
+            fig, p = plt.subplots(num_rows + 1, 2, sharey="row", sharex="row")
+
+            l1_err = np.zeros((len(int_cfgs), len(res_cfgs)))
+            for step in range(len(all_cfgs)):
+                diff, res, int_cfg = all_cfgs[step]
+                # only plot results for this diffusion process
+                if idx_diff != diff_cfgs.index(diff):
+                    continue
+                outname = get_outname(all_cfgs[step])
+                data_filename = (
+                    f"{parameters.output_path}/parthenon.{outname}.final.phdf"
+                )
+                data_file = phdf.phdf(data_filename)
+                # Flatten=true (default) is currently (Sep 24) broken so we manually flatten
+                components = data_file.GetComponents(
+                    data_file.Info["ComponentNames"], flatten=False
+                )
+                zz, yy, xx = data_file.GetVolumeLocations()
+                mask = yy == yy[0]
+                if diff == "visc":
+                    var_name = "prim_velocity_2"
+                elif diff == "ohm":
+                    var_name = "prim_magnetic_field_2"
+                else:
+                    print("Unknon diffusion type to process test results!")
+                    return False
+
+                v2 = components[var_name].ravel()[mask]
+                x = xx[mask]
+                row = res_cfgs.index(res)
+                col = int_cfgs.index(int_cfg)
+
+                v2_ref = get_ref(x)
+                l1 = np.average(np.abs(v2 - v2_ref))
+                l1_err[
+                    int_cfgs.index(int_cfg),
+                    res_cfgs.index(res),
+                ] = l1
+                p[row, col].plot(x, v2, label=f"N={res} L$_1$={l1:.2g}")
+
+            # Plot convergence
+            for j, int_cfg in enumerate(int_cfgs):
+                p[0, j].set_title(f"Integrator: {int_cfg}")
+
+                p[-1, j].plot(
+                    res_cfgs,
+                    l1_err[j, :],
+                    label=f"data",
+                )
+
+                # Simple convergence estimator
+                conv_model = lambda log_n, log_a, conv_rate: conv_rate * log_n + log_a
+                popt, pconv = curve_fit(
+                    conv_model, np.log(res_cfgs), np.log(l1_err[j, :])
+                )
+                conv_a, conv_measured = popt
+                # Note that the RKL2 convergence on the plots is currently significantly better
+                # than expected (<-3) though the L1 errors themself are larger than the unsplit
+                # integrator (as expected).
+                # For a more reasonable test (which would take longer), reduce the RKL2 ratio to,
+                # say, 200 and extend the resolution grid to 1024 (as the first data point at N=128
+                # is comparatively worse than at N>128).
+                if conv_measured > -1.95:
+                    print(
+                        f"!!!\nConvergence for test with {int_cfg} integrator "
+                        f"is worse ({conv_measured}) than expected (-1.95).\n!!!"
+                    )
+                    tests_passed = False
+                p[-1, j].plot(
+                    res_cfgs,
+                    np.exp(conv_a) * res_cfgs**conv_measured,
+                    ":",
+                    lw=0.75,
+                    label=f"Measured conv: {conv_measured:.2f}",
+                )
+
+            p[-1, 0].set_xscale("log")
+            p[-1, 0].set_yscale("log")
+            p[-1, 0].legend(fontsize=6)
+            p[-1, 1].legend(fontsize=6)
+
+            # Plot reference lines
+            x = np.linspace(-6, 6, 400)
+            for i in range(num_rows):
+                for j in range(num_cols):
+                    y = get_ref(x)
+                    p[i, j].plot(x, y, "-", lw=0.5, color="black", alpha=0.8)
+                    p[i, j].grid()
+                    p[i, j].legend(fontsize=6)
+
+            fig.tight_layout()
+            fig.savefig(
+                os.path.join(parameters.output_path, f"{diff_cfgs[idx_diff]}.png"),
+                bbox_inches="tight",
+                dpi=300,
+            )
+
+        return tests_passed
diff --git a/tst/regression/test_suites/diffusion_linwave3d/__init__.py b/tst/regression/test_suites/diffusion_linwave3d/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tst/regression/test_suites/diffusion_linwave3d/diffusion_linwave3d.py b/tst/regression/test_suites/diffusion_linwave3d/diffusion_linwave3d.py
new file mode 100644
index 00000000..d684a55c
--- /dev/null
+++ b/tst/regression/test_suites/diffusion_linwave3d/diffusion_linwave3d.py
@@ -0,0 +1,183 @@
+# ========================================================================================
+# AthenaPK - a performance portable block structured AMR MHD code
+# Copyright (c) 2023-2024, Athena Parthenon Collaboration. All rights reserved.
+# Licensed under the 3-clause BSD License, see LICENSE file for details
+# ========================================================================================
+
+# Modules
+import math
+import numpy as np
+import matplotlib
+
+matplotlib.use("agg")
+import matplotlib.pylab as plt
+import sys
+import os
+import utils.test_case
+from numpy.polynomial import Polynomial
+
+""" To prevent littering up imported folders with .pyc files or __pycache_ folder"""
+sys.dont_write_bytecode = True
+
+# if this is updated make sure to update the assert statements for the number of MPI ranks, too
+lin_res = [16, 32]  # resolution for linear convergence
+
+# Upper bound on relative L1 error for each above nx1:
+error_rel_tols = [0.22, 0.05]
+# lower bound on convergence rate at final (Nx1=64) asymptotic convergence regime
+rate_tols = [2.0]  # convergence rate > 3.0 for this particular resolution, sovler
+
+method = "explicit"
+
+_nu = 0.01
+_kappa = _nu * 2.0
+_eta = _kappa
+_c_s = 0.5  # slow mode wave speed of AthenaPK linear wave configuration
+
+
+class TestCase(utils.test_case.TestCaseAbs):
+    def Prepare(self, parameters, step):
+        res = lin_res[(step - 1)]
+        # make sure we can evenly distribute the MeshBlock sizes
+        err_msg = "Num ranks must be multiples of 2 for test."
+        assert parameters.num_ranks == 1 or parameters.num_ranks % 2 == 0, err_msg
+        # ensure a minimum block size of 4
+        assert 2 * res / parameters.num_ranks >= 8, "Use <= 8 ranks for test."
+
+        mb_nx1 = (2 * res) // parameters.num_ranks
+        # ensure that nx1 is <= 128 when using scratch (V100 limit on test system)
+        while mb_nx1 > 128:
+            mb_nx1 //= 2
+
+        parameters.driver_cmd_line_args = [
+            f"parthenon/mesh/nx1={2 * res}",
+            f"parthenon/meshblock/nx1={mb_nx1}",
+            f"parthenon/mesh/nx2={res}",
+            f"parthenon/meshblock/nx2={res}",
+            f"parthenon/mesh/nx3={res}",
+            f"parthenon/meshblock/nx3={res}",
+            "parthenon/mesh/nghost=2",
+            "parthenon/time/integrator=vl2",
+            "parthenon/time/tlim=3.0",
+            "hydro/reconstruction=plm",
+            "hydro/fluid=glmmhd",
+            "hydro/riemann=hlld",
+            # enable history dump to track decay of v2 component
+            "parthenon/output2/file_type=hst",
+            "parthenon/output2/dt=0.03",
+            "problem/linear_wave/dump_max_v2=true",
+            f"parthenon/job/problem_id={res}",  # hack to rename parthenon.hst to res.hst
+            # setup linear wave (L slow mode)
+            "job/problem_id=linear_wave_mhd",
+            "problem/linear_wave/amp=1e-4",
+            "problem/linear_wave/wave_flag=2",
+            "problem/linear_wave/compute_error=false",  # done here below, not in the pgen
+            # setup diffusive processes
+            "diffusion/integrator=unsplit",
+            "diffusion/conduction=isotropic",
+            "diffusion/conduction_coeff=fixed",
+            f"diffusion/thermal_diff_coeff_code={_kappa}",
+            "diffusion/viscosity=isotropic",
+            "diffusion/viscosity_coeff=fixed",
+            f"diffusion/mom_diff_coeff_code={_nu}",
+            "diffusion/resistivity=ohmic",
+            "diffusion/resistivity_coeff=fixed",
+            f"diffusion/ohm_diff_coeff_code={_eta}",
+        ]
+
+        return parameters
+
+    def Analyse(self, parameters):
+        analyze_status = True
+
+        # Following test evaluation is adapted from the one in Athena++.
+        # This also includes the limits/tolerances set above, which are identical to Athena++.
+
+        # Lambda=1 for Athena++'s linear wave setups in 1D, 2D, and 3D:
+        L = 1.0
+        ksqr = (2.0 * np.pi / L) ** 2
+        # Equation 3.13 from Ryu, et al. (modified to add thermal conduction term)
+        # fast mode decay rate = (19\nu/4 + 3\eta + 3(\gamma-1)^2*kappa/gamma/4)*(2/15)*k^2
+        # Equation 3.14 from Ryu, et al. (modified to add thermal conduction term)
+        # slow mode decay rate = (4\nu + 3\eta/4 + 3(\gamma-1)^2*kappa/gamma)*(2/15)*k^2
+        slow_mode_rate = (
+            (4.0 * _nu + 3.0 * _eta / 4.0 + _kappa * 4.0 / 5.0) * (2.0 / 15.0) * ksqr
+        )
+
+        # Equation 3.16
+        re_num = (4.0 * np.pi**2 * _c_s) / (L * slow_mode_rate)
+        analyze_status = True
+        errors_abs = []
+
+        for nx, err_tol in zip(lin_res, error_rel_tols):
+            print(
+                "[Decaying 3D Linear Wave]: "
+                "Mesh size {} x {} x {}".format(2 * nx, nx, nx)
+            )
+            filename = os.path.join(parameters.output_path, f"{nx}.out2.hst")
+            hst_data = np.genfromtxt(filename, names=True, skip_header=1)
+
+            tt = hst_data["1time"]
+            max_vy = hst_data["13MaxAbsV2"]
+            # estimate the decay rate from simulation, using weighted least-squares (WLS)
+            yy = np.log(np.abs(max_vy))
+            plt.plot(tt, yy)
+            plt.show()
+            p, [resid, rank, sv, rcond] = Polynomial.fit(
+                tt, yy, 1, w=np.sqrt(max_vy), full=True
+            )
+            resid_normal = np.sum((yy - p(tt)) ** 2)
+            r2 = 1 - resid_normal / (yy.size * yy.var())
+            pnormal = p.convert(domain=(-1, 1))
+            fit_rate = -pnormal.coef[-1]
+
+            error_abs = np.fabs(slow_mode_rate - fit_rate)
+            errors_abs += [error_abs]
+            error_rel = np.fabs(slow_mode_rate / fit_rate - 1.0)
+            err_rel_tol_percent = err_tol * 100.0
+
+            print(
+                "[Decaying 3D Linear Wave {}]: Reynolds number of slow mode: {}".format(
+                    method, re_num
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: R-squared of WLS regression = {}".format(
+                    method, r2
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Analytic decay rate = {}".format(
+                    method, slow_mode_rate
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Measured decay rate = {}".format(
+                    method, fit_rate
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Decay rate absolute error = {}".format(
+                    method, error_abs
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Decay rate relative error = {}".format(
+                    method, error_rel
+                )
+            )
+
+            if error_rel > err_tol:
+                print(
+                    "WARN [Decaying 3D Linear Wave {}]: decay rate disagrees"
+                    " with prediction by >{}%".format(method, err_rel_tol_percent)
+                )
+                analyze_status = False
+            else:
+                print(
+                    "[Decaying 3D Linear Wave {}]: decay rate is within "
+                    "{}% of analytic value".format(method, err_rel_tol_percent)
+                )
+                print("")
+
+        return analyze_status