From 3daf5c9108261ac4c168c253d5fdc345f4c03bb1 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 21 Aug 2024 21:34:41 +0200 Subject: [PATCH] Add isotropic thermal conduction and RKL2 supertimestepping (#1) * 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 1d0cb198c92b9eee579d3ecb5617d0629b1f9941. * Revert "Add perturb to cloud pgen" This reverts commit 6df018fe852fe4c285ac5929ed61d5fc67b17ddf. * 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 changelog * Address comments --- .github/workflows/ci.yml | 1 + CHANGELOG.md | 1 + README.md | 3 +- docs/input.md | 50 +- inputs/diffusion.in | 8 +- src/CMakeLists.txt | 1 + src/hydro/diffusion/conduction.cpp | 468 +++++++++++++----- src/hydro/diffusion/diffusion.cpp | 31 ++ src/hydro/diffusion/diffusion.hpp | 26 +- src/hydro/hydro.cpp | 147 ++++-- src/hydro/hydro.hpp | 1 + src/hydro/hydro_driver.cpp | 361 +++++++++++++- src/main.cpp | 1 + src/main.hpp | 4 +- src/pgen/cloud.cpp | 17 +- src/pgen/diffusion.cpp | 22 +- src/units.hpp | 2 + tst/regression/CMakeLists.txt | 3 + .../aniso_therm_cond_gauss_conv/__init__.py | 0 .../aniso_therm_cond_gauss_conv.py | 221 +++++++++ .../aniso_therm_cond_ring_conv/README.md | 2 +- 21 files changed, 1168 insertions(+), 202 deletions(-) create mode 100644 src/hydro/diffusion/diffusion.cpp create mode 100644 tst/regression/test_suites/aniso_therm_cond_gauss_conv/__init__.py create mode 100644 tst/regression/test_suites/aniso_therm_cond_gauss_conv/aniso_therm_cond_gauss_conv.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e24e197b..05049ba7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -53,6 +53,7 @@ jobs: build/tst/regression/outputs/cluster_hse/analytic_comparison.png build/tst/regression/outputs/cluster_tabular_cooling/convergence.png build/tst/regression/outputs/aniso_therm_cond_ring_conv/ring_convergence.png + build/tst/regression/outputs/aniso_therm_cond_gauss_conv/cond.png build/tst/regression/outputs/field_loop/field_loop.png build/tst/regression/outputs/riemann_hydro/shock_tube.png build/tst/regression/outputs/turbulence/parthenon.hst diff --git a/CHANGELOG.md b/CHANGELOG.md index 83bc715d..47b1e93d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop (i.e., `main` branch) ### Added (new features/APIs/variables/...) +- [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping ### Changed (changing behavior/API/variables/...) - [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15) diff --git a/README.md b/README.md index 04881cd0..13d7280a 100644 --- a/README.md +++ b/README.md @@ -14,7 +14,8 @@ Current features include - HLLE (hydro and MHD), HLLC (hydro), and HLLD (MHD) Riemann solvers - adiabatic equation of state - MHD based on hyperbolic divergence cleaning following Dedner+ 2002 - - anisotropic thermal conduction + - isotropic and anisotropic thermal conduction + - operator-split, second-order RKL2 supertimestepping for diffusive terms - optically thin cooling based on tabulated cooling tables with either Townsend 2009 exact integration or operator-split subcycling - static and adaptive mesh refinement - problem generators for diff --git a/docs/input.md b/docs/input.md index f7023463..11dae8c5 100644 --- a/docs/input.md +++ b/docs/input.md @@ -69,15 +69,22 @@ conserved to primitive conversion if both are defined. #### Diffusive processes -##### Anisotropic thermal conduction (required MHD) +##### Isotropic (hydro and MHD) and anisotropic thermal conduction (only MHD) In the presence of magnetic fields thermal conduction is becoming anisotropic with the flux along the local magnetic field direction typically being much stronger than the flux perpendicular to the magnetic field. From a theoretical point of view, thermal conduction is included in the system of MHD equations by an additional term in the total energy equation: ```math -\delta_t E + \nabla \cdot (... + \mathbf{F}) \quad \mathrm{with}\\ -\mathbf{F} = - \kappa \mathbf{\hat b} (\mathbf{\hat b \cdot \nabla T}) +\delta_t E + \nabla \cdot (... + \mathbf{F}_\mathrm{c}) +``` +where the full thermal conduction flux $`\mathbf{F}_\mathrm{c}`$ contains both the classic thermal conduction +```math +\mathbf{F}_\mathrm{classic} = - \kappa \mathbf{\hat b} (\mathbf{\hat b \cdot \nabla T}) +``` +as well as the saturated flux (as introduced by [^CM77]) +```math +\mathbf{F}_\mathrm{sat} = - 5 \phi \rho^{-1/2} p^{3/2} \mathrm{sgn}(\mathbf{\hat b \cdot \nabla T}) \mathbf{\hat b} ``` From an implementation point of view, two options implemented and can be configured within a `` block in the input file. @@ -86,23 +93,52 @@ the integration step (before flux correction in case of AMR, and calculating the Moreover, they are implemented explicitly, i.e., they add a (potentially very restrictive) constraint to the timestep due to the scaling with $`\propto \Delta_x^2`$. Finally, we employ limiters for calculating the temperature gradients following Sharma & Hammett (2007)[^SH07]. This prevents unphysical conduction against the gradient, which may be introduced because the off-axis gradients are not centered on the interfaces. +Similarly, to account for the different nature of classic and saturated fluxes (parabolic and hyperbolic, respectively), +we follow [^M+12] and use a smooth transition +```math +\mathbf{F}_\mathrm{c} = \frac{q}{q + F_\mathrm{classic}} \mathbf{F}_\mathrm{classic} \quad \mathrm{with} \quad q = 5 \phi \rho^{-1/2} p^{3/2} +``` +and upwinding of the hyperbolic, saturated fluxes. -To enable conduction, set +To enable thermal conduction, set Parameter: `conduction` (string) - `none` : No thermal conduction +- `isotropic` : Isotropic thermal conduction +- `anisotropic` : Anisotropic thermal conduction + +In addition the coefficient (or diffusivity) needs to be set + +Parameter: `conduction_coeff` (string) - `spitzer` : Anisotropic thermal conduction with a temperature dependent classic Spitzer thermal conductivity $`\kappa (T) = c_\kappa T^{5/2} \mathrm{erg/s/K/cm}`$ and - $`c_\kappa`$ being constant prefactor (set via `diffusion/spitzer_cond_in_erg_by_s_K_cm` with a default value of $`4.6\times10^{-7}`$). Note, as indicated by the units in the input parameter name, this kind of thermal conductivity requires a full set of units + $`c_\kappa`$ being constant prefactor (set via the additional `diffusion/spitzer_cond_in_erg_by_s_K_cm` parameter with a default value of $`4.6\times10^{-7}`$ which assumes a fully ionized hydrogen plasma [^CM77] with $`\ln \lambda = 40`$ approximating ICM conditions). Note, as indicated by the units in the input parameter name, this kind of thermal conductivity requires a full set of units to be defined for the simulation. -- `thermal_diff` : Contrary to a temperature dependent conductivity, a simple thermal diffusivity can be used instead for which +- `fixed` : Contrary to a temperature dependent conductivity, a simple thermal diffusivity can be used instead for which the conduction flux is $`\mathbf{F} = - \chi \rho \mathbf{\hat b} (\mathbf{\hat b \cdot \nabla \frac{p_\mathrm{th}}{\rho}})`$ -Here, the strength, $`\chi`$, is controlled via the `thermal_diff_coeff_code` parameter in code units. +Here, the strength, $`\chi`$, is controlled via the additional `thermal_diff_coeff_code` parameter in code units. Given the dimensions of $`L^2/T`$ it is referred to a thermal diffusivity rather than thermal conductivity. +Parameter: `conduction_sat_phi` (float) +- Default value 0.3\ +Factor to account for the uncertainty in the estimated of saturated fluxes, see [^CM77]. +Default value corresponds to the typical value used in literature and goes back to [^MMM80] and [^BM82]. + + [^SH07]: P. Sharma and G. W. Hammett, "Preserving monotonicity in anisotropic diffusion," Journal of Computational Physics, vol. 227, no. 1, Art. no. 1, 2007, doi: https://doi.org/10.1016/j.jcp.2007.07.026. +[^M+12]: + A. Mignone, C. Zanni, P. Tzeferacos, B. van Straalen, P. Colella, and G. Bodo, “THE PLUTO CODE FOR ADAPTIVE MESH COMPUTATIONS IN ASTROPHYSICAL FLUID DYNAMICS,” The Astrophysical Journal Supplement Series, vol. 198, Art. no. 1, Dec. 2011, doi: https://doi.org/10.1088/0067-0049/198/1/7 + +[^CM77]: + L. Cowie and C. F. McKee, “The evaporation of spherical clouds in a hot gas. I. Classical and saturated mass loss rates.,” , vol. 211, pp. 135–146, Jan. 1977, doi: https://doi.org/10.1086/154911 + +[^MMM80]: + C. E. Max, C. F. McKee, and W. C. Mead, “A model for laser driven ablative implosions,” The Physics of Fluids, vol. 23, Art. no. 8, 1980, doi: https://doi.org/10.1063/1.863183 + +[^BM82]: + S. A. Balbus and C. F. McKee, “The evaporation of spherical clouds in a hot gas. III - Suprathermal evaporation,” , vol. 252, pp. 529–552, Jan. 1982, doi: https://doi.org/10.1086/159581 ### Additional MHD options in `` block diff --git a/inputs/diffusion.in b/inputs/diffusion.in index a1a3f723..44d02bcb 100644 --- a/inputs/diffusion.in +++ b/inputs/diffusion.in @@ -14,7 +14,8 @@ Bx = 1.0 # Bx for x1 step function (permutated for iprobs in other direction By = 0.0 # By for x1 step function (permutated for iprobs in other directions) #iprob = 10 # Diffusion of Gaussian profile in x1 direction -sigma = 0.1 # standard deviation of Gaussian for iprob=10 +t0 = 0.5 # Temporal offset for initial Gaussian profile +amp = 1e-6 # Amplitude of Gaussian profile iprob = 20 # ring diffusion in x1-x2 plane; 21 for x2-x3 plane; 22 for x3-x1 plane @@ -59,8 +60,11 @@ reconstruction = dc gamma = 2.0 -conduction = thermal_diff +integrator = unsplit +conduction = anisotropic +conduction_coeff = fixed thermal_diff_coeff_code = 0.01 +rkl2_max_dt_ratio = 400.0 file_type = hdf5 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5d8fd375..ddcf09b7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,7 @@ add_executable( eos/adiabatic_glmmhd.cpp units.hpp eos/adiabatic_hydro.cpp + hydro/diffusion/diffusion.cpp hydro/diffusion/diffusion.hpp hydro/diffusion/conduction.cpp hydro/hydro_driver.cpp diff --git a/src/hydro/diffusion/conduction.cpp b/src/hydro/diffusion/conduction.cpp index 1e2bf0d6..d6e57d30 100644 --- a/src/hydro/diffusion/conduction.cpp +++ b/src/hydro/diffusion/conduction.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 //======================================================================================== // Athena++ astrophysical MHD code @@ -12,31 +12,29 @@ //! \brief // Parthenon headers +#include #include // AthenaPK headers #include "../../main.hpp" +#include "config.hpp" #include "diffusion.hpp" +#include "utils/error_checking.hpp" using namespace parthenon::package::prelude; +// Calculate the thermal *diffusivity*, \chi, in code units as the energy flux itself +// is calculated from -\chi \rho \nabla (p/\rho). KOKKOS_INLINE_FUNCTION -Real ThermalDiffusivity::Get(const Real pres, const Real rho, const Real gradTmag) const { - if (conduction_ == Conduction::thermal_diff) { +Real ThermalDiffusivity::Get(const Real pres, const Real rho) const { + if (conduction_coeff_type_ == ConductionCoeff::fixed) { return coeff_; - } else if (conduction_ == Conduction::spitzer) { - const Real T = mbar_over_kb_ * pres / rho; - const Real kappa = coeff_ * std::pow(T, 5. / 2.); // Full spitzer - const Real chi_spitzer = kappa * mbar_over_kb_ / rho; - - // Saturated total flux: fac * \rho * c_{s,isoth}^3 - // In practice: fac * \rho * c_{s,isoth}^3 * (gradT / gradTmag) - // where T is calculated based on p/rho in the code. - // Thus, everything is in code units and no conversion is required. - // The \rho above is cancelled as we convert the condution above to a diffusvity here. - const Real chi_sat = - 0.34 * std::pow(pres / rho, 3.0 / 2.0) / (gradTmag + TINY_NUMBER); - return std::min(chi_spitzer, chi_sat); + } else if (conduction_coeff_type_ == ConductionCoeff::spitzer) { + const Real T_cgs = mbar_ / kb_ * pres / rho; + const Real kappa_spitzer = coeff_ * std::pow(T_cgs, 5. / 2.); // Full spitzer + + // Convert conductivity to diffusivity + return kappa_spitzer * mbar_ / kb_ / rho; } else { return 0.0; @@ -64,74 +62,207 @@ Real EstimateConductionTimestep(MeshData *md) { const auto gm1 = hydro_pkg->Param("AdiabaticIndex"); const auto &thermal_diff = hydro_pkg->Param("thermal_diff"); + const auto &flux_sat_prefac = hydro_pkg->Param("conduction_sat_prefac"); + + if (thermal_diff.GetType() == Conduction::isotropic && + thermal_diff.GetCoeffType() == ConductionCoeff::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 thermal_diff_coeff = thermal_diff.Get(0.0, 0.0); + Kokkos::parallel_reduce( + "EstimateConductionTimestep (iso fixed)", + Kokkos::MDRangePolicy>( + 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)) / (thermal_diff_coeff + TINY_NUMBER)); + if (ndim >= 2) { + min_dt = fmin(min_dt, SQR(coords.Dxc<2>(k, j, i)) / + (thermal_diff_coeff + TINY_NUMBER)); + } + if (ndim >= 3) { + min_dt = fmin(min_dt, SQR(coords.Dxc<3>(k, j, i)) / + (thermal_diff_coeff + TINY_NUMBER)); + } + }, + Kokkos::Min(min_dt_cond)); + } else { + Kokkos::parallel_reduce( + "EstimateConductionTimestep (general)", + Kokkos::MDRangePolicy>( + 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); + const auto &prim = prim_pack(b); + const auto &rho = prim(IDN, k, j, i); + const auto &p = prim(IPR, k, j, i); + + const auto dTdx = 0.5 * + (prim(IPR, k, j, i + 1) / prim(IDN, k, j, i + 1) - + prim(IPR, k, j, i - 1) / prim(IDN, k, j, i - 1)) / + coords.Dxc<1>(i); + + const auto dTdy = ndim >= 2 + ? 0.5 * + (prim(IPR, k, j + 1, i) / prim(IDN, k, j + 1, i) - + prim(IPR, k, j - 1, i) / prim(IDN, k, j - 1, i)) / + coords.Dxc<2>(j) + : 0.0; + + const auto dTdz = ndim >= 3 + ? 0.5 * + (prim(IPR, k + 1, j, i) / prim(IDN, k + 1, j, i) - + prim(IPR, k - 1, j, i) / prim(IDN, k - 1, j, i)) / + coords.Dxc<3>(k) + : 0.0; + const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); + + // No temperature gradient -> no thermal conduction-> no timestep restriction + if (gradTmag == 0.0) { + return; + } + auto thermal_diff_coeff = thermal_diff.Get(p, rho); + + if (thermal_diff.GetType() == Conduction::isotropic) { + min_dt = fmin(min_dt, SQR(coords.Dxc<1>(k, j, i)) / thermal_diff_coeff); + if (ndim >= 2) { + min_dt = fmin(min_dt, SQR(coords.Dxc<2>(k, j, i)) / thermal_diff_coeff); + } + if (ndim >= 3) { + min_dt = fmin(min_dt, SQR(coords.Dxc<3>(k, j, i)) / thermal_diff_coeff); + } + return; + } + const auto &Bx = prim(IB1, k, j, i); + const auto &By = prim(IB2, k, j, i); + const auto &Bz = prim(IB3, k, j, i); + const auto Bmag = sqrt(SQR(Bx) + SQR(By) + SQR(Bz)); + // Need to have some local field for anisotropic conduction + if (Bmag == 0.0) { + return; + } + + // In the saturated regime, i.e., when the ratio of classic to saturated fluxes + // is large, the equation becomes hyperbolic with the signal speed of the + // conduction front being comparable to the sound speed, see [Balsara, Tilley, + // and Howk MANRAS 2008]. Therefore, we don't need to contrain the "parabolic" + // timestep here (and the hyperbolic one is constrained automatically by the + // fluid EstimateTimestep call). + auto const flux_sat = flux_sat_prefac * std::sqrt(p / rho) * p; + auto const flux_classic = thermal_diff_coeff * rho * gradTmag; + if (flux_classic / flux_sat > 100.) { + return; + } + + const auto costheta = + fabs(Bx * dTdx + By * dTdy + Bz * dTdz) / (Bmag * gradTmag); + + min_dt = fmin(min_dt, SQR(coords.Dxc<1>(k, j, i)) / + (thermal_diff_coeff * fabs(Bx) / Bmag * costheta + + TINY_NUMBER)); + if (ndim >= 2) { + min_dt = fmin(min_dt, SQR(coords.Dxc<2>(k, j, i)) / + (thermal_diff_coeff * fabs(By) / Bmag * costheta + + TINY_NUMBER)); + } + if (ndim >= 3) { + min_dt = fmin(min_dt, SQR(coords.Dxc<3>(k, j, i)) / + (thermal_diff_coeff * fabs(Bz) / Bmag * costheta + + TINY_NUMBER)); + } + }, + Kokkos::Min(min_dt_cond)); + } + + return fac * min_dt_cond; +} + +//--------------------------------------------------------------------------------------- +//! Calculate isotropic thermal conduction with fixed coefficient + +void ThermalFluxIsoFixed(MeshData *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 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{"prim"}); + + const int ndim = pmb->pmy_mesh->ndim; + + const auto &thermal_diff = hydro_pkg->Param("thermal_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 thermal_diff_coeff = thermal_diff.Get(0.0, 0.0); - Kokkos::parallel_reduce( - "EstimateConductionTimestep", - Kokkos::MDRangePolicy>( - 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) { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Thermal conduction X1 fluxes (iso)", + parthenon::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); - const auto &rho = prim(IDN, k, j, i); - const auto &p = prim(IPR, k, j, i); - // TODO(pgrete) when we introduce isotropic thermal conduction a lot of the - // following machinery should be hidden behind conditionals - const auto &Bx = prim(IB1, k, j, i); - const auto &By = prim(IB2, k, j, i); - const auto &Bz = prim(IB3, k, j, i); - const auto Bmag = sqrt(SQR(Bx) + SQR(By) + SQR(Bz)); - - const auto dTdx = 0.5 * - (prim(IPR, k, j, i + 1) / prim(IDN, k, j, i + 1) - - prim(IPR, k, j, i - 1) / prim(IDN, k, j, i - 1)) / - coords.Dxc<1>(i); - - const auto dTdy = 0.5 * - (prim(IPR, k, j + 1, i) / prim(IDN, k, j + 1, i) - - prim(IPR, k, j - 1, i) / prim(IDN, k, j - 1, i)) / - coords.Dxc<2>(j); - - const auto dTdz = ndim >= 3 - ? 0.5 * - (prim(IPR, k + 1, j, i) / prim(IDN, k + 1, j, i) - - prim(IPR, k - 1, j, i) / prim(IDN, k - 1, j, i)) / - coords.Dxc<3>(k) - : 0.0; - const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); - auto thermal_diff_coeff = thermal_diff.Get(p, rho, gradTmag); + const auto T_i = prim(IPR, k, j, i) / prim(IDN, k, j, i); + const auto T_im1 = prim(IPR, k, j, i - 1) / prim(IDN, k, j, i - 1); + const auto dTdx = (T_i - T_im1) / coords.Dxc<1>(k, j, i); + const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j, i - 1)); + cons.flux(X1DIR, IEN, k, j, i) -= thermal_diff_coeff * denf * dTdx; + }); - const auto denom = Bmag * gradTmag; - // if either Bmag or gradTmag are 0, no anisotropic thermal conduction - if (denom == 0.0) { - return; - } - const auto costheta = fabs(Bx * dTdx + By * dTdy + Bz * dTdz) / denom; - - min_dt = fmin( - min_dt, SQR(coords.Dxc<1>(k, j, i)) / - (thermal_diff_coeff * fabs(Bx) / Bmag * costheta + TINY_NUMBER)); - if (ndim >= 2) { - min_dt = fmin(min_dt, SQR(coords.Dxc<2>(k, j, i)) / - (thermal_diff_coeff * fabs(By) / Bmag * costheta + - TINY_NUMBER)); - } - if (ndim >= 3) { - min_dt = fmin(min_dt, SQR(coords.Dxc<3>(k, j, i)) / - (thermal_diff_coeff * fabs(Bz) / Bmag * costheta + - TINY_NUMBER)); - } - }, - Kokkos::Min(min_dt_cond)); + if (ndim < 2) { + return; + } + /* Compute heat fluxes in 2-direction --------------------------------------*/ + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Thermal conduction X2 fluxes (iso)", + 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); - return fac * min_dt_cond; + const auto T_j = prim(IPR, k, j, i) / prim(IDN, k, j, i); + const auto T_jm1 = prim(IPR, k, j - 1, i) / prim(IDN, k, j - 1, i); + const auto dTdy = (T_j - T_jm1) / coords.Dxc<2>(k, j, i); + const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j - 1, i)); + cons.flux(X2DIR, IEN, k, j, i) -= thermal_diff_coeff * denf * dTdy; + }); + /* Compute heat fluxes in 3-direction, 3D problem ONLY ---------------------*/ + if (ndim < 3) { + return; + } + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Thermal conduction X3 fluxes (iso)", + 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); + + const auto T_k = prim(IPR, k, j, i) / prim(IDN, k, j, i); + const auto T_km1 = prim(IPR, k - 1, j, i) / prim(IDN, k - 1, j, i); + const auto dTdz = (T_k - T_km1) / coords.Dxc<3>(k, j, i); + const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k - 1, j, i)); + cons.flux(X3DIR, IEN, k, j, i) -= thermal_diff_coeff * denf * dTdz; + }); } //--------------------------------------------------------------------------------------- -//! Calculate anisotropic thermal conduction +//! Calculate thermal conduction, general case, i.e., anisotropic and/or with varying +//! (incl. saturated) coefficient -void ThermalFluxAniso(MeshData *md) { +void ThermalFluxGeneral(MeshData *md) { auto pmb = md->GetBlockData(0)->GetBlockPointer(); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); @@ -146,18 +277,18 @@ void ThermalFluxAniso(MeshData *md) { const int ndim = pmb->pmy_mesh->ndim; const auto &thermal_diff = hydro_pkg->Param("thermal_diff"); + const auto &flux_sat_prefac = hydro_pkg->Param("conduction_sat_prefac"); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "Thermal conduction X1 fluxes", parthenon::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) { + DEFAULT_LOOP_PATTERN, "Thermal conduction X1 fluxes (general)", + parthenon::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); // Variables only required in 3D case Real dTdz = 0.0; - Real Bz = 0.0; // clang-format off /* Monotonized temperature difference dT/dy */ @@ -170,7 +301,7 @@ void ThermalFluxAniso(MeshData *md) { prim(IPR, k, j , i - 1) / prim(IDN, k, j , i - 1), prim(IPR, k, j , i - 1) / prim(IDN, k, j , i - 1) - prim(IPR, k, j - 1, i - 1) / prim(IDN, k, j - 1, i - 1)) / - coords.Dxc<2>(k, j, i); + coords.Dxc<2>( k, j, i); if (ndim >= 3) { /* Monotonized temperature difference dT/dz, 3D problem ONLY */ @@ -182,8 +313,7 @@ void ThermalFluxAniso(MeshData *md) { prim(IPR, k , j, i - 1) / prim(IDN, k , j, i - 1), prim(IPR, k , j, i - 1) / prim(IDN, k , j, i - 1) - prim(IPR, k - 1, j, i - 1) / prim(IDN, k - 1, j, i - 1)) / - coords.Dxc<3>(k, j, i); - Bz = 0.5 * (prim(IB3, k, j, i - 1) + prim(IB3, k, j, i)); + coords.Dxc<3>( k, j, i); } // clang-format on @@ -191,34 +321,69 @@ void ThermalFluxAniso(MeshData *md) { const auto T_im1 = prim(IPR, k, j, i - 1) / prim(IDN, k, j, i - 1); const auto dTdx = (T_i - T_im1) / coords.Dxc<1>(k, j, i); - // Calc interface values - const auto Bx = 0.5 * (prim(IB1, k, j, i - 1) + prim(IB1, k, j, i)); - const auto By = 0.5 * (prim(IB2, k, j, i - 1) + prim(IB2, k, j, i)); - auto B02 = SQR(Bx) + SQR(By) + SQR(Bz); - B02 = std::max(B02, TINY_NUMBER); /* limit in case B=0 */ - const auto bDotGradT = Bx * dTdx + By * dTdy + Bz * dTdz; - const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j, i - 1)); - const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); const auto thermal_diff_f = - 0.5 * - (thermal_diff.Get(prim(IPR, k, j, i), prim(IDN, k, j, i), gradTmag) + - thermal_diff.Get(prim(IPR, k, j, i - 1), prim(IDN, k, j, i - 1), gradTmag)); - cons.flux(X1DIR, IEN, k, j, i) -= thermal_diff_f * denf * (Bx * bDotGradT) / B02; + 0.5 * (thermal_diff.Get(prim(IPR, k, j, i), prim(IDN, k, j, i)) + + thermal_diff.Get(prim(IPR, k, j, i - 1), prim(IDN, k, j, i - 1))); + const auto gradTmag = std::sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); + + // Calculate "classic" fluxes + Real flux_classic = 0.0; + Real flux_classic_mag = 0.0; + if (thermal_diff.GetType() == Conduction::anisotropic) { + const auto Bx = 0.5 * (prim(IB1, k, j, i - 1) + prim(IB1, k, j, i)); + const auto By = 0.5 * (prim(IB2, k, j, i - 1) + prim(IB2, k, j, i)); + const auto Bz = + ndim >= 3 ? 0.5 * (prim(IB3, k, j, i - 1) + prim(IB3, k, j, i)) : 0.0; + auto Bmag = std::sqrt(SQR(Bx) + SQR(By) + SQR(Bz)); + Bmag = std::max(Bmag, TINY_NUMBER); /* limit in case B=0 */ + const auto bx = Bx / Bmag; // unit vector component + const auto bDotGradT = (Bx * dTdx + By * dTdy + Bz * dTdz) / Bmag; + flux_classic = -thermal_diff_f * denf * bDotGradT * bx; + flux_classic_mag = std::abs(thermal_diff_f * denf * bDotGradT); + } else if (thermal_diff.GetType() == Conduction::isotropic) { + flux_classic = -thermal_diff_f * denf * dTdx; + flux_classic_mag = thermal_diff_f * denf * gradTmag; + } else { + PARTHENON_FAIL("Unknown thermal diffusion flux."); + } + + // Calculate saturated fluxes using upwinding, see (A3) in Mignone+12. + // Note that we are not concerned about the sign of flux_sat here. The way it is + // calculated it's always positive because we use it in the harmonic mean with + // the flux_classic_mag below. The correct sign is eventually picked up again from + // flux_classic. + Real flux_sat; + // Use first order limiting for now. + if (flux_classic > 0.0) { + flux_sat = flux_sat_prefac * std::sqrt(prim(IPR, k, j, i - 1) / denf) * + prim(IPR, k, j, i - 1); + } else if (flux_classic < 0.0) { + flux_sat = + flux_sat_prefac * std::sqrt(prim(IPR, k, j, i) / denf) * prim(IPR, k, j, i); + } else { + const auto presf = 0.5 * (prim(IPR, k, j, i) + prim(IPR, k, j, i - 1)); + flux_sat = flux_sat_prefac * std::sqrt(presf / denf) * presf; + } + + cons.flux(X1DIR, IEN, k, j, i) += + (flux_sat / (flux_sat + flux_classic_mag)) * flux_classic; }); + if (ndim < 2) { + return; + } /* Compute heat fluxes in 2-direction --------------------------------------*/ parthenon::par_for( - DEFAULT_LOOP_PATTERN, "Thermal conduction X2 fluxes", 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) { + DEFAULT_LOOP_PATTERN, "Thermal conduction X2 fluxes (general)", + 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); // Variables only required in 3D case Real dTdz = 0.0; - Real Bz = 0.0; // clang-format off /* Monotonized temperature difference dT/dx */ @@ -245,7 +410,6 @@ void ThermalFluxAniso(MeshData *md) { prim(IPR, k - 1, j - 1, i) / prim(IDN, k - 1, j - 1, i)) / coords.Dxc<3>(k, j, i); - Bz = 0.5 * (prim(IB3, k, j - 1, i) + prim(IB3, k, j, i)); } // clang-format on @@ -253,20 +417,50 @@ void ThermalFluxAniso(MeshData *md) { const auto T_jm1 = prim(IPR, k, j - 1, i) / prim(IDN, k, j - 1, i); const auto dTdy = (T_j - T_jm1) / coords.Dxc<2>(k, j, i); - // Calc interface values - const auto Bx = 0.5 * (prim(IB1, k, j - 1, i) + prim(IB1, k, j, i)); - const auto By = 0.5 * (prim(IB2, k, j - 1, i) + prim(IB2, k, j, i)); - Real B02 = SQR(Bx) + SQR(By) + SQR(Bz); - B02 = std::max(B02, TINY_NUMBER); /* limit in case B=0 */ - const auto bDotGradT = Bx * dTdx + By * dTdy + Bz * dTdz; - const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k, j - 1, i)); const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); const auto thermal_diff_f = - 0.5 * - (thermal_diff.Get(prim(IPR, k, j, i), prim(IDN, k, j, i), gradTmag) + - thermal_diff.Get(prim(IPR, k, j - 1, i), prim(IDN, k, j - 1, i), gradTmag)); - cons.flux(X2DIR, IEN, k, j, i) -= thermal_diff_f * denf * (By * bDotGradT) / B02; + 0.5 * (thermal_diff.Get(prim(IPR, k, j, i), prim(IDN, k, j, i)) + + thermal_diff.Get(prim(IPR, k, j - 1, i), prim(IDN, k, j - 1, i))); + + // Calculate "classic" fluxes + Real flux_classic = 0.0; + Real flux_classic_mag = 0.0; + if (thermal_diff.GetType() == Conduction::anisotropic) { + const auto Bx = 0.5 * (prim(IB1, k, j - 1, i) + prim(IB1, k, j, i)); + const auto By = 0.5 * (prim(IB2, k, j - 1, i) + prim(IB2, k, j, i)); + const auto Bz = + ndim >= 3 ? 0.5 * (prim(IB3, k, j - 1, i) + prim(IB3, k, j, i)) : 0.0; + auto Bmag = std::sqrt(SQR(Bx) + SQR(By) + SQR(Bz)); + Bmag = std::max(Bmag, TINY_NUMBER); /* limit in case B=0 */ + const auto by = By / Bmag; // unit vector component + const auto bDotGradT = (Bx * dTdx + By * dTdy + Bz * dTdz) / Bmag; + flux_classic = -thermal_diff_f * denf * bDotGradT * by; + flux_classic_mag = std::abs(thermal_diff_f * denf * bDotGradT); + } else if (thermal_diff.GetType() == Conduction::isotropic) { + flux_classic = -thermal_diff_f * denf * dTdy; + flux_classic_mag = thermal_diff_f * denf * gradTmag; + } else { + PARTHENON_FAIL("Unknown thermal diffusion flux."); + } + + // Calculate saturated fluxes,see comment above. + Real flux_sat; + // Use first order limiting for now. + if (flux_classic > 0.0) { + flux_sat = flux_sat_prefac * std::sqrt(prim(IPR, k, j - 1, i) / denf) * + prim(IPR, k, j - 1, i); + } else if (flux_classic < 0.0) { + flux_sat = + flux_sat_prefac * std::sqrt(prim(IPR, k, j, i) / denf) * prim(IPR, k, j, i); + } else { + const auto presf = 0.5 * (prim(IPR, k, j, i) + prim(IPR, k, j - 1, i)); + flux_sat = flux_sat_prefac * std::sqrt(presf / denf) * presf; + } + + // Calc interface values + cons.flux(X2DIR, IEN, k, j, i) += + (flux_sat / (flux_sat + flux_classic_mag)) * flux_classic; }); /* Compute heat fluxes in 3-direction, 3D problem ONLY ---------------------*/ if (ndim < 3) { @@ -274,9 +468,9 @@ void ThermalFluxAniso(MeshData *md) { } parthenon::par_for( - DEFAULT_LOOP_PATTERN, "Thermal conduction X3 fluxes", 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) { + DEFAULT_LOOP_PATTERN, "Thermal conduction X3 fluxes (general)", + 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); @@ -311,20 +505,46 @@ void ThermalFluxAniso(MeshData *md) { const auto T_km1 = prim(IPR, k - 1, j, i) / prim(IDN, k - 1, j, i); const auto dTdz = (T_k - T_km1) / coords.Dxc<3>(k, j, i); - const auto Bx = 0.5 * (prim(IB1, k - 1, j, i) + prim(IB1, k, j, i)); - const auto By = 0.5 * (prim(IB2, k - 1, j, i) + prim(IB2, k, j, i)); - const auto Bz = 0.5 * (prim(IB3, k - 1, j, i) + prim(IB3, k, j, i)); - Real B02 = SQR(Bx) + SQR(By) + SQR(Bz); - B02 = std::max(B02, TINY_NUMBER); /* limit in case B=0 */ - const auto bDotGradT = Bx * dTdx + By * dTdy + Bz * dTdz; - const auto denf = 0.5 * (prim(IDN, k, j, i) + prim(IDN, k - 1, j, i)); const auto gradTmag = sqrt(SQR(dTdx) + SQR(dTdy) + SQR(dTdz)); const auto thermal_diff_f = - 0.5 * - (thermal_diff.Get(prim(IPR, k, j, i), prim(IDN, k, j, i), gradTmag) + - thermal_diff.Get(prim(IPR, k - 1, j, i), prim(IDN, k - 1, j, i), gradTmag)); + 0.5 * (thermal_diff.Get(prim(IPR, k, j, i), prim(IDN, k, j, i)) + + thermal_diff.Get(prim(IPR, k - 1, j, i), prim(IDN, k - 1, j, i))); + + // Calculate "classic" fluxes + Real flux_classic = 0.0; + Real flux_classic_mag = 0.0; + if (thermal_diff.GetType() == Conduction::anisotropic) { + const auto Bx = 0.5 * (prim(IB1, k - 1, j, i) + prim(IB1, k, j, i)); + const auto By = 0.5 * (prim(IB2, k - 1, j, i) + prim(IB2, k, j, i)); + const auto Bz = 0.5 * (prim(IB3, k - 1, j, i) + prim(IB3, k, j, i)); + auto Bmag = std::sqrt(SQR(Bx) + SQR(By) + SQR(Bz)); + Bmag = std::max(Bmag, TINY_NUMBER); /* limit in case B=0 */ + const auto bz = Bz / Bmag; // unit vector component + const auto bDotGradT = (Bx * dTdx + By * dTdy + Bz * dTdz) / Bmag; + flux_classic = -thermal_diff_f * denf * bDotGradT * bz; + flux_classic_mag = std::abs(thermal_diff_f * denf * bDotGradT); + } else if (thermal_diff.GetType() == Conduction::isotropic) { + flux_classic = -thermal_diff_f * denf * dTdz; + flux_classic_mag = thermal_diff_f * denf * gradTmag; + } else { + PARTHENON_FAIL("Unknown thermal diffusion flux."); + } + // Calculate saturated fluxes,see comment above. + Real flux_sat; + // Use first order limiting for now. + if (flux_classic > 0.0) { + flux_sat = flux_sat_prefac * std::sqrt(prim(IPR, k - 1, j, i) / denf) * + prim(IPR, k - 1, j, i); + } else if (flux_classic < 0.0) { + flux_sat = + flux_sat_prefac * std::sqrt(prim(IPR, k, j, i) / denf) * prim(IPR, k, j, i); + } else { + const auto presf = 0.5 * (prim(IPR, k, j, i) + prim(IPR, k - 1, j, i)); + flux_sat = flux_sat_prefac * std::sqrt(presf / denf) * presf; + } - cons.flux(X3DIR, IEN, k, j, i) -= thermal_diff_f * denf * (Bz * bDotGradT) / B02; + cons.flux(X3DIR, IEN, k, j, i) += + (flux_sat / (flux_sat + flux_classic_mag)) * flux_classic; }); } diff --git a/src/hydro/diffusion/diffusion.cpp b/src/hydro/diffusion/diffusion.cpp new file mode 100644 index 00000000..20617ca2 --- /dev/null +++ b/src/hydro/diffusion/diffusion.cpp @@ -0,0 +1,31 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file diffusion.cpp +//! \brief + +// Parthenon headers +#include + +// AthenaPK headers +#include "../../main.hpp" +#include "diffusion.hpp" + +using namespace parthenon::package::prelude; + +TaskStatus CalcDiffFluxes(StateDescriptor *hydro_pkg, MeshData *md) { + const auto &conduction = hydro_pkg->Param("conduction"); + if (conduction != Conduction::none) { + const auto &thermal_diff = hydro_pkg->Param("thermal_diff"); + + if (conduction == Conduction::isotropic && + thermal_diff.GetCoeffType() == ConductionCoeff::fixed) { + ThermalFluxIsoFixed(md); + } else { + ThermalFluxGeneral(md); + } + } + return TaskStatus::complete; +} diff --git a/src/hydro/diffusion/diffusion.hpp b/src/hydro/diffusion/diffusion.hpp index 273464d7..7d522902 100644 --- a/src/hydro/diffusion/diffusion.hpp +++ b/src/hydro/diffusion/diffusion.hpp @@ -69,23 +69,37 @@ KOKKOS_INLINE_FUNCTION Real lim4(const Real A, const Real B, const Real C, const struct ThermalDiffusivity { private: - Real mbar_over_kb_; + Real mbar_, me_, kb_; Conduction conduction_; + ConductionCoeff conduction_coeff_type_; // "free" coefficient/prefactor. Value depends on conduction is set in the constructor. Real coeff_; public: KOKKOS_INLINE_FUNCTION - ThermalDiffusivity(Conduction conduction, Real coeff, Real mbar_over_kb) - : coeff_(coeff), conduction_(conduction), mbar_over_kb_(mbar_over_kb) {} + ThermalDiffusivity(Conduction conduction, ConductionCoeff conduction_coeff_type, + Real coeff, Real mbar, Real me, Real kb) + : conduction_(conduction), conduction_coeff_type_(conduction_coeff_type), + coeff_(coeff), mbar_(mbar), me_(me), kb_(kb) {} KOKKOS_INLINE_FUNCTION - Real Get(const Real pres, const Real rho, const Real gradTmag) const; + Real Get(const Real pres, const Real rho) const; + + KOKKOS_INLINE_FUNCTION + Conduction GetType() const { return conduction_; } + + KOKKOS_INLINE_FUNCTION + ConductionCoeff GetCoeffType() const { return conduction_coeff_type_; } }; Real EstimateConductionTimestep(MeshData *md); -//! Calculate anisotropic thermal conduction -void ThermalFluxAniso(MeshData *md); +//! Calculate isotropic thermal conduction with fixed coefficient +void ThermalFluxIsoFixed(MeshData *md); +//! Calculate thermal conduction (general case incl. anisotropic and saturated) +void ThermalFluxGeneral(MeshData *md); + +// Calculate all diffusion fluxes, i.e., update the .flux views in md +TaskStatus CalcDiffFluxes(StateDescriptor *hydro_pkg, MeshData *md); #endif // HYDRO_DIFFUSION_DIFFUSION_HPP_ diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index e5402454..54a97b7f 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -55,6 +55,35 @@ parthenon::Packages_t ProcessPackages(std::unique_ptr &pin) { return packages; } +// Using this per cycle function to populate various variables in +// Params that require global reduction *and* need to be set/known when +// the task list is constructed (versus when the task list is being executed). +// TODO(next person touching this function): If more/separate feature are required +// please separate concerns. +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::rkl2) && + (hydro_pkg->Param("conduction") != Conduction::none)) { + auto dt_diff = std::numeric_limits::max(); + 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())); + } +#ifdef MPI_PARALLEL + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &dt_diff, 1, MPI_PARTHENON_REAL, + MPI_MIN, MPI_COMM_WORLD)); +#endif + hydro_pkg->UpdateParam("dt_diff", dt_diff); + const auto max_dt_ratio = hydro_pkg->Param("rkl2_max_dt_ratio"); + if (max_dt_ratio > 0.0 && tm.dt / dt_diff > max_dt_ratio) { + tm.dt = max_dt_ratio * dt_diff; + } + } +} + template Real HydroHst(MeshData *md) { auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); @@ -404,6 +433,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { pkg->AddParam<>("mu", mu); pkg->AddParam<>("mu_e", mu_e); pkg->AddParam<>("He_mass_fraction", He_mass_fraction); + pkg->AddParam<>("mbar", mu * units.atomic_mass_unit()); // Following convention in the astro community, we're using mh as unit for the mean // molecular weight pkg->AddParam<>("mbar_over_kb", mu * units.mh() / units.k_boltzmann()); @@ -444,36 +474,93 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto conduction = Conduction::none; auto conduction_str = pin->GetOrAddString("diffusion", "conduction", "none"); - if (conduction_str == "spitzer") { - if (!pkg->AllParams().hasKey("mbar_over_kb")) { - PARTHENON_FAIL("Spitzer thermal conduction requires units and gas composition. " - "Please set a 'units' block and the 'hydro/He_mass_fraction' in " - "the input file."); - } - conduction = Conduction::spitzer; - - Real spitzer_coeff = - pin->GetOrAddReal("diffusion", "spitzer_cond_in_erg_by_s_K_cm", 4.6e-7); - // Convert to code units. No temp conversion as [T_phys] = [T_code]. - auto units = pkg->Param("units"); - spitzer_coeff *= units.erg() / (units.s() * units.cm()); - - auto mbar_over_kb = pkg->Param("mbar_over_kb"); - auto thermal_diff = ThermalDiffusivity(conduction, spitzer_coeff, mbar_over_kb); - pkg->AddParam<>("thermal_diff", thermal_diff); - - } else if (conduction_str == "thermal_diff") { - conduction = Conduction::thermal_diff; - Real thermal_diff_coeff_code = pin->GetReal("diffusion", "thermal_diff_coeff_code"); - auto thermal_diff = ThermalDiffusivity(conduction, thermal_diff_coeff_code, 0.0); - pkg->AddParam<>("thermal_diff", thermal_diff); - + if (conduction_str == "isotropic") { + conduction = Conduction::isotropic; + } else if (conduction_str == "anisotropic") { + conduction = Conduction::anisotropic; } else if (conduction_str != "none") { PARTHENON_FAIL( - "AthenaPK unknown conduction method. Options are: spitzer, thermal_diff"); + "Unknown conduction method. Options are: none, isotropic, anisotropic"); + } + // If conduction is enabled, process supported coefficients + if (conduction != Conduction::none) { + auto conduction_coeff_str = + pin->GetOrAddString("diffusion", "conduction_coeff", "none"); + auto conduction_coeff = ConductionCoeff::none; + + // Saturated conduction factor to account for "uncertainty", see + // Cowie & McKee 77 and a value of 0.3 is typical chosen (though using "weak + // evidence", see Balbus & MacKee 1982 and Max, McKee, and Mead 1980). + const auto conduction_sat_phi = + pin->GetOrAddReal("diffusion", "conduction_sat_phi", 0.3); + Real conduction_sat_prefac = 0.0; + + if (conduction_coeff_str == "spitzer") { + if (!pkg->AllParams().hasKey("mbar")) { + PARTHENON_FAIL("Spitzer thermal conduction requires units and gas composition. " + "Please set a 'units' block and the 'hydro/He_mass_fraction' in " + "the input file."); + } + conduction_coeff = ConductionCoeff::spitzer; + + // Default value assume fully ionized hydrogen plasma with Coulomb logarithm of 40 + // to approximate ICM conditions, i.e., 1.84e-5/ln Lambda = 4.6e-7. + Real spitzer_coeff = + pin->GetOrAddReal("diffusion", "spitzer_cond_in_erg_by_s_K_cm", 4.6e-7); + // Convert to code units. No temp conversion as [T_phys] = [T_code]. + auto units = pkg->Param("units"); + spitzer_coeff *= units.erg() / (units.s() * units.cm()); + + const auto mbar = pkg->Param("mbar"); + auto thermal_diff = + ThermalDiffusivity(conduction, conduction_coeff, spitzer_coeff, mbar, + units.electron_mass(), units.k_boltzmann()); + pkg->AddParam<>("thermal_diff", thermal_diff); + + const auto mu = pkg->Param("mu"); + // 6.86 again assumes a fully ionized hydrogen plasma in agreement with + // the assumptions above (technically this means mu = 0.5) and can be derived + // from eq (7) in CM77 assuming T_e = T_i. + conduction_sat_prefac = 6.86 * std::sqrt(mu) * conduction_sat_phi; + + } else if (conduction_coeff_str == "fixed") { + conduction_coeff = ConductionCoeff::fixed; + Real thermal_diff_coeff_code = + pin->GetReal("diffusion", "thermal_diff_coeff_code"); + auto thermal_diff = ThermalDiffusivity(conduction, conduction_coeff, + thermal_diff_coeff_code, 0.0, 0.0, 0.0); + pkg->AddParam<>("thermal_diff", thermal_diff); + // 5.0 prefactor comes from eq (8) in Cowie & McKee 1977 + // https://doi.org/10.1086/154911 + conduction_sat_prefac = 5.0 * conduction_sat_phi; + + } else { + PARTHENON_FAIL("Thermal conduction is enabled but no coefficient is set. Please " + "set diffusion/conduction_coeff to either 'spitzer' or 'fixed'"); + } + PARTHENON_REQUIRE(conduction_sat_prefac != 0.0, + "Saturated thermal conduction prefactor uninitialized."); + pkg->AddParam<>("conduction_sat_prefac", conduction_sat_prefac); } pkg->AddParam<>("conduction", conduction); + auto diffint_str = pin->GetOrAddString("diffusion", "integrator", "none"); + auto diffint = DiffInt::none; + if (diffint_str == "unsplit") { + diffint = DiffInt::unsplit; + } else if (diffint_str == "rkl2") { + diffint = DiffInt::rkl2; + auto rkl2_dt_ratio = pin->GetOrAddReal("diffusion", "rkl2_max_dt_ratio", -1.0); + pkg->AddParam<>("rkl2_max_dt_ratio", rkl2_dt_ratio); + } else if (diffint_str != "none") { + PARTHENON_FAIL("AthenaPK unknown integration method for diffusion processes. " + "Options are: none, unsplit, rkl2"); + } + if (diffint != DiffInt::none) { + pkg->AddParam("dt_diff", 0.0, true); // diffusive timestep constraint + } + pkg->AddParam<>("diffint", diffint); + if (fluid == Fluid::euler) { AdiabaticHydroEOS eos(pfloor, dfloor, efloor, vceil, eceil, gamma); pkg->AddParam<>("eos", eos); @@ -702,7 +789,9 @@ Real EstimateTimestep(MeshData *md) { min_dt = std::min(min_dt, tabular_cooling.EstimateTimeStep(md)); } - if (hydro_pkg->Param("conduction") != Conduction::none) { + // For RKL2 STS, the diffusive timestep is calculated separately in the driver + if ((hydro_pkg->Param("diffint") == DiffInt::unsplit) && + (hydro_pkg->Param("conduction") != Conduction::none)) { min_dt = std::min(min_dt, EstimateConductionTimestep(md)); } @@ -943,9 +1032,9 @@ TaskStatus CalculateFluxes(std::shared_ptr> &md) { }); } - const auto &conduction = pkg->Param("conduction"); - if (conduction != Conduction::none) { - ThermalFluxAniso(md.get()); + const auto &diffint = pkg->Param("diffint"); + if (diffint == DiffInt::unsplit) { + CalcDiffFluxes(pkg.get(), md.get()); } return TaskStatus::complete; diff --git a/src/hydro/hydro.hpp b/src/hydro/hydro.hpp index 0967d12e..d54e96b1 100644 --- a/src/hydro/hydro.hpp +++ b/src/hydro/hydro.hpp @@ -16,6 +16,7 @@ using namespace parthenon::package::prelude; namespace Hydro { parthenon::Packages_t ProcessPackages(std::unique_ptr &pin); +void PreStepMeshUserWorkInLoop(Mesh *pmesh, ParameterInput *pin, parthenon::SimTime &tm); std::shared_ptr Initialize(ParameterInput *pin); template diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index e33b1b0d..36bdce91 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -1,6 +1,6 @@ //======================================================================================== // AthenaPK - a performance portable block structured AMR astrophysical MHD code. -// Copyright (c) 2020-2022, Athena-Parthenon Collaboration. All rights reserved. +// Copyright (c) 2020-2023, Athena-Parthenon Collaboration. All rights reserved. // Licensed under the BSD 3-Clause License (the "LICENSE"). //======================================================================================== @@ -19,6 +19,7 @@ #include "../eos/adiabatic_hydro.hpp" #include "../pgen/cluster/agn_triggering.hpp" #include "../pgen/cluster/magnetic_tower.hpp" +#include "diffusion/diffusion.hpp" #include "glmmhd/glmmhd.hpp" #include "hydro.hpp" #include "hydro_driver.hpp" @@ -76,6 +77,310 @@ TaskStatus CalculateGlobalMinDx(MeshData *md) { return TaskStatus::complete; } +// Sets all fluxes to 0 +TaskStatus ResetFluxes(MeshData *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); + + // In principle, we'd only need to pack Metadata::WithFluxes here, but + // choosing to mirror other use in the code so that the packs are already cached. + std::vector flags_ind({Metadata::Independent}); + auto cons_pack = md->PackVariablesAndFluxes(flags_ind); + + const int ndim = pmb->pmy_mesh->ndim; + // Using separate loops for each dim as the launch overhead should be hidden + // by enough work over the entire pack and it allows to not use any conditionals. + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ResetFluxes X1", parthenon::DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, 0, cons_pack.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, + ib.e + 1, + KOKKOS_LAMBDA(const int b, const int v, const int k, const int j, const int i) { + auto &cons = cons_pack(b); + cons.flux(X1DIR, v, k, j, i) = 0.0; + }); + + if (ndim < 2) { + return TaskStatus::complete; + } + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ResetFluxes X2", parthenon::DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, 0, cons_pack.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e + 1, + ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int v, const int k, const int j, const int i) { + auto &cons = cons_pack(b); + cons.flux(X2DIR, v, k, j, i) = 0.0; + }); + + if (ndim < 3) { + return TaskStatus::complete; + } + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ResetFluxes X3", parthenon::DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, 0, cons_pack.GetDim(4) - 1, kb.s, kb.e + 1, jb.s, jb.e, + ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int v, const int k, const int j, const int i) { + auto &cons = cons_pack(b); + cons.flux(X3DIR, v, k, j, i) = 0.0; + }); + return TaskStatus::complete; +} + +TaskStatus RKL2StepFirst(MeshData *md_Y0, MeshData *md_Yjm1, + MeshData *md_Yjm2, MeshData *md_MY0, const int s_rkl, + const Real tau) { + auto pmb = md_Y0->GetBlockData(0)->GetBlockPointer(); + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + + // Compute coefficients. Meyer+2014 eq. (18) + Real mu_tilde_1 = 4. / 3. / + (static_cast(s_rkl) * static_cast(s_rkl) + + static_cast(s_rkl) - 2.); + + // In principle, we'd only need to pack Metadata::WithFluxes here, but + // choosing to mirror other use in the code so that the packs are already cached. + std::vector flags_ind({Metadata::Independent}); + auto Y0 = md_Y0->PackVariablesAndFluxes(flags_ind); + auto Yjm1 = md_Yjm1->PackVariablesAndFluxes(flags_ind); + auto Yjm2 = md_Yjm2->PackVariablesAndFluxes(flags_ind); + auto MY0 = md_MY0->PackVariablesAndFluxes(flags_ind); + + const int ndim = pmb->pmy_mesh->ndim; + // Using separate loops for each dim as the launch overhead should be hidden + // by enough work over the entire pack and it allows to not use any conditionals. + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "RKL first step", parthenon::DevExecSpace(), 0, + Y0.GetDim(5) - 1, 0, Y0.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int v, const int k, const int j, const int i) { + Yjm1(b, v, k, j, i) = + Y0(b, v, k, j, i) + mu_tilde_1 * tau * MY0(b, v, k, j, i); // Y_1 + Yjm2(b, v, k, j, i) = Y0(b, v, k, j, i); // Y_0 + }); + + return TaskStatus::complete; +} + +TaskStatus RKL2StepOther(MeshData *md_Y0, MeshData *md_Yjm1, + MeshData *md_Yjm2, MeshData *md_MY0, const Real mu_j, + const Real nu_j, const Real mu_tilde_j, const Real gamma_tilde_j, + const Real tau) { + auto pmb = md_Y0->GetBlockData(0)->GetBlockPointer(); + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + + // In principle, we'd only need to pack Metadata::WithFluxes here, but + // choosing to mirror other use in the code so that the packs are already cached. + std::vector flags_ind({Metadata::Independent}); + auto Y0 = md_Y0->PackVariablesAndFluxes(flags_ind); + auto Yjm1 = md_Yjm1->PackVariablesAndFluxes(flags_ind); + auto Yjm2 = md_Yjm2->PackVariablesAndFluxes(flags_ind); + auto MY0 = md_MY0->PackVariablesAndFluxes(flags_ind); + + const int ndim = pmb->pmy_mesh->ndim; + // Using separate loops for each dim as the launch overhead should be hidden + // by enough work over the entire pack and it allows to not use any conditionals. + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "RKL other step", parthenon::DevExecSpace(), 0, + Y0.GetDim(5) - 1, 0, Y0.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int v, const int k, const int j, const int i) { + // First calc this step + const auto &coords = Yjm1.GetCoords(b); + const Real MYjm1 = + parthenon::Update::FluxDivHelper(v, k, j, i, ndim, coords, Yjm1(b)); + const Real Yj = mu_j * Yjm1(b, v, k, j, i) + nu_j * Yjm2(b, v, k, j, i) + + (1.0 - mu_j - nu_j) * Y0(b, v, k, j, i) + + mu_tilde_j * tau * MYjm1 + + gamma_tilde_j * tau * MY0(b, v, k, j, i); + // Then shuffle vars for next step + Yjm2(b, v, k, j, i) = Yjm1(b, v, k, j, i); + Yjm1(b, v, k, j, i) = Yj; + }); + + return TaskStatus::complete; +} + +// Assumes that prim and cons are in sync initially. +// Guarantees that prim and cons are in sync at the end. +void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks, + const Real tau) { + + auto hydro_pkg = blocks[0]->packages.Get("Hydro"); + auto mindt_diff = hydro_pkg->Param("dt_diff"); + + // get number of RKL steps + // eq (21) using half hyperbolic timestep due to Strang split + int s_rkl = + static_cast(0.5 * (std::sqrt(9.0 + 16.0 * tau / mindt_diff) - 1.0)) + 1; + // ensure odd number of stages + if (s_rkl % 2 == 0) s_rkl += 1; + + if (parthenon::Globals::my_rank == 0) { + const auto ratio = 2.0 * tau / mindt_diff; + std::cout << "STS ratio: " << ratio << " Taking " << s_rkl << " steps." << std::endl; + if (ratio > 400.1) { + std::cout << "WARNING: ratio is > 400. Proceed at own risk." << std::endl; + } + } + + TaskID none(0); + + // Store initial u0 in u1 as "base" will continuously be updated but initial state Y0 is + // required for each stage. + TaskRegion ®ion_copy_out = ptask_coll->AddRegion(blocks.size()); + for (int i = 0; i < blocks.size(); i++) { + auto &tl = region_copy_out[i]; + auto &Y0 = blocks[i]->meshblock_data.Get("u1"); + auto &base = blocks[i]->meshblock_data.Get(); + tl.AddTask( + none, + [](MeshBlockData *dst, MeshBlockData *src) { + dst->Get("cons").data.DeepCopy(src->Get("cons").data); + dst->Get("prim").data.DeepCopy(src->Get("prim").data); + return TaskStatus::complete; + }, + Y0.get(), base.get()); + } + + TaskRegion ®ion_init = ptask_coll->AddRegion(blocks.size()); + for (int i = 0; i < blocks.size(); i++) { + auto &pmb = blocks[i]; + auto &tl = region_init[i]; + auto &base = pmb->meshblock_data.Get(); + + // Add extra registers. No-op for existing variables so it's safe to call every + // time. + // TODO(pgrete) this allocates all Variables, i.e., prim and cons vector, but only a + // subset is actually needed. Streamline to allocate only required vars. + pmb->meshblock_data.Add("MY0", base); + pmb->meshblock_data.Add("Yjm2", base); + } + + const int num_partitions = pmesh->DefaultNumPartitions(); + TaskRegion ®ion_calc_fluxes_step_init = ptask_coll->AddRegion(num_partitions); + for (int i = 0; i < num_partitions; i++) { + auto &tl = region_calc_fluxes_step_init[i]; + auto &base = pmesh->mesh_data.GetOrAdd("base", i); + const auto any = parthenon::BoundaryType::any; + auto start_bnd = tl.AddTask(none, parthenon::StartReceiveBoundBufs, base); + auto start_flxcor_recv = + tl.AddTask(none, parthenon::StartReceiveFluxCorrections, base); + + // Reset flux arrays (not guaranteed to be zero) + auto reset_fluxes = tl.AddTask(none, ResetFluxes, base.get()); + + // Calculate the diffusive fluxes for Y0 (here still "base" as nothing has been + // updated yet) so that we can store the result as MY0 and reuse later + // (in every subsetp). + auto hydro_diff_fluxes = + tl.AddTask(reset_fluxes, CalcDiffFluxes, hydro_pkg.get(), base.get()); + + auto send_flx = + tl.AddTask(hydro_diff_fluxes, parthenon::LoadAndSendFluxCorrections, base); + auto recv_flx = + tl.AddTask(start_flxcor_recv, parthenon::ReceiveFluxCorrections, base); + auto set_flx = + tl.AddTask(recv_flx | hydro_diff_fluxes, parthenon::SetFluxCorrections, base); + + auto &Y0 = pmesh->mesh_data.GetOrAdd("u1", i); + auto &MY0 = pmesh->mesh_data.GetOrAdd("MY0", i); + auto &Yjm2 = pmesh->mesh_data.GetOrAdd("Yjm2", i); + + auto init_MY0 = tl.AddTask(set_flx, parthenon::Update::FluxDivergence>, + base.get(), MY0.get()); + + // Initialize Y0 and Y1 and the recursion relation starting with j = 2 needs data from + // the two preceeding stages. + auto rkl2_step_first = tl.AddTask(init_MY0, RKL2StepFirst, Y0.get(), base.get(), + Yjm2.get(), MY0.get(), s_rkl, tau); + + // Update ghost cells of Y1 (as MY1 is calculated for each Y_j). + // Y1 stored in "base", see rkl2_step_first task. + // Update ghost cells (local and non local), prolongate and apply bound cond. + // TODO(someone) experiment with split (local/nonlocal) comms with respect to + // performance for various tests (static, amr, block sizes) and then decide on the + // best impl. Go with default call (split local/nonlocal) for now. + // TODO(pgrete) optimize (in parthenon) to only send subset of updated vars + auto bounds_exchange = parthenon::AddBoundaryExchangeTasks( + rkl2_step_first | start_bnd, tl, base, pmesh->multilevel); + + tl.AddTask(bounds_exchange, parthenon::Update::FillDerived>, + base.get()); + } + + // Compute coefficients. Meyer+2012 eq. (16) + Real b_j = 1. / 3.; + Real b_jm1 = 1. / 3.; + Real b_jm2 = 1. / 3.; + Real w1 = 4. / (static_cast(s_rkl) * static_cast(s_rkl) + + static_cast(s_rkl) - 2.); + Real mu_j, nu_j, j, mu_tilde_j, gamma_tilde_j; + + // RKL loop + for (int jj = 2; jj <= s_rkl; jj++) { + j = static_cast(jj); + b_j = (j * j + j - 2.0) / (2 * j * (j + 1.0)); + mu_j = (2.0 * j - 1.0) / j * b_j / b_jm1; + nu_j = -(j - 1.0) / j * b_j / b_jm2; + mu_tilde_j = mu_j * w1; + gamma_tilde_j = -(1.0 - b_jm1) * mu_tilde_j; // -a_jm1*mu_tilde_j + + TaskRegion ®ion_calc_fluxes_step_other = ptask_coll->AddRegion(num_partitions); + for (int i = 0; i < num_partitions; i++) { + auto &tl = region_calc_fluxes_step_other[i]; + auto &base = pmesh->mesh_data.GetOrAdd("base", i); + + // Only need boundaries for base as it's the only "active" container exchanging + // data/fluxes with neighbors. All other containers are passive (i.e., data is only + // used but not exchanged). + const auto any = parthenon::BoundaryType::any; + auto start_bnd = tl.AddTask(none, parthenon::StartReceiveBoundBufs, base); + auto start_flxcor_recv = + tl.AddTask(none, parthenon::StartReceiveFluxCorrections, base); + + // Reset flux arrays (not guaranteed to be zero) + auto reset_fluxes = tl.AddTask(none, ResetFluxes, base.get()); + + // Calculate the diffusive fluxes for Yjm1 (here u1) + auto hydro_diff_fluxes = + tl.AddTask(reset_fluxes, CalcDiffFluxes, hydro_pkg.get(), base.get()); + + auto send_flx = + tl.AddTask(hydro_diff_fluxes, parthenon::LoadAndSendFluxCorrections, base); + auto recv_flx = + tl.AddTask(start_flxcor_recv, parthenon::ReceiveFluxCorrections, base); + auto set_flx = + tl.AddTask(recv_flx | hydro_diff_fluxes, parthenon::SetFluxCorrections, base); + + auto &Y0 = pmesh->mesh_data.GetOrAdd("u1", i); + auto &MY0 = pmesh->mesh_data.GetOrAdd("MY0", i); + auto &Yjm2 = pmesh->mesh_data.GetOrAdd("Yjm2", i); + + auto rkl2_step_other = + tl.AddTask(set_flx, RKL2StepOther, Y0.get(), base.get(), Yjm2.get(), MY0.get(), + mu_j, nu_j, mu_tilde_j, gamma_tilde_j, tau); + + // update ghost cells of base (currently storing Yj) + // Update ghost cells (local and non local), prolongate and apply bound cond. + // TODO(someone) experiment with split (local/nonlocal) comms with respect to + // performance for various tests (static, amr, block sizes) and then decide on the + // best impl. Go with default call (split local/nonlocal) for now. + // TODO(pgrete) optimize (in parthenon) to only send subset of updated vars + auto bounds_exchange = parthenon::AddBoundaryExchangeTasks( + rkl2_step_other | start_bnd, tl, base, pmesh->multilevel); + + tl.AddTask(bounds_exchange, parthenon::Update::FillDerived>, + base.get()); + } + + b_jm2 = b_jm1; + b_jm1 = b_j; + } +} + // See the advection.hpp declaration for a description of how this function gets called. TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { TaskCollection tc; @@ -232,6 +537,12 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { // First add split sources before the main time integration if (stage == 1) { + // If any tasks modify the conserved variables before this place, then + // the STS tasks should be updated to not assume prim and cons are in sync. + const auto &diffint = hydro_pkg->Param("diffint"); + if (diffint == DiffInt::rkl2) { + AddSTSTasks(&tc, pmesh, blocks, 0.5 * tm.dt); + } TaskRegion &strang_init_region = tc.AddRegion(num_partitions); for (int i = 0; i < num_partitions; i++) { auto &tl = strang_init_region[i]; @@ -277,7 +588,11 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { auto &tl = single_tasklist_per_pack_region[i]; auto &mu0 = pmesh->mesh_data.GetOrAdd("base", i); auto &mu1 = pmesh->mesh_data.GetOrAdd("u1", i); - tl.AddTask(none, parthenon::StartReceiveFluxCorrections, mu0); + + const auto any = parthenon::BoundaryType::any; + auto start_bnd = tl.AddTask(none, parthenon::StartReceiveBoundBufs, mu0); + auto start_flxcor_recv = + tl.AddTask(none, parthenon::StartReceiveFluxCorrections, mu0); const auto flux_str = (stage == 1) ? "flux_first_stage" : "flux_other_stage"; FluxFun_t *calc_flux_fun = hydro_pkg->Param(flux_str); @@ -298,9 +613,9 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { auto send_flx = tl.AddTask(first_order_flux_correct, parthenon::LoadAndSendFluxCorrections, mu0); - auto recv_flx = - tl.AddTask(first_order_flux_correct, parthenon::ReceiveFluxCorrections, mu0); - auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, mu0); + auto recv_flx = tl.AddTask(start_flxcor_recv, parthenon::ReceiveFluxCorrections, mu0); + auto set_flx = tl.AddTask(recv_flx | first_order_flux_correct, + parthenon::SetFluxCorrections, mu0); // compute the divergence of fluxes of conserved variables auto update = tl.AddTask( @@ -332,29 +647,14 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { tl.AddTask(source_split_strang_final, AddSplitSourcesFirstOrder, mu0.get(), tm); } - // Update ghost cells (local and non local) + // Update ghost cells (local and non local), prolongate and apply bound cond. // TODO(someone) experiment with split (local/nonlocal) comms with respect to // performance for various tests (static, amr, block sizes) and then decide on the // best impl. Go with default call (split local/nonlocal) for now. - parthenon::AddBoundaryExchangeTasks(source_split_first_order, tl, mu0, + parthenon::AddBoundaryExchangeTasks(source_split_first_order | start_bnd, tl, mu0, pmesh->multilevel); } - TaskRegion &async_region_3 = tc.AddRegion(num_task_lists_executed_independently); - for (int i = 0; i < blocks.size(); i++) { - auto &tl = async_region_3[i]; - auto &u0 = blocks[i]->meshblock_data.Get("base"); - auto prolongBound = none; - // Currently taken care of by AddBoundaryExchangeTasks above. - // Needs to be reintroduced once we reintroduce split (local/nonlocal) communication. - // if (pmesh->multilevel) { - // prolongBound = tl.AddTask(none, parthenon::ProlongateBoundaries, u0); - //} - - // set physical boundaries - auto set_bc = tl.AddTask(prolongBound, parthenon::ApplyBoundaryConditions, u0); - } - // Single task in single (serial) region to reset global vars used in reductions in the // first stage. if (stage == integrator->nstages && hydro_pkg->Param("calc_c_h")) { @@ -376,10 +676,21 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { auto &mu0 = pmesh->mesh_data.GetOrAdd("base", i); auto fill_derived = tl.AddTask(none, parthenon::Update::FillDerived>, mu0.get()); + } + const auto &diffint = hydro_pkg->Param("diffint"); + // If any tasks modify the conserved variables before this place and after FillDerived, + // then the STS tasks should be updated to not assume prim and cons are in sync. + if (diffint == DiffInt::rkl2 && stage == integrator->nstages) { + AddSTSTasks(&tc, pmesh, blocks, 0.5 * tm.dt); + } - if (stage == integrator->nstages) { - auto new_dt = tl.AddTask( - fill_derived, parthenon::Update::EstimateTimestep>, mu0.get()); + if (stage == integrator->nstages) { + TaskRegion &tr = tc.AddRegion(num_partitions); + for (int i = 0; i < num_partitions; i++) { + auto &tl = tr[i]; + auto &mu0 = pmesh->mesh_data.GetOrAdd("base", i); + auto new_dt = tl.AddTask(none, parthenon::Update::EstimateTimestep>, + mu0.get()); } } diff --git a/src/main.cpp b/src/main.cpp index 87d82021..93f984bb 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -44,6 +44,7 @@ int main(int argc, char *argv[]) { // Redefine defaults pman.app_input->ProcessPackages = Hydro::ProcessPackages; + pman.app_input->PreStepMeshUserWorkInLoop = Hydro::PreStepMeshUserWorkInLoop; const auto problem = pman.pinput->GetOrAddString("job", "problem_id", "unset"); if (problem == "linear_wave") { diff --git a/src/main.hpp b/src/main.hpp index bbcd0786..033118ad 100644 --- a/src/main.hpp +++ b/src/main.hpp @@ -35,7 +35,9 @@ enum class Reconstruction { undefined, dc, plm, ppm, wenoz, weno3, limo3 }; enum class Integrator { undefined, rk1, rk2, vl2, rk3 }; enum class Fluid { undefined, euler, glmmhd }; enum class Cooling { none, tabular }; -enum class Conduction { none, spitzer, thermal_diff }; +enum class Conduction { none, isotropic, anisotropic }; +enum class ConductionCoeff { none, fixed, spitzer }; +enum class DiffInt { none, unsplit, rkl2 }; enum class Hst { idx, ekin, emag, divb }; diff --git a/src/pgen/cloud.cpp b/src/pgen/cloud.cpp index b496dd2d..db53990a 100644 --- a/src/pgen/cloud.cpp +++ b/src/pgen/cloud.cpp @@ -32,6 +32,7 @@ using namespace parthenon::driver::prelude; Real rho_wind, mom_wind, rhoe_wind, r_cloud, rho_cloud; Real Bx = 0.0; Real By = 0.0; +Real Bz = 0.0; //======================================================================================== //! \fn void InitUserMeshData(Mesh *mesh, ParameterInput *pin) @@ -77,9 +78,13 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) { By = std::sqrt(2.0 * pressure / plasma_beta); } else if (mag_field_angle_str == "transverse") { Bx = std::sqrt(2.0 * pressure / plasma_beta); + } else if (mag_field_angle_str == "oblique") { + const auto B = std::sqrt(2.0 * pressure / plasma_beta); + Bx = B / std::sqrt(5.0); + Bz = 2 * Bx; } else { PARTHENON_FAIL("Unsupported problem/cloud/mag_field_angle. Please use either " - "'aligned' or 'transverse'."); + "'aligned', 'transverse', or 'oblique'."); } } @@ -152,7 +157,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const auto nscalars = hydro_pkg->Param("nscalars"); const bool mhd_enabled = hydro_pkg->Param("fluid") == Fluid::glmmhd; - if (((Bx != 0.0) || (By != 0.0)) && !mhd_enabled) { + if (((Bx != 0.0) || (By != 0.0) || (Bz != 0.0)) && !mhd_enabled) { PARTHENON_FAIL("Requested to initialize magnetic fields by `cloud/plasma_beta > 0`, " "but `hydro/fluid` is not supporting MHD."); } @@ -195,7 +200,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (mhd_enabled) { u(IB1, k, j, i) = Bx; u(IB2, k, j, i) = By; - u(IEN, k, j, i) += 0.5 * (Bx * Bx + By * By); + u(IB3, k, j, i) = Bz; + u(IEN, k, j, i) += 0.5 * (Bx * Bx + By * By + Bz * Bz); } // Init passive scalars @@ -222,6 +228,7 @@ void InflowWindX2(std::shared_ptr> &mbd, bool coarse) { const auto rhoe_wind_ = rhoe_wind; const auto Bx_ = Bx; const auto By_ = By; + const auto Bz_ = Bz; pmb->par_for_bndry( "InflowWindX2", nb, IndexDomain::inner_x2, parthenon::TopologicalElement::CC, coarse, KOKKOS_LAMBDA(const int, const int &k, const int &j, const int &i) { @@ -236,6 +243,10 @@ void InflowWindX2(std::shared_ptr> &mbd, bool coarse) { cons(IB2, k, j, i) = By_; cons(IEN, k, j, i) += 0.5 * By_ * By_; } + if (Bz_ != 0.0) { + cons(IB3, k, j, i) = Bz_; + cons(IEN, k, j, i) += 0.5 * Bz_ * Bz_; + } }); } diff --git a/src/pgen/diffusion.cpp b/src/pgen/diffusion.cpp index 85ac913a..24ae60f2 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, Athena Parthenon Collaboration. All rights reserved. +// Copyright (c) 2021-2023, Athena Parthenon Collaboration. All rights reserved. // Licensed under the 3-Clause License (the "LICENSE") // Parthenon headers @@ -27,7 +27,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const auto By = pin->GetOrAddReal("problem/diffusion", "By", 0.0); const auto iprob = pin->GetInteger("problem/diffusion", "iprob"); - const auto sigma = pin->GetOrAddReal("problem/diffusion", "sigma", 0.1); + 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"); + t0 = pin->GetOrAddReal("problem/diffusion", "t0", t0); + amp = pin->GetOrAddReal("problem/diffusion", "amp", amp); + } auto &coords = pmb->coords; @@ -64,7 +72,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } else if (iprob == 10) { u(IB1, k, j, i) = Bx; u(IB2, k, j, i) = By; - eint = 1 + std::exp(-SQR(coords.Xc<1>(i) / sigma) / 2.0); + // Adjust for anisotropic thermal conduction. + // If there's no conduction for the setup (because the field is perp.) + // treat as 1 (also in analysis) to prevent division by 0. + // Note, this is very constructed and needs to be updated/adjusted for isotropic + // conduction, other directions, and Bfield configs with |B| != 1 + Real eff_diff_coeff = Bx == 0.0 ? diff_coeff : diff_coeff * Bx * Bx; + eint = 1 + amp / std::sqrt(4. * M_PI * eff_diff_coeff * t0) * + std::exp(-(std::pow(coords.Xc<1>(i), 2.)) / + (4. * eff_diff_coeff * t0)); // Ring diffusion in x1-x2 plane } else if (iprob == 20) { const auto x = coords.Xc<1>(i); diff --git a/src/units.hpp b/src/units.hpp index a089e632..3290e87e 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -30,6 +30,7 @@ class Units { static constexpr parthenon::Real dyne_cm2_cgs = 1.0; // dyne/cm^2 static constexpr parthenon::Real msun_cgs = 1.98841586e+33; // g static constexpr parthenon::Real atomic_mass_unit_cgs = 1.660538921e-24; // g + static constexpr parthenon::Real electron_mass_cgs = 9.1093837015e-28; // g static constexpr parthenon::Real g_cm3_cgs = 1.0; // gcm**3 static constexpr parthenon::Real erg_cgs = 1; // erg static constexpr parthenon::Real gauss_cgs = 1; // gauss @@ -124,6 +125,7 @@ class Units { parthenon::Real atomic_mass_unit() const { return atomic_mass_unit_cgs / code_mass_cgs(); } + parthenon::Real electron_mass() const { return electron_mass_cgs / code_mass_cgs(); } parthenon::Real mh() const { return mh_cgs / code_mass_cgs(); } parthenon::Real erg() const { return erg_cgs / code_energy_cgs(); } parthenon::Real gauss() const { return gauss_cgs / code_magnetic_cgs(); } diff --git a/tst/regression/CMakeLists.txt b/tst/regression/CMakeLists.txt index 3461a817..c43b4de5 100644 --- a/tst/regression/CMakeLists.txt +++ b/tst/regression/CMakeLists.txt @@ -45,6 +45,9 @@ setup_test_both("aniso_therm_cond_ring_conv" "--driver ${PROJECT_BINARY_DIR}/bin setup_test_both("aniso_therm_cond_ring_multid" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \ --driver_input ${PROJECT_SOURCE_DIR}/inputs/diffusion.in --num_steps 4" "convergence") + +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("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_gauss_conv/__init__.py b/tst/regression/test_suites/aniso_therm_cond_gauss_conv/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tst/regression/test_suites/aniso_therm_cond_gauss_conv/aniso_therm_cond_gauss_conv.py b/tst/regression/test_suites/aniso_therm_cond_gauss_conv/aniso_therm_cond_gauss_conv.py new file mode 100644 index 00000000..9f002723 --- /dev/null +++ b/tst/regression/test_suites/aniso_therm_cond_gauss_conv/aniso_therm_cond_gauss_conv.py @@ -0,0 +1,221 @@ +# ======================================================================================== +# AthenaPK - a performance portable block structured AMR MHD code +# Copyright (c) 2020-2021, Athena Parthenon Collaboration. All rights reserved. +# Licensed under the 3-clause BSD License, see LICENSE file for details +# ======================================================================================== +# (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 for Los +# Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +# for the U.S. Department of Energy/National Nuclear Security Administration. All rights +# in the program are reserved by Triad National Security, LLC, and the U.S. Department +# of Energy/National Nuclear Security Administration. The Government is granted for +# itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +# license in this material to reproduce, prepare derivative works, distribute copies to +# the public, perform publicly and display publicly, and to permit others to do so. +# ======================================================================================== + +# 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 + +int_cfgs = ["unsplit", "rkl2"] +res_cfgs = [128, 256, 512] +field_cfgs = ["none", "aligned", "angle", "perp"] +tlim = 2.0 + +all_cfgs = list(itertools.product(res_cfgs, field_cfgs, int_cfgs)) + + +def get_outname(all_cfg): + res, field_cfg, int_cfg = all_cfg + return f"{res}_{field_cfg}_{int_cfg}" + + +def get_B(field_cfg): + if field_cfg == "aligned": + Bx = 1.0 + By = 0.0 + elif field_cfg == "perp": + Bx = 0.0 + By = 1.0 + elif field_cfg == "angle": + Bx = 1 / np.sqrt(2) + By = 1 / np.sqrt(2) + # isotropic case + elif field_cfg == "none": + Bx = 0.0 + By = 0.0 + else: + raise "Unknown field_cfg: %s" % field_cfg + + return Bx, By + + +class TestCase(utils.test_case.TestCaseAbs): + def Prepare(self, parameters, step): + assert parameters.num_ranks <= 4, "Use <= 4 ranks for diffusion test." + + res, field_cfg, int_cfg = all_cfgs[step - 1] + + Bx, By = get_B(field_cfg) + + outname = get_outname(all_cfgs[step - 1]) + + if field_cfg == "none": + conduction = "isotropic" + else: + conduction = "anisotropic" + + 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/nx2=32", + "parthenon/meshblock/nx2=32", + "parthenon/mesh/x2min=-1.0", + "parthenon/mesh/x2max=1.0", + "parthenon/mesh/nx3=1", + "parthenon/meshblock/nx3=1", + "parthenon/time/integrator=%s" + % ("rk2" if (int_cfg == "unsplit") else "rk1"), + "problem/diffusion/Bx=%f" % Bx, + "problem/diffusion/By=%f" % By, + "problem/diffusion/iprob=10", + "parthenon/output0/id=%s" % outname, + "hydro/gamma=2.0", + "parthenon/time/tlim=%f" % tlim, + "diffusion/conduction=%s" % conduction, + "diffusion/thermal_diff_coeff_code=0.25", + "diffusion/integrator=%s" % int_cfg, + ] + + 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, Bx, field_cfg): + eff_diff_coeff = 0.25 if Bx == 0.0 else 0.25 * Bx * Bx + tlim_ = 0.0 if field_cfg == "perp" else tlim + return 1.0 + 1e-6 / ( + np.sqrt(4 * np.pi * eff_diff_coeff * (0.5 + tlim_)) + / np.exp(-(x**2) / (4.0 * eff_diff_coeff * (0.5 + tlim_))) + ) + + 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(field_cfgs), len(int_cfgs), len(res_cfgs))) + for step in range(len(all_cfgs)): + outname = get_outname(all_cfgs[step]) + data_filename = f"{parameters.output_path}/parthenon.{outname}.final.phdf" + data_file = phdf.phdf(data_filename) + prim = data_file.Get("prim") + zz, yy, xx = data_file.GetVolumeLocations() + mask = yy == yy[0] + temp = prim[4][mask] + x = xx[mask] + res, field_cfg, int_cfg = all_cfgs[step] + row = res_cfgs.index(res) + col = int_cfgs.index(int_cfg) + + Bx, By = get_B(field_cfg) + temp_ref = get_ref(x, Bx, field_cfg) + l1 = np.average(np.abs(temp - temp_ref)) + l1_err[ + field_cfgs.index(field_cfg), + int_cfgs.index(int_cfg), + res_cfgs.index(res), + ] = l1 + p[row, col].plot(x, temp, label=f"{field_cfg} N={res} L$_1$={l1:.2g}") + + # Plot convergence + for i, field_cfg in enumerate(field_cfgs): + for j, int_cfg in enumerate(int_cfgs): + p[0, j].set_title(f"Integrator: {int_cfg}") + if field_cfg == "perp": + continue + + p[-1, j].plot( + res_cfgs, + l1_err[i, j, :], + label=f"{field_cfg} 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[i, 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.98: + print( + f"!!!\nConvergence for {field_cfg} test with {int_cfg} integrator " + f"is worse ({conv_measured}) than expected (-1.98).\n!!!" + ) + tests_passed = False + p[-1, j].plot( + res_cfgs, + np.exp(conv_a) * res_cfgs**conv_measured, + ":", + lw=0.75, + label=f"{field_cfg} 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 field_cfg in field_cfgs: + Bx, By = get_B(field_cfg) + for i in range(num_rows): + for j in range(num_cols): + y = get_ref(x, Bx, field_cfg) + 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, "cond.png"), + bbox_inches="tight", + dpi=300, + ) + + return tests_passed diff --git a/tst/regression/test_suites/aniso_therm_cond_ring_conv/README.md b/tst/regression/test_suites/aniso_therm_cond_ring_conv/README.md index bfd8226d..8222148a 100644 --- a/tst/regression/test_suites/aniso_therm_cond_ring_conv/README.md +++ b/tst/regression/test_suites/aniso_therm_cond_ring_conv/README.md @@ -2,5 +2,5 @@ Executes 2D ring diffusion problem following Sharma & Hammett (2007) and calculates convergence rate. Errors are calculated based on comparison to steady state solution. -Convergence for this problem is not great, but matches other numbers reported in literate, e.g., Balsara, Tilley & Howk MNRAS (2008) doi:10.1111/j.1365-2966.2008.13085.x . +Convergence for this problem is not great, but matches other numbers reported in literature, e.g., Balsara, Tilley & Howk MNRAS (2008) doi:10.1111/j.1365-2966.2008.13085.x . Also the minium temperature is checked to ensure that limiting is working (i.e., the temperature is nowhere below the initial background temperature).