diff --git a/inputs/diffusion/conduction.in b/inputs/diffusion/conduction.in index b59714cb..dc2a9555 100644 --- a/inputs/diffusion/conduction.in +++ b/inputs/diffusion/conduction.in @@ -61,6 +61,11 @@ gas = true gravity = true conduction = true drag = true +sts = false + + +integrator = rkl1 +sts_max_dt_ratio = 200.0 cfl = 0.3 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 95133ac3..85fb5f6a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -74,6 +74,9 @@ set (SRC_LIST radiation/imc/imc.hpp + sts/sts.cpp + sts/sts.hpp + utils/artemis_utils.cpp utils/artemis_utils.hpp utils/history.hpp diff --git a/src/artemis.cpp b/src/artemis.cpp index 782a6863..ffab8fd2 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -22,6 +22,7 @@ #include "gravity/gravity.hpp" #include "nbody/nbody.hpp" #include "rotating_frame/rotating_frame.hpp" +#include "sts/sts.hpp" #include "utils/artemis_utils.hpp" #include "utils/history.hpp" #include "utils/units.hpp" @@ -70,6 +71,9 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { const bool do_viscosity = pin->GetOrAddBoolean("physics", "viscosity", false); const bool do_conduction = pin->GetOrAddBoolean("physics", "conduction", false); const bool do_radiation = pin->GetOrAddBoolean("physics", "radiation", false); + const bool do_sts = pin->GetOrAddBoolean("physics", "sts", false); + + artemis->AddParam("do_sts", do_sts); artemis->AddParam("do_gas", do_gas); artemis->AddParam("do_dust", do_dust); artemis->AddParam("do_gravity", do_gravity); @@ -81,6 +85,11 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { artemis->AddParam("do_conduction", do_conduction); artemis->AddParam("do_diffusion", do_conduction || do_viscosity); artemis->AddParam("do_radiation", do_radiation); + + PARTHENON_REQUIRE(!(do_sts) || (do_sts && do_gas), + "STS requires the gas package, but there is not gas!"); + PARTHENON_REQUIRE(!(do_sts) || (do_sts && (do_conduction || do_viscosity)), + "STS requires diffusion to be enabled!"); PARTHENON_REQUIRE(!(do_cooling) || (do_cooling && do_gas), "Cooling requires the gas package, but there is not gas!"); PARTHENON_REQUIRE(!(do_viscosity) || (do_viscosity && do_gas), @@ -105,6 +114,8 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { if (do_rotating_frame) packages.Add(RotatingFrame::Initialize(pin.get())); if (do_cooling) packages.Add(Gas::Cooling::Initialize(pin.get())); if (do_drag) packages.Add(Drag::Initialize(pin.get())); + if (do_sts) packages.Add(STS::Initialize(pin.get())); + if (do_radiation) { auto eos_h = packages.Get("gas")->Param("eos_h"); auto opacity_h = packages.Get("gas")->Param("opacity_h"); diff --git a/src/artemis.hpp b/src/artemis.hpp index 5dc7d4a9..4bda42d4 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -104,6 +104,9 @@ enum class ArtemisBC { none }; +// ...STS integrator types +enum class STSInt { rkl1, rkl2, null }; + // Floating point limits template KOKKOS_FORCEINLINE_FUNCTION constexpr auto Big() { diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index e910b836..358bfb39 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -29,6 +29,7 @@ #include "nbody/nbody.hpp" #include "radiation/imc/imc.hpp" #include "rotating_frame/rotating_frame.hpp" +#include "sts/sts.hpp" #include "utils/integrators/artemis_integrator.hpp" using namespace parthenon::driver::prelude; @@ -65,6 +66,7 @@ ArtemisDriver::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in do_conduction = artemis_pkg->template Param("do_conduction"); do_nbody = artemis_pkg->template Param("do_nbody"); do_diffusion = do_viscosity || do_conduction; + do_sts = artemis_pkg->template Param("do_sts"); do_radiation = artemis_pkg->template Param("do_radiation"); // NBody initialization tasks @@ -103,10 +105,14 @@ TaskListStatus ArtemisDriver::Step() { // Prepare registers PreStepTasks(); + // Execute STS first stage + if (do_sts) STSFirstStage(); + // Execute explicit, unsplit physics auto status = StepTasks().Execute(); if (status != TaskListStatus::complete) return status; + // STS_second_stage(); // Execute operator split physics if (do_radiation) status = IMC::JaybenneIMC(pmesh, tm.time, tm.dt); if (status != TaskListStatus::complete) return status; @@ -138,6 +144,31 @@ void ArtemisDriver::PreStepTasks() { auto &u1 = pmesh->mesh_data.Add("u1", u0); } +//---------------------------------------------------------------------------------------- +//! \fn TaskCollection ArtemisDriver::STS_first_stage +//! \brief Define the tasks for the first stage of the STS integrator +template +void ArtemisDriver::STSFirstStage() { + + // Assign sts registers and First stage of STS integration + auto &gas_pkg = pmesh->packages.Get("gas"); + auto min_diff_dt = gas_pkg->template Param("diff_dt"); + auto &sts_pkg = pmesh->packages.Get("STS"); + auto info_output = sts_pkg->template Param("info_output"); + // compute the number of stages needed for the STS integrator (for rkl1 only) + int s_sts = static_cast(0.5 * (-1. + std::sqrt(1. + 8. * tm.dt / min_diff_dt))); + if (s_sts % 2 == 0) s_sts += 1; + + if (parthenon::Globals::my_rank == 0 && info_output) { + Real ratio = tm.dt / min_diff_dt; + std::cout << "STS ratio: " << ratio << ", Taking " << s_sts << " steps." << std::endl; + if (ratio > 200.0) { + std::cout << "WARNING: ratio is > 200. Proceed at own risk." << std::endl; + } + } + STS::PreStepSTSTasks(pmesh, tm.time, tm.dt, s_sts); +} + //---------------------------------------------------------------------------------------- //! \fn TaskCollection ArtemisDriver::PreStepTasks //! \brief Defines the main integrator's TaskCollection for the ArtemisDriver @@ -186,7 +217,7 @@ TaskCollection ArtemisDriver::StepTasks() { // Compute (gas) diffusive fluxes TaskID diff_flx = none; - if (do_diffusion && do_gas) { + if (do_diffusion && do_gas && !(do_sts)) { auto zf = tl.AddTask(none, Gas::ZeroDiffusionFlux, u0.get()); TaskID vflx = zf, tflx = zf; if (do_viscosity) vflx = tl.AddTask(zf, Gas::ViscousFlux, u0.get()); @@ -215,7 +246,8 @@ TaskCollection ArtemisDriver::StepTasks() { // NOTE(@pdmullen): I believe set_flx dependency implicitly inside gas_coord_src, // but included below explicitly for posterity TaskID gas_diff_src = gas_coord_src | diff_flx | set_flx; - if (do_diffusion && do_gas) { + + if (do_diffusion && do_gas && !(do_sts)) { gas_diff_src = tl.AddTask(gas_coord_src | diff_flx | set_flx, Gas::DiffusionUpdate, u0.get(), bdt); } @@ -252,18 +284,18 @@ TaskCollection ArtemisDriver::StepTasks() { tl.AddTask(cooling_src, ArtemisDerived::SetAuxillaryFields, u0.get()); // Set (remaining) fields to be communicated - auto c2p = tl.AddTask(set_aux, PreCommFillDerived>, u0.get()); + auto pre_comm = tl.AddTask(set_aux, PreCommFillDerived>, u0.get()); // Set boundary conditions (both physical and logical) - auto bcs = parthenon::AddBoundaryExchangeTasks(c2p, tl, u0, pmesh->multilevel); + auto bcs = parthenon::AddBoundaryExchangeTasks(pre_comm, tl, u0, pmesh->multilevel); - // Sync fields - auto p2c = tl.AddTask(TQ::local_sync, bcs, FillDerived>, u0.get()); + // Update primitive variables + auto c2p = tl.AddTask(TQ::local_sync, bcs, FillDerived>, u0.get()); // Advance nbody integrator - TaskID nbadv = p2c; + TaskID nbadv = c2p; if (do_nbody) { - nbadv = tl.AddTask(TQ::once_per_region, p2c, NBody::Advance, pmesh, time, stage, + nbadv = tl.AddTask(TQ::once_per_region, c2p, NBody::Advance, pmesh, time, stage, nbody_integrator.get()); } } @@ -281,6 +313,8 @@ TaskCollection ArtemisDriver::PostStepTasks() { TaskCollection tc; TaskID none(0); + // TODO: Implement RKL2 STS integration + const int num_partitions = pmesh->DefaultNumPartitions(); auto &post_region = tc.AddRegion(num_partitions); for (int i = 0; i < num_partitions; i++) { @@ -305,31 +339,37 @@ typedef ApplicationInput AI; template ArtemisDriver::ArtemisDriver(PI *p, AI *a, M *m, const bool r); template TaskListStatus ArtemisDriver::Step(); template void ArtemisDriver::PreStepTasks(); +template void ArtemisDriver::STSFirstStage(); template TaskCollection ArtemisDriver::StepTasks(); template TaskCollection ArtemisDriver::PostStepTasks(); template ArtemisDriver::ArtemisDriver(PI *p, AI *a, M *m, const bool r); template TaskListStatus ArtemisDriver::Step(); template void ArtemisDriver::PreStepTasks(); +template void ArtemisDriver::STSFirstStage(); template TaskCollection ArtemisDriver::StepTasks(); template TaskCollection ArtemisDriver::PostStepTasks(); template ArtemisDriver::ArtemisDriver(PI *p, AI *a, M *m, const bool r); template TaskListStatus ArtemisDriver::Step(); template void ArtemisDriver::PreStepTasks(); +template void ArtemisDriver::STSFirstStage(); template TaskCollection ArtemisDriver::StepTasks(); template TaskCollection ArtemisDriver::PostStepTasks(); template ArtemisDriver::ArtemisDriver(PI *p, AI *a, M *m, const bool r); template TaskListStatus ArtemisDriver::Step(); template void ArtemisDriver::PreStepTasks(); +template void ArtemisDriver::STSFirstStage(); template TaskCollection ArtemisDriver::StepTasks(); template TaskCollection ArtemisDriver::PostStepTasks(); template ArtemisDriver::ArtemisDriver(PI *p, AI *a, M *m, const bool r); template TaskListStatus ArtemisDriver::Step(); template void ArtemisDriver::PreStepTasks(); +template void ArtemisDriver::STSFirstStage(); template TaskCollection ArtemisDriver::StepTasks(); template TaskCollection ArtemisDriver::PostStepTasks(); template ArtemisDriver::ArtemisDriver(PI *p, AI *a, M *m, const bool r); template TaskListStatus ArtemisDriver::Step(); template void ArtemisDriver::PreStepTasks(); +template void ArtemisDriver::STSFirstStage(); template TaskCollection ArtemisDriver::StepTasks(); template TaskCollection ArtemisDriver::PostStepTasks(); diff --git a/src/artemis_driver.hpp b/src/artemis_driver.hpp index f3116a38..90ce9ce0 100644 --- a/src/artemis_driver.hpp +++ b/src/artemis_driver.hpp @@ -46,6 +46,7 @@ class ArtemisDriver : public EvolutionDriver { const bool is_restart_in); TaskListStatus Step(); void PreStepTasks(); + void STSFirstStage(); TaskCollection StepTasks(); TaskCollection PostStepTasks(); @@ -53,7 +54,7 @@ class ArtemisDriver : public EvolutionDriver { IntegratorPtr_t integrator, nbody_integrator; StateDescriptor *artemis_pkg; bool do_gas, do_dust, do_gravity, do_rotating_frame, do_cooling, do_drag, do_viscosity, - do_nbody, do_conduction, do_diffusion, do_radiation; + do_nbody, do_conduction, do_diffusion, do_sts, do_radiation; const bool is_restart; }; diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index ffbef0de..bbedfe94 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -18,6 +18,7 @@ #include "artemis.hpp" #include "gas.hpp" #include "geometry/geometry.hpp" +#include "sts/sts.hpp" #include "utils/artemis_utils.hpp" #include "utils/diffusion/diffusion.hpp" #include "utils/diffusion/diffusion_coeff.hpp" @@ -183,6 +184,15 @@ std::shared_ptr Initialize(ParameterInput *pin, const bool do_conduction = pin->GetOrAddBoolean("physics", "conduction", false); params.Add("do_conduction", do_conduction); + const bool do_sts = pin->GetOrAddBoolean("physics", "sts", false); + params.Add("do_sts", do_sts); + + if (do_sts) { + params.Add("diff_dt", std::numeric_limits::max(), Params::Mutability::Mutable); + Real sts_max_dt_ratio = pin->GetOrAddReal("sts", "sts_max_dt_ratio", -1.0); + params.Add("sts_max_dt_ratio", sts_max_dt_ratio); + } + const bool do_diffusion = do_viscosity || do_conduction; params.Add("do_diffusion", do_diffusion); @@ -388,6 +398,7 @@ std::shared_ptr Initialize(ParameterInput *pin, //---------------------------------------------------------------------------------------- //! \fn Real Gas::EstimateTimestepMesh //! \brief Compute gas hydrodynamics timestep +//! STS_Flag template Real EstimateTimestepMesh(MeshData *md) { using parthenon::MakePackDescriptor; @@ -464,7 +475,23 @@ Real EstimateTimestepMesh(MeshData *md) { Real diff_dt = std::min(visc_dt, cond_dt); const auto cfl_number = params.template Get("cfl"); - return cfl_number * std::min(min_dt, diff_dt); + + // STS Time Stepping Control + const auto do_sts = params.template Get("do_sts"); + if (do_sts) { + const auto sts_max_dt_ratio = params.template Get("sts_max_dt_ratio"); + auto dt_ratio = min_dt / diff_dt; + // limit the timestep within the STS ratio, otherwise use the hyperbolic timestep + if (sts_max_dt_ratio > 0.0 && dt_ratio > sts_max_dt_ratio) { + min_dt = sts_max_dt_ratio * diff_dt; + } + // update the parabolic timestep + gas_pkg->UpdateParam("diff_dt", cfl_number * diff_dt); + } else { + min_dt = std::min(min_dt, diff_dt); + } + + return cfl_number * min_dt; } //---------------------------------------------------------------------------------------- diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp new file mode 100644 index 00000000..4bc52d77 --- /dev/null +++ b/src/sts/sts.cpp @@ -0,0 +1,152 @@ +//======================================================================================== +// (C) (or copyright) 2023-2024. 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. +//======================================================================================== + +// C++ headers +#include +#include + +// Parthenon includes +#include + +// Artemis includes +#include "artemis.hpp" +#include "gas/gas.hpp" +#include "sts.hpp" +#include "utils/artemis_utils.hpp" + +using ArtemisUtils::VI; + +namespace STS { + +IntegratorPtr_t sts_integrator; + +//---------------------------------------------------------------------------------------- +//! \fn StateDescriptor STS ::Initialize +//! \brief Adds intialization function for STS package +std::shared_ptr Initialize(ParameterInput *pin) { + + auto STS = std::make_shared("STS"); + Params ¶ms = STS->AllParams(); + + // initial check of getting the physics needed for the sts integrator + // Determine input file specified physics + const bool do_gas = pin->GetOrAddBoolean("physics", "gas", true); + const bool do_viscosity = pin->GetOrAddBoolean("physics", "viscosity", false); + const bool do_conduction = pin->GetOrAddBoolean("physics", "conduction", false); + const bool do_sts = pin->GetOrAddBoolean("physics", "sts", false); + const bool do_diffusion = do_conduction || do_viscosity; + + params.Add("do_sts", do_sts); + params.Add("do_gas", do_gas); + params.Add("do_viscosity", do_viscosity); + params.Add("do_conduction", do_conduction); + params.Add("do_diffusion", do_diffusion); + + if (!do_diffusion) { + PARTHENON_FAIL("STS integrator requires diffusion to be enabled!"); + } + + // Getting the integrator & time ratio between hyperbolic and parabolic terms + std::string sts_intg_mothod = pin->GetOrAddString("sts", "integrator", "none"); + Real sts_max_dt_ratio = pin->GetOrAddReal("sts", "sts_max_dt_ratio", -1.0); + const bool info_output = pin->GetOrAddBoolean("sts", "info_output", false); + + STSInt sts_intg_mothod_param = STSInt::null; + if (sts_intg_mothod == "rkl1") { + sts_intg_mothod_param = STSInt::rkl1; + sts_integrator = std::make_unique("rk1"); + } else if (sts_intg_mothod == "rkl2") { + PARTHENON_FAIL("rkl2 STS integrator not implemented!"); + sts_intg_mothod_param = STSInt::rkl2; + } else { + PARTHENON_FAIL("STS integrator not recognized!"); + } + + params.Add("sts_intg_mothod", sts_intg_mothod_param); + params.Add("sts_max_dt_ratio", sts_max_dt_ratio); + params.Add("info_output", info_output); + + return STS; +} + +//---------------------------------------------------------------------------------------- +//! \fn PreStepSTSTasks +//! \brief Executes the pre-step tasks for the STS integrator +template +void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int nstages) { + + // Getting the integrator & time ratio between hyperbolic and parabolic terms + auto &sts = pmesh->packages.Get("STS"); + const STSInt sts_intg_mothod = sts->Param("sts_intg_mothod"); + + // Check if the integrator is set + if (sts_intg_mothod == STSInt::null) { + PARTHENON_FAIL("STS integrator not set!"); + } + + // Execute the integrator tasks + if (sts_intg_mothod == STSInt::rkl1) { + // RKL1 : Full timestep dt_sts + for (int stage = 1; stage <= nstages; ++stage) { + //----------------------------------------------------------------- + // RKL1 STS update + // Y_{j} = nuj*Y_{j-2} + muj*Y_{j-1} + dt_sts*muj_tilde*F_diff(Y_{j-1}), + // where F_diff(Y_{j-1}) = divf/vol + // + // Set up the STS stage coefficients + // gam1 = nuj = (2.*j - 1.)/j; + // gam0 = muj = (1. - j)/j; + // beta_dt/dt = muj_tilde = pm->muj*2./(std::pow(s, 2.) + s); + // + // Update strategy for the registers + // Let u0 = Y_{j-1}, u1 = Y_{j-2} + // 1. After ApplyUpdate: u1 = nuj*Y_{j-2} + muj*Y_{j-1}, u1 -> Y'_{j} + // 2. swap u0 <-> u1, u0 -> Y'_{j}, u0 -> Y_{j-1} + // 3. After DiffusionUpdate: u0 = Y_{j}, u1 = Y_{j-1} + + Real muj = (2. * stage - 1.) / stage; + Real nuj = (1. - stage) / stage; + Real muj_tilde = muj * 2. / (std::pow(nstages, 2.) + nstages); + Real bdt = muj_tilde * dt; + + // We always use the stage 1 state for the coefficients + sts_integrator->beta[0] = 0.0; + sts_integrator->gam0[0] = nuj; // since we swap u0 and u1 + sts_integrator->gam1[0] = muj; + STSRKL1(pmesh, time, bdt, stage, nstages).Execute(); + } + } else if (sts_intg_mothod == STSInt::rkl2) { + PARTHENON_FAIL("STS rkl2 integrator not implemented!"); + // (TODO) RKL2 : // eq (21) using half hyperbolic timestep + // due to Strang split + // STSRKL2FirstStage(pmesh, time, 0.5*dt, nstages); + } +} + +//---------------------------------------------------------------------------------------- +//! template instantiations +typedef Coordinates C; +typedef Mesh M; +// PreStepSTSTasks template instantiations +template void PreStepSTSTasks(M *m, const Real time, Real dt, int nstages); +template void PreStepSTSTasks(M *m, const Real time, Real dt, + int nstages); +template void PreStepSTSTasks(M *m, const Real time, Real dt, + int nstages); +template void PreStepSTSTasks(M *m, const Real time, Real dt, + int nstages); +template void PreStepSTSTasks(M *m, const Real time, Real dt, + int nstages); +template void PreStepSTSTasks(M *m, const Real time, Real dt, + int nstages); +} // namespace STS \ No newline at end of file diff --git a/src/sts/sts.hpp b/src/sts/sts.hpp new file mode 100644 index 00000000..6eccbe9e --- /dev/null +++ b/src/sts/sts.hpp @@ -0,0 +1,131 @@ +//======================================================================================== +// (C) (or copyright) 2023-2024. 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. +//======================================================================================== +#ifndef STS_STS_HPP_ +#define STS_STS_HPP_ + +#include "artemis.hpp" +#include "derived/fill_derived.hpp" +#include "gas/gas.hpp" +#include "geometry/geometry.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/integrators/artemis_integrator.hpp" +#include "utils/units.hpp" + +namespace STS { + +using Integrator_t = parthenon::LowStorageIntegrator; +using IntegratorPtr_t = std::unique_ptr; + +// Global variable declaration +extern IntegratorPtr_t sts_integrator; + +std::shared_ptr Initialize(ParameterInput *pin); + +//---------------------------------------------------------------------------------------- +//! \fn STSRKL1 +//! \brief Assembles the tasks for the STS RKL1 integrator +// comment: Maybe it should be moved back to artemis_driver.cpp +template +TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nstages) { + + using TQ = TaskQualifier; + TaskCollection tc; + + using namespace ::parthenon::Update; + TaskID none(0); + const auto any = parthenon::BoundaryType::any; + const int num_partitions = pmesh->DefaultNumPartitions(); + + auto &pkg = pmesh->packages.Get("artemis"); + const auto do_viscosity = pkg->template Param("do_viscosity"); + const auto do_conduction = pkg->template Param("do_conduction"); + const auto do_diffusion = pkg->template Param("do_diffusion"); + const auto do_gas = pkg->template Param("do_gas"); + + // Deep copy u0 into u1 for integrator logic + if (stage == 1) { + auto &init_region = tc.AddRegion(num_partitions); + for (int i = 0; i < num_partitions; i++) { + auto &tl = init_region[i]; + auto &u0 = pmesh->mesh_data.GetOrAdd("u0", i); + auto &u1 = pmesh->mesh_data.GetOrAdd("u1", i); + tl.AddTask(none, ArtemisUtils::DeepCopyConservedData, u1.get(), u0.get()); + } + } + + TaskRegion &tr = tc.AddRegion(num_partitions); + for (int i = 0; i < num_partitions; i++) { + auto &tl = tr[i]; + auto &u0 = pmesh->mesh_data.GetOrAdd("u0", i); + auto &u1 = pmesh->mesh_data.GetOrAdd("u1", i); + + // Start looking for incoming messages (including for flux correction) + auto start_recv_u0 = tl.AddTask(none, parthenon::StartReceiveBoundBufs, u0); + auto start_flx_recv_u0 = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, u0); + + // Compute (gas) diffusive fluxes + TaskID diff_flx = none; + auto zf = tl.AddTask(none, Gas::ZeroDiffusionFlux, u0.get()); + TaskID vflx = zf, tflx = zf; + if (do_viscosity) vflx = tl.AddTask(zf, Gas::ViscousFlux, u0.get()); + if (do_conduction) tflx = tl.AddTask(zf | vflx, Gas::ThermalFlux, u0.get()); + diff_flx = vflx | tflx; + + // TODO(KWHO) Dust diffusion fluxes in the future + + // Communicate and set fluxes + auto send_flx = tl.AddTask( + diff_flx, parthenon::SendBoundBufs, u0); + auto recv_flx_u0 = + tl.AddTask(start_flx_recv_u0, parthenon::ReceiveFluxCorrections, u0); + auto set_flx_u0 = tl.AddTask(recv_flx_u0, parthenon::SetFluxCorrections, u0); + + // Apply flux divergence + auto update = none; + update = tl.AddTask(diff_flx | set_flx_u0, ArtemisUtils::ApplyUpdate, u1.get(), + u0.get(), 1, sts_integrator.get()); + + // swap u0 <-> u1 + auto swap_data_1 = tl.AddTask(update, ArtemisUtils::SwapData, u0.get(), u1.get()); + + TaskID gas_diff_src = swap_data_1 | diff_flx | set_flx_u0; + gas_diff_src = tl.AddTask(swap_data_1 | diff_flx | set_flx_u0, + Gas::DiffusionUpdate, u0.get(), dt); + + // Set auxillary fields + auto set_aux_u0 = + tl.AddTask(gas_diff_src, ArtemisDerived::SetAuxillaryFields, u0.get()); + + // Set (remaining) fields to be communicated + auto pre_comm_u0 = tl.AddTask(set_aux_u0, // update, + PreCommFillDerived>, u0.get()); + + // Set boundary conditions (both physical and logical) + auto bcs_u0 = + parthenon::AddBoundaryExchangeTasks(pre_comm_u0, tl, u0, pmesh->multilevel); + + // Update primitive variables + auto c2p_u0 = + tl.AddTask(TQ::local_sync, bcs_u0, FillDerived>, u0.get()); + } + + return tc; +} + +// STS integrator functions +template +void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int nstages); + +} // namespace STS + +#endif // STS_STS_HPP_ \ No newline at end of file diff --git a/src/utils/integrators/artemis_integrator.hpp b/src/utils/integrators/artemis_integrator.hpp index 537615d1..915823e6 100644 --- a/src/utils/integrators/artemis_integrator.hpp +++ b/src/utils/integrators/artemis_integrator.hpp @@ -50,6 +50,34 @@ inline TaskStatus DeepCopyConservedData(MeshData *to, MeshData *from return TaskStatus::complete; } +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus ArtemisUtils::SwapData +//! \brief swap data between two MeshData objects +inline TaskStatus SwapData(MeshData *u0, MeshData *u1) { + using parthenon::MakePackDescriptor; + using parthenon::variable_names::any; + + std::vector flags({Metadata::Conserved}); + static auto desc = MakePackDescriptor(u0, flags); + const auto vt = desc.GetPack(u0); + const auto vf = desc.GetPack(u1); + const auto ibe = u0->GetBoundsI(IndexDomain::entire); + const auto jbe = u0->GetBoundsJ(IndexDomain::entire); + const auto kbe = u0->GetBoundsK(IndexDomain::entire); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "SwapData", parthenon::DevExecSpace(), 0, u0->NumBlocks() - 1, + kbe.s, kbe.e, jbe.s, jbe.e, ibe.s, ibe.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + for (int n = vt.GetLowerBound(b); n <= vt.GetUpperBound(b); ++n) { + auto swap_data = vt(b, n, k, j, i); + vt(b, n, k, j, i) = vf(b, n, k, j, i); + vf(b, n, k, j, i) = swap_data; + } + }); + return TaskStatus::complete; +} + //---------------------------------------------------------------------------------------- //! \fn TaskStatus ArtemisUtils::ApplyUpdate //! \brief @@ -75,7 +103,6 @@ TaskStatus ApplyUpdate(MeshData *u0, MeshData *u1, const int stage, const auto kb = u0->GetBoundsK(IndexDomain::interior); const bool multi_d = (pm->ndim > 1); const bool three_d = (pm->ndim > 2); - parthenon::par_for( DEFAULT_LOOP_PATTERN, "ApplyUpdate", parthenon::DevExecSpace(), 0, u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, diff --git a/tst/scripts/diffusion/thermal_diffusion.py b/tst/scripts/diffusion/thermal_diffusion.py index 15c4ac60..e721da82 100644 --- a/tst/scripts/diffusion/thermal_diffusion.py +++ b/tst/scripts/diffusion/thermal_diffusion.py @@ -51,6 +51,7 @@ def run(**kwargs): "parthenon/job/problem_id=" + name, "artemis/coordinates=" + g, "parthenon/time/tlim=50.0", + "physics/sts=true", "gas/conductivity/cond={:.8f}".format(_kcond), "gravity/uniform/gx1={:.8f}".format(_gx1), "problem/flux={:.8f}".format(_flux),