From 06ddc808d5cb9f4935802e3e64f71bbda2cd5b36 Mon Sep 17 00:00:00 2001 From: Ka Wai Ho Date: Mon, 6 Jan 2025 14:52:10 -0700 Subject: [PATCH 01/11] Layout the Infrastructure for STS algorthim 1. Added sts folder and files to handle the update. 2. artemis.cpp/hpp modified a) added sts flag b) added Initialize for sts 3. artemis_driver.hpp a) added sts flag b) added the entry for sts in PreStepTasks() fn c) added the entry for sts in PostStepTasks() fn d) Added a check to remove diffusion from main step if sts is enabled 4. gas.cpp a) Added the sts flag b) Added the diffusion timetep for paramters c) Modified the timestep compuation --- src/artemis.cpp | 6 +- src/artemis.hpp | 3 + src/artemis_driver.cpp | 47 ++++++++++--- src/artemis_driver.hpp | 2 +- src/gas/gas.cpp | 27 +++++++- src/sts/sts.cpp | 154 +++++++++++++++++++++++++++++++++++++++++ src/sts/sts.hpp | 23 ++++++ 7 files changed, 251 insertions(+), 11 deletions(-) create mode 100644 src/sts/sts.cpp create mode 100644 src/sts/sts.hpp diff --git a/src/artemis.cpp b/src/artemis.cpp index 79da0a81..26eaadac 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -70,6 +70,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); @@ -86,7 +89,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { PARTHENON_REQUIRE(!(do_viscosity) || (do_viscosity && do_gas), "Viscosity requires the gas package, but there is not gas!"); PARTHENON_REQUIRE(!(do_conduction) || (do_conduction && do_gas), - "Conduction requires the gas package, but there is not gas!"); + "Conduction requires the gas package, but there is not gas!"); PARTHENON_REQUIRE(!(do_radiation) || (do_radiation && do_gas), "Radiation requires the gas package, but there is not gas!"); @@ -105,6 +108,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { if (do_cooling) packages.Add(Gas::Cooling::Initialize(pin.get())); if (do_drag) packages.Add(Drag::Initialize(pin.get())); if (do_nbody) packages.Add(NBody::Initialize(pin.get(), constants)); + 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..a6602c41 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -30,6 +30,7 @@ #include "radiation/imc/imc.hpp" #include "rotating_frame/rotating_frame.hpp" #include "utils/integrators/artemis_integrator.hpp" +#include "sts/sts.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 @@ -136,6 +138,33 @@ void ArtemisDriver::PreStepTasks() { auto &base = pmesh->mesh_data.Get(); auto &u0 = pmesh->mesh_data.AddShallow("u0", base, names); auto &u1 = pmesh->mesh_data.Add("u1", u0); + + // 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"); + if (do_sts) { + // compute the number of stages needed for the STS integrator + int s_sts = + static_cast(0.5 * (std::sqrt(9.0 + 16.0 * tau / min_diff_dt) - 1.0)) + 1; + if (s_sts % 2 == 0) s_sts += 1; + + if (parthenon::Globals::my_rank == 0) { + const auto ratio = 2.0 * tau / mindt_diff; + std::cout << "STS ratio: " << ratio << ", Taking " << s_sts << " steps." << std::endl; + if (ratio > 400.1) { + std::cout << "WARNING: ratio is > 400. Proceed at own risk." << std::endl; + } + } + + if (STSInt::rkl1){ + // (TODO) RKL1 : Full timestep dt_sts + //STSRKL1(pmesh, tm.time, tm.dt, s_sts); + }else if (STSInt::rkl2){ + // (TODO) RKL2 : // eq (21) using half hyperbolic timestep + // due to Strang split + //STSRKL2(pmesh, tm.time, 0.5*tm.dt, s_sts); + } + } } //---------------------------------------------------------------------------------------- @@ -186,7 +215,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 +244,7 @@ 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 +281,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 +310,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++) { diff --git a/src/artemis_driver.hpp b/src/artemis_driver.hpp index f3116a38..20e839c2 100644 --- a/src/artemis_driver.hpp +++ b/src/artemis_driver.hpp @@ -53,7 +53,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 92d1c78e..cb1b8b03 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" @@ -182,6 +183,14 @@ 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) { + Real diff_dt = Big(); + params.Add("diff_dt", diff_dt); + } + const bool do_diffusion = do_viscosity || do_conduction; params.Add("do_diffusion", do_diffusion); @@ -386,6 +395,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; @@ -462,7 +472,22 @@ 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"); + // limit the timestep within the STS ratio, otherwise use the hyperbolic timestep + if (sts_max_dt_ratio > 0.0 && dt_hyp > sts_max_dt_ratio*diff_dt) { + min_dt = std::min(min_dt, max_dt_ratio * diff_dt); + } + // update the parabolic timestep + gas_pkg->UpdateParam("diff_dt", cfl*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..dee8a676 --- /dev/null +++ b/src/sts/sts.cpp @@ -0,0 +1,154 @@ +//======================================================================================== +// (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 "sts.hpp" +#include "utils/artemis_utils.hpp" + + +namespace STS{ +//---------------------------------------------------------------------------------------- +//! \fn StateDescriptor STS ::Initialize +//! \brief Adds intialization function for STS package +std::shared_ptr Initialize(ParameterInput *pin) { + + // Extract artemis package + artemis_pkg = pm->packages.Get("artemis").get(); + + // initial check of getting the physics needed for the sts integrator + do_viscosity = artemis_pkg->template Param("do_viscosity"); + do_conduction = artemis_pkg->template Param("do_conduction"); + bool do_diffusion = do_viscosity || do_conduction; + + 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_integrator = pin->GetOrAddString("sts", "integrator", "none"); + Real sts_max_dt_ratio = pin->GetOrAddReal("sts","sts_max_dt_ratio", -1.0); + + STSInt Diff_integrator = STSInt::null; + if (sts_integrator == "rkl1") { + Diff_integrator = STSInt::rkl1; + } else if (sts_integrator == "rkl2") { + PARTHENON_FAIL("rkl2 STS integrator not implemented!"); + Diff_integrator = STSInt::rkl2; + } else { + PARTHENON_FAIL("STS integrator not recognized!"); + } + + artemis_pkg->AddParam("sts_integrator", sts_integrator); + artemis_pkg->AddParam("sts_max_dt_ratio", sts_max_dt_ratio); + +} + +//---------------------------------------------------------------------------------------- +//! \fn STSRKL1 +//! \brief Assembles the tasks for the STS RKL1 integrator +// comment: Maybe it should be moved back to artemis_driver.cpp +void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages) { + using namespace ::parthenon::Update; + const auto any = parthenon::BoundaryType::any; + const int num_partitions = pmesh->DefaultNumPartitions(); + + // Deep copy u0 into u1 for integrator logic + 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()); + } + + //---------------------------------------------------------------------------------------- + // gam0 = nuj + // gam1 = muj + // beta_dt = dt_sts*muj_tilde ( but be careful with the puting v1 in divf calculation) + for (int stage = 1; stage <= nstages; stage++) { + // Set up the STS stage coefficients + // v0 = Y_{n-2} + // v1 = Y_{n-1} + // gam1 = muj = (2.*j - 1.)/j; + // gam0 = nuj = (1. - j)/j; + // beta_dt/dt = muj_tilde = pm->muj*2./(std::pow(s, 2.) + s); + Real muj = (1. - stage)/stage; + Real nuj = (2.*stage - 1.)/stage; + Real muj_tilde = (2.*stage - 1.)/stage * 2./(std::pow(nstages, 2.) + nstages); + + const Real bdt = dt; + + 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 = tl.AddTask(none, parthenon::StartReceiveBoundBufs, u0); + auto start_flx_recv = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, u0); + + // Compute (gas) diffusive fluxes + TaskID diff_flx = none; + if (do_diffusion && do_gas) { + 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; + } + + // Communicate and set fluxes + auto send_flx = + tl.AddTask(gas_flx | dust_flx | diff_flx, + parthenon::SendBoundBufs, u0); + auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0); + auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0); + + // Apply flux divergence, STS need stage to 0 for the sts ceofficients + auto update = + tl.AddTask(gas_flx | dust_flx | set_flx, ArtemisUtils::ApplyUpdate, + u0.get(), u1.get(), 0, integrator.get()); + } + + } + +} + +//---------------------------------------------------------------------------------------- +//! \fn STSRKL2FirstStage +//! \brief Assembles the tasks for first stage of the STS RKL2 integrator +void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages) { + // TODO: Implement RKL2 STS integration +} + + + +//---------------------------------------------------------------------------------------- +//! \fn STSRKL2SecondStage +//! \brief Assembles the tasks for first stage of the STS RKL2 integrator +void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages) { + // TODO: Implement RKL2 STS integration +} + + + +} \ No newline at end of file diff --git a/src/sts/sts.hpp b/src/sts/sts.hpp new file mode 100644 index 00000000..c19bb7e2 --- /dev/null +++ b/src/sts/sts.hpp @@ -0,0 +1,23 @@ +//======================================================================================== +// (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 "utils/units.hpp" + +namespace STS { + void STSRKL1( Mesh *pm, const Real time, Real dt, int nstages); + void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages); + void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages); +} \ No newline at end of file From 9dc6ac65dc45a69dad15ec86d315b95fa09ac4ad Mon Sep 17 00:00:00 2001 From: Ka Wai Ho Date: Tue, 7 Jan 2025 16:55:33 -0700 Subject: [PATCH 02/11] STS RKL1 Integrator implmentation - changed `sts.cpp/hpp` for implmenting the RKL1 - Bug Fixed in enterting sts_integrator section in `artemis_driver.cpp` --- src/artemis_driver.cpp | 14 ++++---- src/sts/sts.cpp | 77 ++++++++++++++++++++++++++++++++++++------ src/sts/sts.hpp | 8 ++++- 3 files changed, 80 insertions(+), 19 deletions(-) diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index a6602c41..a94a1b4f 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -145,21 +145,21 @@ void ArtemisDriver::PreStepTasks() { if (do_sts) { // compute the number of stages needed for the STS integrator int s_sts = - static_cast(0.5 * (std::sqrt(9.0 + 16.0 * tau / min_diff_dt) - 1.0)) + 1; + static_cast(0.5 * (std::sqrt(9.0 + 16.0 * tm.dt / min_diff_dt) - 1.0)) + 1; if (s_sts % 2 == 0) s_sts += 1; if (parthenon::Globals::my_rank == 0) { - const auto ratio = 2.0 * tau / mindt_diff; + const auto ratio = 2.0 * tm.dt / min_diff_dt; std::cout << "STS ratio: " << ratio << ", Taking " << s_sts << " steps." << std::endl; if (ratio > 400.1) { std::cout << "WARNING: ratio is > 400. Proceed at own risk." << std::endl; } } - - if (STSInt::rkl1){ + auto sts_integrator = artemis_pkg->template Param("sts_integrator"); + if (sts_integrator == STSInt::rkl1){ // (TODO) RKL1 : Full timestep dt_sts //STSRKL1(pmesh, tm.time, tm.dt, s_sts); - }else if (STSInt::rkl2){ + }else if (sts_integrator == STSInt::rkl2){ // (TODO) RKL2 : // eq (21) using half hyperbolic timestep // due to Strang split //STSRKL2(pmesh, tm.time, 0.5*tm.dt, s_sts); @@ -215,7 +215,7 @@ TaskCollection ArtemisDriver::StepTasks() { // Compute (gas) diffusive fluxes TaskID diff_flx = none; - if (do_diffusion && do_gas && !do_sts) { + 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()); @@ -244,7 +244,7 @@ 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 && !do_sts) { + if (do_diffusion && do_gas && !(do_sts) { gas_diff_src = tl.AddTask(gas_coord_src | diff_flx | set_flx, Gas::DiffusionUpdate, u0.get(), bdt); } diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp index dee8a676..19f07d12 100644 --- a/src/sts/sts.cpp +++ b/src/sts/sts.cpp @@ -79,10 +79,6 @@ void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages) { tl.AddTask(none, ArtemisUtils::DeepCopyConservedData, u1.get(), u0.get()); } - //---------------------------------------------------------------------------------------- - // gam0 = nuj - // gam1 = muj - // beta_dt = dt_sts*muj_tilde ( but be careful with the puting v1 in divf calculation) for (int stage = 1; stage <= nstages; stage++) { // Set up the STS stage coefficients // v0 = Y_{n-2} @@ -94,8 +90,6 @@ void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages) { Real nuj = (2.*stage - 1.)/stage; Real muj_tilde = (2.*stage - 1.)/stage * 2./(std::pow(nstages, 2.) + nstages); - const Real bdt = dt; - TaskRegion &tr = tc.AddRegion(num_partitions); for (int i = 0; i < num_partitions; i++) { auto &tl = tr[i]; @@ -118,19 +112,18 @@ void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages) { // Communicate and set fluxes auto send_flx = - tl.AddTask(gas_flx | dust_flx | diff_flx, + tl.AddTask(diff_flx, parthenon::SendBoundBufs, u0); auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0); auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0); // Apply flux divergence, STS need stage to 0 for the sts ceofficients auto update = - tl.AddTask(gas_flx | dust_flx | set_flx, ArtemisUtils::ApplyUpdate, - u0.get(), u1.get(), 0, integrator.get()); + tl.AddTask(set_flx, RKL1FluxUpadte, + u0.get(), u1.get(), + dt, muj, nuj, muj_tilde); } - } - } //---------------------------------------------------------------------------------------- @@ -150,5 +143,67 @@ void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages) { } +//---------------------------------------------------------------------------------------- +//! \fn RKL1FluxUpadte +//! \brief Applies the STS RKL1 update to the conserved variables +template +TaskStatus RKL1FluxUpadte(MeshData *u0, MeshData *u1, const Real dt, + const Real muj, const Real nuj, const Real muj_tilde) { + using parthenon::MakePackDescriptor; + using parthenon::variable_names::any; + auto pm = u0->GetParentPointer(); + + // Packing and indexing + std::vector flags({Metadata::Conserved, Metadata::WithFluxes}); + static auto desc = MakePackDescriptor(u0, flags, {parthenon::PDOpt::WithFluxes}); + const auto v0 = desc.GetPack(u0); + const auto v1 = desc.GetPack(u1); + const auto ib = u0->GetBoundsI(IndexDomain::interior); + const auto jb = u0->GetBoundsJ(IndexDomain::interior); + 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, "RKL1FluxUpadte", parthenon::DevExecSpace(), 0, + u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + // Extract coordinates + using parthenon::TopologicalElement; + geometry::Coords coords(v0.GetCoordinates(b), k, j, i); + + const auto ax1 = coords.GetFaceAreaX1(); + const auto ax2 = (multi_d) ? coords.GetFaceAreaX2() : NewArray(0.0); + const auto ax3 = (three_d) ? coords.GetFaceAreaX3() : NewArray(0.0); + + const Real vol = coords.Volume(); + + for (int n = v0.GetLowerBound(b); n <= v0.GetUpperBound(b); ++n) { + // compute flux divergence + Real divf = (ax1[0] * v0.flux(b, X1DIR, n, k, j, i) - + ax1[1] * v0.flux(b, X1DIR, n, k, j, i + 1)); + if (multi_d) + divf += (ax2[0] * v0.flux(b, X2DIR, n, k, j, i) - + ax2[1] * v0.flux(b, X2DIR, n, k, j + 1, i)); + if (three_d) + divf += (ax3[0] * v0.flux(b, X3DIR, n, k, j, i) - + ax3[1] * v0.flux(b, X3DIR, n, k + 1, j, i)); + + //---------------------------------------------------------------------------------------- + // Apply STS RKL1 update + // Y_{m} = nuj*Y_{j-2} + muj*Y_{j-1} + dt_sts*muj_tilde*M(Y_{j-1}) + // = nuj*Y_{j-2} + muj*Y_{j-1} + dt_sts*muj_tilde*(divf/vol) + // v0 = Y_{j-1}, v1 = Y_{j-2} + Real Y_jm2 = v0(b, n, k, j, i); + v0(b, n, k, j, i) = + nuj * v1(b, n, k, j, i) + muj * v0(b, n, k, j, i) + divf * dt * muj_tilde/ vol; + + // ---------------------------------------------------------------------------------------- + // Rearrange the variables for the next step + v1(b, n, k, j, i) = Y_jm2; + } + }); + return TaskStatus::complete; +} } \ No newline at end of file diff --git a/src/sts/sts.hpp b/src/sts/sts.hpp index c19bb7e2..181fa52a 100644 --- a/src/sts/sts.hpp +++ b/src/sts/sts.hpp @@ -17,7 +17,13 @@ #include "utils/units.hpp" namespace STS { - void STSRKL1( Mesh *pm, const Real time, Real dt, int nstages); + std::shared_ptr Initialize(ParameterInput *pin); + // task list functions + void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages); void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages); void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages); + + // task status functions + TaskStatus RKL1FluxUpadte(MeshData *u0, MeshData *u1, const Real dt, + const Real muj, const Real nuj, const Real muj_tilde); } \ No newline at end of file From 99a7de7b111e92ed7570e6b8c747b063bff09e53 Mon Sep 17 00:00:00 2001 From: Ka Wai HO Date: Thu, 9 Jan 2025 01:31:47 -0700 Subject: [PATCH 03/11] minor bug fix & added TODO --- src/artemis_driver.cpp | 2 +- src/sts/sts.cpp | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index a94a1b4f..0e763e68 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -244,7 +244,7 @@ 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 && !(do_sts) { + if (do_diffusion && do_gas && !(do_sts)) { gas_diff_src = tl.AddTask(gas_coord_src | diff_flx | set_flx, Gas::DiffusionUpdate, u0.get(), bdt); } diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp index 19f07d12..bbf3994a 100644 --- a/src/sts/sts.cpp +++ b/src/sts/sts.cpp @@ -81,8 +81,8 @@ void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages) { for (int stage = 1; stage <= nstages; stage++) { // Set up the STS stage coefficients - // v0 = Y_{n-2} - // v1 = Y_{n-1} + // v0 = Y_{j-1} + // v1 = Y_{j-2} // gam1 = muj = (2.*j - 1.)/j; // gam0 = nuj = (1. - j)/j; // beta_dt/dt = muj_tilde = pm->muj*2./(std::pow(s, 2.) + s); @@ -110,6 +110,8 @@ void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages) { diff_flx = vflx | tflx; } + // TODO Dust diffusion fluxes + // Communicate and set fluxes auto send_flx = tl.AddTask(diff_flx, From 76969348898922f6b96c7e2638ca4aed3c46042f Mon Sep 17 00:00:00 2001 From: Ka Wai HO Date: Mon, 13 Jan 2025 00:37:02 -0700 Subject: [PATCH 04/11] RKL1 implmentation improved and bug fix --- src/CMakeLists.txt | 3 + src/artemis.cpp | 7 ++ src/artemis_driver.cpp | 11 +-- src/gas/gas.cpp | 6 +- src/sts/sts.cpp | 178 ++++++++++++++++++++--------------------- src/sts/sts.hpp | 116 ++++++++++++++++++++++++--- 6 files changed, 208 insertions(+), 113 deletions(-) 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 26eaadac..9d987380 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" @@ -84,6 +85,12 @@ 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), diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index 0e763e68..616ad05e 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -155,15 +155,8 @@ void ArtemisDriver::PreStepTasks() { std::cout << "WARNING: ratio is > 400. Proceed at own risk." << std::endl; } } - auto sts_integrator = artemis_pkg->template Param("sts_integrator"); - if (sts_integrator == STSInt::rkl1){ - // (TODO) RKL1 : Full timestep dt_sts - //STSRKL1(pmesh, tm.time, tm.dt, s_sts); - }else if (sts_integrator == STSInt::rkl2){ - // (TODO) RKL2 : // eq (21) using half hyperbolic timestep - // due to Strang split - //STSRKL2(pmesh, tm.time, 0.5*tm.dt, s_sts); - } + + STS::PreStepSTSTasks(pmesh, tm.time, tm.dt, s_sts); } } diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index cb1b8b03..8378a919 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -478,11 +478,11 @@ Real EstimateTimestepMesh(MeshData *md) { if (do_sts) { const auto sts_max_dt_ratio = params.template Get("sts_max_dt_ratio"); // limit the timestep within the STS ratio, otherwise use the hyperbolic timestep - if (sts_max_dt_ratio > 0.0 && dt_hyp > sts_max_dt_ratio*diff_dt) { - min_dt = std::min(min_dt, max_dt_ratio * diff_dt); + if (sts_max_dt_ratio > 0.0 && min_dt > sts_max_dt_ratio*diff_dt) { + min_dt = std::min(min_dt, sts_max_dt_ratio * diff_dt); } // update the parabolic timestep - gas_pkg->UpdateParam("diff_dt", cfl*diff_dt); + gas_pkg->UpdateParam("diff_dt", cfl_number*diff_dt); }else{ min_dt = std::min(min_dt, diff_dt); } diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp index bbf3994a..4aa414c4 100644 --- a/src/sts/sts.cpp +++ b/src/sts/sts.cpp @@ -20,23 +20,32 @@ // Artemis includes #include "artemis.hpp" +#include "gas/gas.hpp" #include "sts.hpp" #include "utils/artemis_utils.hpp" - namespace STS{ //---------------------------------------------------------------------------------------- //! \fn StateDescriptor STS ::Initialize //! \brief Adds intialization function for STS package std::shared_ptr Initialize(ParameterInput *pin) { - // Extract artemis package - artemis_pkg = pm->packages.Get("artemis").get(); + auto STS = std::make_shared("STS"); + Params ¶ms = STS->AllParams(); // initial check of getting the physics needed for the sts integrator - do_viscosity = artemis_pkg->template Param("do_viscosity"); - do_conduction = artemis_pkg->template Param("do_conduction"); - bool do_diffusion = do_viscosity || do_conduction; + // 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!"); @@ -46,105 +55,47 @@ std::shared_ptr Initialize(ParameterInput *pin) { std::string sts_integrator = pin->GetOrAddString("sts", "integrator", "none"); Real sts_max_dt_ratio = pin->GetOrAddReal("sts","sts_max_dt_ratio", -1.0); - STSInt Diff_integrator = STSInt::null; + STSInt sts_integrator_param = STSInt::null; if (sts_integrator == "rkl1") { - Diff_integrator = STSInt::rkl1; + sts_integrator_param = STSInt::rkl1; } else if (sts_integrator == "rkl2") { PARTHENON_FAIL("rkl2 STS integrator not implemented!"); - Diff_integrator = STSInt::rkl2; + sts_integrator_param = STSInt::rkl2; } else { PARTHENON_FAIL("STS integrator not recognized!"); } - artemis_pkg->AddParam("sts_integrator", sts_integrator); - artemis_pkg->AddParam("sts_max_dt_ratio", sts_max_dt_ratio); + params.Add("sts_integrator", sts_integrator_param); + params.Add("sts_max_dt_ratio", sts_max_dt_ratio); } //---------------------------------------------------------------------------------------- -//! \fn STSRKL1 -//! \brief Assembles the tasks for the STS RKL1 integrator -// comment: Maybe it should be moved back to artemis_driver.cpp -void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages) { - using namespace ::parthenon::Update; - const auto any = parthenon::BoundaryType::any; - const int num_partitions = pmesh->DefaultNumPartitions(); - - // Deep copy u0 into u1 for integrator logic - 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()); - } - - for (int stage = 1; stage <= nstages; stage++) { - // Set up the STS stage coefficients - // v0 = Y_{j-1} - // v1 = Y_{j-2} - // gam1 = muj = (2.*j - 1.)/j; - // gam0 = nuj = (1. - j)/j; - // beta_dt/dt = muj_tilde = pm->muj*2./(std::pow(s, 2.) + s); - Real muj = (1. - stage)/stage; - Real nuj = (2.*stage - 1.)/stage; - Real muj_tilde = (2.*stage - 1.)/stage * 2./(std::pow(nstages, 2.) + nstages); - - 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 = tl.AddTask(none, parthenon::StartReceiveBoundBufs, u0); - auto start_flx_recv = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, u0); - - // Compute (gas) diffusive fluxes - TaskID diff_flx = none; - if (do_diffusion && do_gas) { - 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 Dust diffusion fluxes - - // Communicate and set fluxes - auto send_flx = - tl.AddTask(diff_flx, - parthenon::SendBoundBufs, u0); - auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0); - auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0); - - // Apply flux divergence, STS need stage to 0 for the sts ceofficients - auto update = - tl.AddTask(set_flx, RKL1FluxUpadte, - u0.get(), u1.get(), - dt, muj, nuj, muj_tilde); - } - } -} - -//---------------------------------------------------------------------------------------- -//! \fn STSRKL2FirstStage -//! \brief Assembles the tasks for first stage of the STS RKL2 integrator -void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages) { - // TODO: Implement RKL2 STS integration -} +//! \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_integrator = sts->Param("sts_integrator"); + // Check if the integrator is set + if (sts_integrator == STSInt::null) { + PARTHENON_FAIL("STS integrator not set!"); + } -//---------------------------------------------------------------------------------------- -//! \fn STSRKL2SecondStage -//! \brief Assembles the tasks for first stage of the STS RKL2 integrator -void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages) { - // TODO: Implement RKL2 STS integration + // Execute the integrator tasks + if (sts_integrator == STSInt::rkl1) { + // (TODO) RKL1 : Full timestep dt_sts + STSRKL1(pmesh, time, dt, nstages); + } else if (sts_integrator == STSInt::rkl2) { + // (TODO) RKL2 : // eq (21) using half hyperbolic timestep + // due to Strang split + //STSRKL2FirstStage(pmesh, time, 0.5*dt, nstages); + } } - //---------------------------------------------------------------------------------------- //! \fn RKL1FluxUpadte //! \brief Applies the STS RKL1 update to the conserved variables @@ -208,4 +159,51 @@ TaskStatus RKL1FluxUpadte(MeshData *u0, MeshData *u1, const Real dt, return TaskStatus::complete; } -} \ No newline at end of file +//---------------------------------------------------------------------------------------- +//! \fn STSRKL2FirstStage +//! \brief Assembles the tasks for first stage of the STS RKL2 integrator +template +void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages) { + // TODO: Implement RKL2 STS integration +} + +//---------------------------------------------------------------------------------------- +//! \fn STSRKL2SecondStage +//! \brief Assembles the tasks for first stage of the STS RKL2 integrator +template +void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages) { + // TODO: Implement RKL2 STS integration +} + +//---------------------------------------------------------------------------------------- +//! template instantiations +typedef Mesh M; +//RK1 template instantiations +template void STSRKL1(M *m, const Real time, Real dt, int nstages); +template void STSRKL1( M *m, const Real time, Real dt, int nstages); +template void STSRKL1( M *m, const Real time, Real dt, int nstages); +template void STSRKL1( M *m, const Real time, Real dt, int nstages); +template void STSRKL1( M *m, const Real time, Real dt, int nstages); +template void STSRKL1( M *m, const Real time, Real dt, int nstages); +//RK2 first stage template instantiations +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +//RK2 second stage template instantiations +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +//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 index 181fa52a..3b9dc884 100644 --- a/src/sts/sts.hpp +++ b/src/sts/sts.hpp @@ -14,16 +14,110 @@ #define STS_STS_HPP_ #include "artemis.hpp" +#include "geometry/geometry.hpp" #include "utils/units.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/integrators/artemis_integrator.hpp" -namespace STS { - std::shared_ptr Initialize(ParameterInput *pin); - // task list functions - void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages); - void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages); - void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages); - - // task status functions - TaskStatus RKL1FluxUpadte(MeshData *u0, MeshData *u1, const Real dt, - const Real muj, const Real nuj, const Real muj_tilde); -} \ No newline at end of file +namespace STS{ + +std::shared_ptr Initialize(ParameterInput *pin); + +// STS integrator functions +template +void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int nstages); + +// task list functions +template +void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages); + +template +void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages); +template +void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages); + +// task status functions +template +TaskStatus RKL1FluxUpadte(MeshData *u0, MeshData *u1, const Real dt, + const Real muj, const Real nuj, const Real muj_tilde); + +//---------------------------------------------------------------------------------------- +//! \fn STSRKL1 +//! \brief Assembles the tasks for the STS RKL1 integrator +// comment: Maybe it should be moved back to artemis_driver.cpp +template +void STSRKL1(Mesh *pmesh, const Real time, Real dt, int nstages) { + + using namespace ::parthenon::Update; + TaskCollection tc; + TaskID none(0); + const auto any = parthenon::BoundaryType::any; + const int num_partitions = pmesh->DefaultNumPartitions(); + + auto &pkg = pmesh->packages.Get("STS"); + 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 + 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()); + } + + for (int stage = 1; stage <= nstages; stage++) { + // Set up the STS stage coefficients + // v0 = Y_{j-1} + // v1 = Y_{j-2} + // gam1 = muj = (2.*j - 1.)/j; + // gam0 = nuj = (1. - j)/j; + // beta_dt/dt = muj_tilde = pm->muj*2./(std::pow(s, 2.) + s); + Real muj = (1. - stage)/stage; + Real nuj = (2.*stage - 1.)/stage; + Real muj_tilde = (2.*stage - 1.)/stage * 2./(std::pow(nstages, 2.) + nstages); + + 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 = tl.AddTask(none, parthenon::StartReceiveBoundBufs, u0); + auto start_flx_recv = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, u0); + + // Compute (gas) diffusive fluxes + TaskID diff_flx = none; + if ((do_diffusion) && (do_gas)) { + 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 + + // Communicate and set fluxes + auto send_flx = + tl.AddTask(diff_flx, + parthenon::SendBoundBufs, u0); + auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0); + auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0); + + // Apply flux divergence, STS need stage to 0 for the sts ceofficients + auto update = + tl.AddTask(set_flx, RKL1FluxUpadte, + u0.get(), u1.get(), + dt, muj, nuj, muj_tilde); + } + } +} + +} + +#endif // STS_STS_HPP_ \ No newline at end of file From b81ece763a40ca8d8975155c55df358909d6856d Mon Sep 17 00:00:00 2001 From: Ka Wai HO Date: Mon, 13 Jan 2025 01:11:09 -0700 Subject: [PATCH 05/11] Move the "min_diff_dt" to inner if loop --- src/artemis_driver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index 616ad05e..b7be7a79 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -140,9 +140,9 @@ void ArtemisDriver::PreStepTasks() { auto &u1 = pmesh->mesh_data.Add("u1", u0); // 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"); if (do_sts) { + auto &gas_pkg = pmesh->packages.Get("gas"); + auto min_diff_dt = gas_pkg->template Param("diff_dt"); // compute the number of stages needed for the STS integrator int s_sts = static_cast(0.5 * (std::sqrt(9.0 + 16.0 * tm.dt / min_diff_dt) - 1.0)) + 1; From 3f489d146d3f73a9c47fc204b807ed1909e185ab Mon Sep 17 00:00:00 2001 From: doraemonho Date: Tue, 18 Feb 2025 10:22:13 -0700 Subject: [PATCH 06/11] Workable STS rkl1 Solver Update --- src/artemis_driver.cpp | 55 ++++-- src/artemis_driver.hpp | 1 + src/gas/gas.cpp | 11 +- src/sts/sts.cpp | 170 ++++++++----------- src/sts/sts.hpp | 147 +++++++++------- src/utils/integrators/artemis_integrator.hpp | 30 +++- 6 files changed, 224 insertions(+), 190 deletions(-) diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index b7be7a79..41af5a6f 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -105,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(); + 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,28 +142,36 @@ void ArtemisDriver::PreStepTasks() { auto &base = pmesh->mesh_data.Get(); auto &u0 = pmesh->mesh_data.AddShallow("u0", base, names); 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 - if (do_sts) { - auto &gas_pkg = pmesh->packages.Get("gas"); - auto min_diff_dt = gas_pkg->template Param("diff_dt"); - // compute the number of stages needed for the STS integrator - int s_sts = - static_cast(0.5 * (std::sqrt(9.0 + 16.0 * tm.dt / min_diff_dt) - 1.0)) + 1; - if (s_sts % 2 == 0) s_sts += 1; - - if (parthenon::Globals::my_rank == 0) { - const auto ratio = 2.0 * tm.dt / min_diff_dt; - std::cout << "STS ratio: " << ratio << ", Taking " << s_sts << " steps." << std::endl; - if (ratio > 400.1) { - std::cout << "WARNING: ratio is > 400. Proceed at own risk." << std::endl; - } + 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) { + const auto 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); } + STS::PreStepSTSTasks(pmesh, tm.time, tm.dt, s_sts); + } + //---------------------------------------------------------------------------------------- //! \fn TaskCollection ArtemisDriver::PreStepTasks //! \brief Defines the main integrator's TaskCollection for the ArtemisDriver @@ -237,6 +249,7 @@ 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 && !(do_sts)) { gas_diff_src = tl.AddTask(gas_coord_src | diff_flx | set_flx, Gas::DiffusionUpdate, u0.get(), bdt); @@ -329,31 +342,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 20e839c2..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(); diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index 95212e99..36377cd1 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -188,8 +188,10 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("do_sts", do_sts); if (do_sts) { - Real diff_dt = Big(); - params.Add("diff_dt", diff_dt); + 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; @@ -479,9 +481,10 @@ Real EstimateTimestepMesh(MeshData *md) { 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 && min_dt > sts_max_dt_ratio*diff_dt) { - min_dt = std::min(min_dt, sts_max_dt_ratio * diff_dt); + 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); diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp index 4aa414c4..60a95580 100644 --- a/src/sts/sts.cpp +++ b/src/sts/sts.cpp @@ -24,7 +24,12 @@ #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 @@ -52,22 +57,26 @@ std::shared_ptr Initialize(ParameterInput *pin) { } // Getting the integrator & time ratio between hyperbolic and parabolic terms - std::string sts_integrator = pin->GetOrAddString("sts", "integrator", "none"); + 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_integrator_param = STSInt::null; - if (sts_integrator == "rkl1") { - sts_integrator_param = STSInt::rkl1; - } else if (sts_integrator == "rkl2") { + 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_integrator_param = STSInt::rkl2; + sts_intg_mothod_param = STSInt::rkl2; } else { PARTHENON_FAIL("STS integrator not recognized!"); } - params.Add("sts_integrator", sts_integrator_param); + 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; } //---------------------------------------------------------------------------------------- @@ -78,87 +87,52 @@ 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_integrator = sts->Param("sts_integrator"); + const STSInt sts_intg_mothod = sts->Param("sts_intg_mothod"); // Check if the integrator is set - if (sts_integrator == STSInt::null) { + if (sts_intg_mothod == STSInt::null) { PARTHENON_FAIL("STS integrator not set!"); } // Execute the integrator tasks - if (sts_integrator == STSInt::rkl1) { - // (TODO) RKL1 : Full timestep dt_sts - STSRKL1(pmesh, time, dt, nstages); - } else if (sts_integrator == STSInt::rkl2) { + 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); } } -//---------------------------------------------------------------------------------------- -//! \fn RKL1FluxUpadte -//! \brief Applies the STS RKL1 update to the conserved variables -template -TaskStatus RKL1FluxUpadte(MeshData *u0, MeshData *u1, const Real dt, - const Real muj, const Real nuj, const Real muj_tilde) { - using parthenon::MakePackDescriptor; - using parthenon::variable_names::any; - auto pm = u0->GetParentPointer(); - - // Packing and indexing - std::vector flags({Metadata::Conserved, Metadata::WithFluxes}); - static auto desc = MakePackDescriptor(u0, flags, {parthenon::PDOpt::WithFluxes}); - const auto v0 = desc.GetPack(u0); - const auto v1 = desc.GetPack(u1); - const auto ib = u0->GetBoundsI(IndexDomain::interior); - const auto jb = u0->GetBoundsJ(IndexDomain::interior); - 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, "RKL1FluxUpadte", parthenon::DevExecSpace(), 0, - u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { - // Extract coordinates - using parthenon::TopologicalElement; - geometry::Coords coords(v0.GetCoordinates(b), k, j, i); - - const auto ax1 = coords.GetFaceAreaX1(); - const auto ax2 = (multi_d) ? coords.GetFaceAreaX2() : NewArray(0.0); - const auto ax3 = (three_d) ? coords.GetFaceAreaX3() : NewArray(0.0); - - const Real vol = coords.Volume(); - - for (int n = v0.GetLowerBound(b); n <= v0.GetUpperBound(b); ++n) { - // compute flux divergence - Real divf = (ax1[0] * v0.flux(b, X1DIR, n, k, j, i) - - ax1[1] * v0.flux(b, X1DIR, n, k, j, i + 1)); - if (multi_d) - divf += (ax2[0] * v0.flux(b, X2DIR, n, k, j, i) - - ax2[1] * v0.flux(b, X2DIR, n, k, j + 1, i)); - if (three_d) - divf += (ax3[0] * v0.flux(b, X3DIR, n, k, j, i) - - ax3[1] * v0.flux(b, X3DIR, n, k + 1, j, i)); - - //---------------------------------------------------------------------------------------- - // Apply STS RKL1 update - // Y_{m} = nuj*Y_{j-2} + muj*Y_{j-1} + dt_sts*muj_tilde*M(Y_{j-1}) - // = nuj*Y_{j-2} + muj*Y_{j-1} + dt_sts*muj_tilde*(divf/vol) - // v0 = Y_{j-1}, v1 = Y_{j-2} - Real Y_jm2 = v0(b, n, k, j, i); - v0(b, n, k, j, i) = - nuj * v1(b, n, k, j, i) + muj * v0(b, n, k, j, i) + divf * dt * muj_tilde/ vol; - - // ---------------------------------------------------------------------------------------- - // Rearrange the variables for the next step - v1(b, n, k, j, i) = Y_jm2; - } - }); - return TaskStatus::complete; -} - //---------------------------------------------------------------------------------------- //! \fn STSRKL2FirstStage //! \brief Assembles the tasks for first stage of the STS RKL2 integrator @@ -177,33 +151,27 @@ void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages) { //---------------------------------------------------------------------------------------- //! template instantiations +typedef Coordinates C; typedef Mesh M; -//RK1 template instantiations -template void STSRKL1(M *m, const Real time, Real dt, int nstages); -template void STSRKL1( M *m, const Real time, Real dt, int nstages); -template void STSRKL1( M *m, const Real time, Real dt, int nstages); -template void STSRKL1( M *m, const Real time, Real dt, int nstages); -template void STSRKL1( M *m, const Real time, Real dt, int nstages); -template void STSRKL1( M *m, const Real time, Real dt, int nstages); //RK2 first stage template instantiations -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); //RK2 second stage template instantiations -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); +template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); //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); +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 index 3b9dc884..4c890092 100644 --- a/src/sts/sts.hpp +++ b/src/sts/sts.hpp @@ -14,6 +14,8 @@ #define STS_STS_HPP_ #include "artemis.hpp" +#include "derived/fill_derived.hpp" +#include "gas/gas.hpp" #include "geometry/geometry.hpp" #include "utils/units.hpp" #include "utils/artemis_utils.hpp" @@ -21,35 +23,29 @@ namespace STS{ -std::shared_ptr Initialize(ParameterInput *pin); - -// STS integrator functions -template -void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int nstages); +using Integrator_t = parthenon::LowStorageIntegrator; +using IntegratorPtr_t = std::unique_ptr; -// task list functions -template -void STSRKL1( Mesh *pmesh, const Real time, Real dt, int nstages); +// Global variable declaration +extern IntegratorPtr_t sts_integrator; +std::shared_ptr Initialize(ParameterInput *pin); template void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages); template void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages); -// task status functions -template -TaskStatus RKL1FluxUpadte(MeshData *u0, MeshData *u1, const Real dt, - const Real muj, const Real nuj, const Real muj_tilde); - //---------------------------------------------------------------------------------------- //! \fn STSRKL1 //! \brief Assembles the tasks for the STS RKL1 integrator // comment: Maybe it should be moved back to artemis_driver.cpp template -void STSRKL1(Mesh *pmesh, const Real time, Real dt, int nstages) { - - using namespace ::parthenon::Update; +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(); @@ -61,63 +57,82 @@ void STSRKL1(Mesh *pmesh, const Real time, Real dt, int nstages) { const auto do_gas = pkg->template Param("do_gas"); // Deep copy u0 into u1 for integrator logic - 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()); - } - - for (int stage = 1; stage <= nstages; stage++) { - // Set up the STS stage coefficients - // v0 = Y_{j-1} - // v1 = Y_{j-2} - // gam1 = muj = (2.*j - 1.)/j; - // gam0 = nuj = (1. - j)/j; - // beta_dt/dt = muj_tilde = pm->muj*2./(std::pow(s, 2.) + s); - Real muj = (1. - stage)/stage; - Real nuj = (2.*stage - 1.)/stage; - Real muj_tilde = (2.*stage - 1.)/stage * 2./(std::pow(nstages, 2.) + nstages); - - TaskRegion &tr = tc.AddRegion(num_partitions); + if (stage == 1) { + auto &init_region = tc.AddRegion(num_partitions); for (int i = 0; i < num_partitions; i++) { - auto &tl = tr[i]; + auto &tl = init_region[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 = tl.AddTask(none, parthenon::StartReceiveBoundBufs, u0); - auto start_flx_recv = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, u0); - - // Compute (gas) diffusive fluxes - TaskID diff_flx = none; - if ((do_diffusion) && (do_gas)) { - 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 - - // Communicate and set fluxes - auto send_flx = - tl.AddTask(diff_flx, - parthenon::SendBoundBufs, u0); - auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0); - auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0); - - // Apply flux divergence, STS need stage to 0 for the sts ceofficients - auto update = - tl.AddTask(set_flx, RKL1FluxUpadte, - u0.get(), u1.get(), - dt, muj, nuj, muj_tilde); + 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()); + + // Apply "coordinate source terms" + TaskID gas_coord_src = swap_data_1, dust_coord_src = swap_data_1; + if (do_gas) gas_coord_src = tl.AddTask(swap_data_1, Gas::FluxSource, u0.get(), dt); + + TaskID gas_diff_src = gas_coord_src | diff_flx | set_flx_u0; + gas_diff_src = tl.AddTask( gas_coord_src | 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); + } #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..a28d7390 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, @@ -105,6 +132,7 @@ TaskStatus ApplyUpdate(MeshData *u0, MeshData *u1, const int stage, v0(b, n, k, j, i) = gam0 * v0(b, n, k, j, i) + gam1 * v1(b, n, k, j, i) + divf * beta_dt / vol; } + }); return TaskStatus::complete; } From d0f484146cff01c3132a28d3568365ac77f5c643 Mon Sep 17 00:00:00 2001 From: doraemonho Date: Tue, 18 Feb 2025 10:51:36 -0700 Subject: [PATCH 07/11] Format Fix --- src/artemis.cpp | 6 +- src/artemis_driver.cpp | 12 ++-- src/gas/gas.cpp | 15 +++-- src/sts/sts.cpp | 66 +++++++------------- src/sts/sts.hpp | 38 ++++++----- src/utils/integrators/artemis_integrator.hpp | 5 +- 6 files changed, 56 insertions(+), 86 deletions(-) diff --git a/src/artemis.cpp b/src/artemis.cpp index 670cab67..f02a2c28 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -85,18 +85,18 @@ 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)), + 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), "Viscosity requires the gas package, but there is not gas!"); PARTHENON_REQUIRE(!(do_conduction) || (do_conduction && do_gas), - "Conduction requires the gas package, but there is not gas!"); + "Conduction requires the gas package, but there is not gas!"); PARTHENON_REQUIRE(!(do_radiation) || (do_radiation && do_gas), "Radiation requires the gas package, but there is not gas!"); diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index 41af5a6f..140ab96a 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -149,26 +149,24 @@ void ArtemisDriver::PreStepTasks() { //! \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))); + 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) { - const auto ratio = tm.dt / min_diff_dt; + 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); - } @@ -249,7 +247,7 @@ 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 && !(do_sts)) { gas_diff_src = tl.AddTask(gas_coord_src | diff_flx | set_flx, Gas::DiffusionUpdate, u0.get(), bdt); diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index 36377cd1..1cb91e68 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -188,9 +188,8 @@ std::shared_ptr Initialize(ParameterInput *pin, 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("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); } @@ -476,7 +475,7 @@ Real EstimateTimestepMesh(MeshData *md) { Real diff_dt = std::min(visc_dt, cond_dt); const auto cfl_number = params.template Get("cfl"); - + // STS Time Stepping Control const auto do_sts = params.template Get("do_sts"); if (do_sts) { @@ -484,15 +483,15 @@ Real EstimateTimestepMesh(MeshData *md) { 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; + min_dt = sts_max_dt_ratio*diff_dt; } // update the parabolic timestep - gas_pkg->UpdateParam("diff_dt", cfl_number*diff_dt); - }else{ + gas_pkg->UpdateParam("diff_dt", cfl_number * diff_dt); + } else { min_dt = std::min(min_dt, diff_dt); } - return cfl_number*min_dt; + return cfl_number * min_dt; } //---------------------------------------------------------------------------------------- diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp index 60a95580..f062c79f 100644 --- a/src/sts/sts.cpp +++ b/src/sts/sts.cpp @@ -58,7 +58,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { // 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); + 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; @@ -75,7 +75,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { 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; } @@ -100,7 +100,7 @@ void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int nstages) { 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}), + // 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 @@ -114,11 +114,11 @@ void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int nstages) { // 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 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 @@ -127,51 +127,27 @@ void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int nstages) { } } else if (sts_intg_mothod == STSInt::rkl2) { PARTHENON_FAIL("STS rkl2 integrator not implemented!"); - // (TODO) RKL2 : // eq (21) using half hyperbolic timestep + // (TODO) RKL2 : // eq (21) using half hyperbolic timestep // due to Strang split - //STSRKL2FirstStage(pmesh, time, 0.5*dt, nstages); + // STSRKL2FirstStage(pmesh, time, 0.5*dt, nstages); } } -//---------------------------------------------------------------------------------------- -//! \fn STSRKL2FirstStage -//! \brief Assembles the tasks for first stage of the STS RKL2 integrator -template -void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages) { - // TODO: Implement RKL2 STS integration -} - -//---------------------------------------------------------------------------------------- -//! \fn STSRKL2SecondStage -//! \brief Assembles the tasks for first stage of the STS RKL2 integrator -template -void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages) { - // TODO: Implement RKL2 STS integration -} - //---------------------------------------------------------------------------------------- //! template instantiations typedef Coordinates C; typedef Mesh M; -//RK2 first stage template instantiations -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2FirstStage(M *m, const Real time, Real dt, int nstages); -//RK2 second stage template instantiations -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); -template void STSRKL2SecondStage(M *m, const Real time, Real dt, int nstages); //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); +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 index 4c890092..d8fa813f 100644 --- a/src/sts/sts.hpp +++ b/src/sts/sts.hpp @@ -21,7 +21,7 @@ #include "utils/artemis_utils.hpp" #include "utils/integrators/artemis_integrator.hpp" -namespace STS{ +namespace STS { using Integrator_t = parthenon::LowStorageIntegrator; using IntegratorPtr_t = std::unique_ptr; @@ -30,10 +30,6 @@ using IntegratorPtr_t = std::unique_ptr; extern IntegratorPtr_t sts_integrator; std::shared_ptr Initialize(ParameterInput *pin); -template -void STSRKL2FirstStage( Mesh *pm, const Real time, Real dt, int nstages); -template -void STSRKL2SecondStage( Mesh *pm, const Real time, Real dt, int nstages); //---------------------------------------------------------------------------------------- //! \fn STSRKL1 @@ -66,13 +62,13 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst 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); @@ -88,16 +84,16 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst // 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); + 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()); + 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()); @@ -107,22 +103,24 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst if (do_gas) gas_coord_src = tl.AddTask(swap_data_1, Gas::FluxSource, u0.get(), dt); TaskID gas_diff_src = gas_coord_src | diff_flx | set_flx_u0; - gas_diff_src = tl.AddTask( gas_coord_src | diff_flx | set_flx_u0, + gas_diff_src = tl.AddTask(gas_coord_src | 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()); + 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, + 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); + 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()); + auto c2p_u0 = + tl.AddTask(TQ::local_sync, bcs_u0, FillDerived>, u0.get()); } @@ -133,6 +131,6 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst 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 a28d7390..915823e6 100644 --- a/src/utils/integrators/artemis_integrator.hpp +++ b/src/utils/integrators/artemis_integrator.hpp @@ -66,8 +66,8 @@ inline TaskStatus SwapData(MeshData *u0, MeshData *u1) { 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, + 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); @@ -132,7 +132,6 @@ TaskStatus ApplyUpdate(MeshData *u0, MeshData *u1, const int stage, v0(b, n, k, j, i) = gam0 * v0(b, n, k, j, i) + gam1 * v1(b, n, k, j, i) + divf * beta_dt / vol; } - }); return TaskStatus::complete; } From 2209813c2ed2a8a744ac5f27f446bc91fa4887a3 Mon Sep 17 00:00:00 2001 From: doraemonho Date: Tue, 18 Feb 2025 11:04:22 -0700 Subject: [PATCH 08/11] Format Fix --- src/artemis.cpp | 3 +-- src/artemis_driver.cpp | 6 +++--- src/gas/gas.cpp | 2 +- src/sts/sts.cpp | 17 ++++++++--------- src/sts/sts.hpp | 13 ++++++------- 5 files changed, 19 insertions(+), 22 deletions(-) diff --git a/src/artemis.cpp b/src/artemis.cpp index f02a2c28..ffab8fd2 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -88,7 +88,6 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { 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), @@ -96,7 +95,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { PARTHENON_REQUIRE(!(do_viscosity) || (do_viscosity && do_gas), "Viscosity requires the gas package, but there is not gas!"); PARTHENON_REQUIRE(!(do_conduction) || (do_conduction && do_gas), - "Conduction requires the gas package, but there is not gas!"); + "Conduction requires the gas package, but there is not gas!"); PARTHENON_REQUIRE(!(do_radiation) || (do_radiation && do_gas), "Radiation requires the gas package, but there is not gas!"); diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index 140ab96a..8d5549bf 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -29,8 +29,8 @@ #include "nbody/nbody.hpp" #include "radiation/imc/imc.hpp" #include "rotating_frame/rotating_frame.hpp" -#include "utils/integrators/artemis_integrator.hpp" #include "sts/sts.hpp" +#include "utils/integrators/artemis_integrator.hpp" using namespace parthenon::driver::prelude; @@ -109,9 +109,9 @@ TaskListStatus ArtemisDriver::Step() { if (do_sts) STSFirstStage(); // Execute explicit, unsplit physics - auto status = StepTasks().Execute(); + 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); diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index 1cb91e68..bbedfe94 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -483,7 +483,7 @@ Real EstimateTimestepMesh(MeshData *md) { 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; + min_dt = sts_max_dt_ratio * diff_dt; } // update the parabolic timestep gas_pkg->UpdateParam("diff_dt", cfl_number * diff_dt); diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp index f062c79f..b0d1d66f 100644 --- a/src/sts/sts.cpp +++ b/src/sts/sts.cpp @@ -26,7 +26,7 @@ using ArtemisUtils::VI; -namespace STS{ +namespace STS { IntegratorPtr_t sts_integrator; @@ -137,17 +137,16 @@ void PreStepSTSTasks(Mesh *pmesh, const Real time, Real dt, int 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, +// 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, +template void PreStepSTSTasks(M *m, const Real time, Real dt, int nstages); -template void PreStepSTSTasks(M *m, const Real time, Real dt, +template void PreStepSTSTasks(M *m, const Real time, Real dt, int nstages); -template void PreStepSTSTasks(M *m, const Real time, Real dt, +template void PreStepSTSTasks(M *m, const Real time, Real dt, int nstages); -template void PreStepSTSTasks(M *m,const Real time, Real dt, +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 index d8fa813f..d9812646 100644 --- a/src/sts/sts.hpp +++ b/src/sts/sts.hpp @@ -17,9 +17,9 @@ #include "derived/fill_derived.hpp" #include "gas/gas.hpp" #include "geometry/geometry.hpp" -#include "utils/units.hpp" #include "utils/artemis_utils.hpp" #include "utils/integrators/artemis_integrator.hpp" +#include "utils/units.hpp" namespace STS { @@ -85,8 +85,8 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst // Communicate and set fluxes auto send_flx = tl.AddTask( - diff_flx, parthenon::SendBoundBufs, u0); - auto recv_flx_u0 = + 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); @@ -94,7 +94,7 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst 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()); @@ -115,13 +115,12 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst PreCommFillDerived>, u0.get()); // Set boundary conditions (both physical and logical) - auto bcs_u0 = + auto bcs_u0 = parthenon::AddBoundaryExchangeTasks(pre_comm_u0, tl, u0, pmesh->multilevel); // Update primitive variables - auto c2p_u0 = + auto c2p_u0 = tl.AddTask(TQ::local_sync, bcs_u0, FillDerived>, u0.get()); - } return tc; From 5dd332ad41386ad432557d13d3b1e095649e162a Mon Sep 17 00:00:00 2001 From: doraemonho Date: Tue, 18 Feb 2025 11:09:20 -0700 Subject: [PATCH 09/11] Format Fix --- src/artemis_driver.cpp | 1 - src/sts/sts.cpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index 8d5549bf..358bfb39 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -169,7 +169,6 @@ void ArtemisDriver::STSFirstStage() { STS::PreStepSTSTasks(pmesh, tm.time, tm.dt, s_sts); } - //---------------------------------------------------------------------------------------- //! \fn TaskCollection ArtemisDriver::PreStepTasks //! \brief Defines the main integrator's TaskCollection for the ArtemisDriver diff --git a/src/sts/sts.cpp b/src/sts/sts.cpp index b0d1d66f..4bc52d77 100644 --- a/src/sts/sts.cpp +++ b/src/sts/sts.cpp @@ -147,6 +147,6 @@ 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, +template void PreStepSTSTasks(M *m, const Real time, Real dt, int nstages); } // namespace STS \ No newline at end of file From 979e35d34a89a8eef92d6b0f95089881ccadcea5 Mon Sep 17 00:00:00 2001 From: doraemonho Date: Tue, 25 Feb 2025 09:43:19 -0700 Subject: [PATCH 10/11] STS test problem update - Bug fixed for sts.hpp - changed duffusion problem py script to sts solver - changed the conduction input file --- inputs/diffusion/conduction.in | 5 +++++ src/sts/sts.hpp | 10 +++------- tst/scripts/diffusion/thermal_diffusion.py | 3 ++- 3 files changed, 10 insertions(+), 8 deletions(-) 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/sts/sts.hpp b/src/sts/sts.hpp index d9812646..6eccbe9e 100644 --- a/src/sts/sts.hpp +++ b/src/sts/sts.hpp @@ -46,7 +46,7 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst const auto any = parthenon::BoundaryType::any; const int num_partitions = pmesh->DefaultNumPartitions(); - auto &pkg = pmesh->packages.Get("STS"); + 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"); @@ -98,12 +98,8 @@ TaskCollection STSRKL1(Mesh *pmesh, const Real time, Real dt, int stage, int nst // swap u0 <-> u1 auto swap_data_1 = tl.AddTask(update, ArtemisUtils::SwapData, u0.get(), u1.get()); - // Apply "coordinate source terms" - TaskID gas_coord_src = swap_data_1, dust_coord_src = swap_data_1; - if (do_gas) gas_coord_src = tl.AddTask(swap_data_1, Gas::FluxSource, u0.get(), dt); - - TaskID gas_diff_src = gas_coord_src | diff_flx | set_flx_u0; - gas_diff_src = tl.AddTask(gas_coord_src | diff_flx | set_flx_u0, + 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 diff --git a/tst/scripts/diffusion/thermal_diffusion.py b/tst/scripts/diffusion/thermal_diffusion.py index 15c4ac60..1de1e9f9 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), @@ -96,7 +97,7 @@ def analyze(): "final", dir=artemis.get_data_dir(), base="{}.out1".format(name) ) xc = 0.5 * (x[1:] + x[:-1]) - ans = Tans(xc.ravel(), f=_flux, T0=_gtemp, x0=1.2, xi=0.2, d=dind[g], k=_kcond) + ans = Tans(xc.ravel(), f=_flux, T0=_gtemp, x0=1.2, xi=0.2, d=())[g], k=_kcond) temp = T[0, :].ravel() err = abs(temp / ans - 1.0) ax.plot(xc, ans, "--k") From f01b76514126cc0104275053b02a4ff4b3ffc5fd Mon Sep 17 00:00:00 2001 From: doraemonho Date: Tue, 25 Feb 2025 09:55:26 -0700 Subject: [PATCH 11/11] typo fixed in thermal_duffusion.py --- tst/scripts/diffusion/thermal_diffusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tst/scripts/diffusion/thermal_diffusion.py b/tst/scripts/diffusion/thermal_diffusion.py index 1de1e9f9..e721da82 100644 --- a/tst/scripts/diffusion/thermal_diffusion.py +++ b/tst/scripts/diffusion/thermal_diffusion.py @@ -97,7 +97,7 @@ def analyze(): "final", dir=artemis.get_data_dir(), base="{}.out1".format(name) ) xc = 0.5 * (x[1:] + x[:-1]) - ans = Tans(xc.ravel(), f=_flux, T0=_gtemp, x0=1.2, xi=0.2, d=())[g], k=_kcond) + ans = Tans(xc.ravel(), f=_flux, T0=_gtemp, x0=1.2, xi=0.2, d=dind[g], k=_kcond) temp = T[0, :].ravel() err = abs(temp / ans - 1.0) ax.plot(xc, ans, "--k")