diff --git a/CHANGELOG.md b/CHANGELOG.md index 152bd814..b1d49a38 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping ### Changed (changing behavior/API/variables/...) +- [[PR 122]](https://github.com/parthenon-hpc-lab/athenapk/pull/122) Fixed sqrt(4pi) factor in CGS Gauss unit and add unit doc - [[PR 119]](https://github.com/parthenon-hpc-lab/athenapk/pull/119) Fixed Athena++ paper test case for KHI pgen. Added turbulence pgen doc. - [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97) Fixed Schure cooling curve. Removed SD one. Added description of cooling function conventions. - [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15) @@ -14,6 +15,7 @@ ### Fixed (not changing behavior/API/variables/...) ### Infrastructure +- [[PR 124]](https://github.com/parthenon-hpc-lab/athenapk/pull/124) Bump Kokkos 4.4.1 (and Parthenon to include view-of-view fix) - [[PR 117]](https://github.com/parthenon-hpc-lab/athenapk/pull/117) Update devcontainer.json to latest CI container - [[PR 114]](https://github.com/parthenon-hpc-lab/athenapk/pull/114) Bump Parthenon 24.08 and Kokkos to 4.4.00 - [[PR 112]](https://github.com/parthenon-hpc-lab/athenapk/pull/112) Add dev container configuration @@ -23,6 +25,8 @@ ### Removed (removing behavior/API/varaibles/...) ### Incompatibilities (i.e. breaking changes) +- [[PR 124]](https://github.com/parthenon-hpc-lab/athenapk/pull/124) Enrolling custom boundary conditions changed + - Boundary conditions can now be enrolled using a string that can be subsequently be used in the input file (see, e.g., cloud problem generator) - [[PR 114]](https://github.com/parthenon-hpc-lab/athenapk/pull/114) Bump Parthenon 24.08 and Kokkos to 4.4.00 - Changed signature of `UserWorkBeforeOutput` to include `SimTime` as last paramter - Fixes bitwise idential restarts for AMR simulations (the derefinement counter is now included) diff --git a/README.md b/README.md index fd66e824..882d9625 100644 --- a/README.md +++ b/README.md @@ -131,21 +131,8 @@ the `file_type = hdf5` format, see [VisIt](https://wci.llnl.gov/simulation/computer-codes/visit/). In ParaView, select the "XDMF Reader" when prompted. -2. With [yt](https://yt-project.org/) -- though currently through a custom frontend -that is not yet part of the main yt branch and, thus, has to be installed manually, e.g., -as follows: -```bash -cd ~/src # or any other folder of choice -git clone https://github.com/forrestglines/yt.git -cd yt -git checkout parthenon-frontend - -# If you're using conda or virtualenv -pip install -e . -# OR alternatively, if you using the plain Python environment -pip install --user -e . -``` -Afterwards, `*.phdf` files can be read as usual with `yt.load()`. +2. With [yt](https://yt-project.org/) +As of versions >=4.4 `*.phdf` files can be read as usual with `yt.load()`. 3. Using [Ascent](https://github.com/Alpine-DAV/ascent) (for in situ visualization and analysis). This requires Ascent to be installed/available at compile time of AthenaPK. diff --git a/docs/README.md b/docs/README.md index 018d4b68..145cc0f8 100644 --- a/docs/README.md +++ b/docs/README.md @@ -13,6 +13,7 @@ The documentation currently includes - [Notebooks to calculate cooling tables from literature](cooling) - [Brief notes on developing code for AthenaPK](development.md) - [How to add a custom/user problem generator](pgen.md) +- [Units](units.md) - Detailed descriptions of more complex problem generators - [Galaxy Cluster and Cluster-like Problem Setup](cluster.md) - [Driven turbulence](turbulence.md) diff --git a/docs/input.md b/docs/input.md index d5e48d4f..0395a448 100644 --- a/docs/input.md +++ b/docs/input.md @@ -67,6 +67,10 @@ To control the floors, following parameters can be set in the `` block: *Note* the pressure floor will take precedence over the temperature floor in the conserved to primitive conversion if both are defined. +#### Units + +See(here)[units.md]. + #### Diffusive processes Diffusive processes in AthenaPK can be configured in the `` block of the input file. @@ -114,9 +118,6 @@ However, as reported by [^V+17] a large number of stages (in combination with anisotropic, limited diffusion) may lead to a loss in accuracy, which is why the difference between hyperbolic and parabolic timesteps can be limited by `diffusion/rkl2_max_dt_ratio=...` and a warning is shown if the ratio is above 400. -Note that if this limit is enforced the `dt=` shown on the terminal might not be the actual -`dt` taken in the code as the limit is currently enforced only after the output -has been printed. [^M+14]: C. D. Meyer, D. S. Balsara, and T. D. Aslam, “A stabilized Runge–Kutta–Legendre method for explicit super-time-stepping of parabolic and mixed equations,” Journal of Computational Physics, vol. 257, pp. 594–626, 2014, doi: https://doi.org/10.1016/j.jcp.2013.08.021. diff --git a/docs/units.md b/docs/units.md new file mode 100644 index 00000000..b26ccf6e --- /dev/null +++ b/docs/units.md @@ -0,0 +1,58 @@ +# AthenaPK units + +## General unit system + +Internally, all calculations are done in "code units" and there are no conversions between +code and physical units during runtime (with excetion like the temperature when cooling is +being used). +Therefore, in general no units need to be prescribed to run a simulation. + +If units are required (e.g., if cooling is used and, thus, a conversion between internal energy +in code units and physical temperature is required) they are configured in the input block +as follows: + +``` + +code_length_cgs = 3.085677580962325e+24 # 1 Mpc in cm +code_mass_cgs = 1.98841586e+47 # 1e14 Msun in g +code_time_cgs = 3.15576e+16 # 1 Gyr in s +``` + +This information will also be used by postprocessing tools (like yt) to convert between +code units and a physical unit system (like cgs). + +Moreover, internally a set of factors from code to cgs units are available to process conversions +if required (e.g., from the input file). + +For example, for an input parameter (in the input file) like + +``` + +r0_cgs = 3.085677580962325e+20 # 100 pc +``` + +the conversion should happen in the problem generator lik + +```c++ + r_cloud = pin->GetReal("problem/cloud", "r0_cgs") / units.code_length_cgs(); +``` + +so that the resulting quantity is internally in code units (here code length). + +It is highly recommended to be *very* explicit/specific about units everywhere (as it is +a common source of confusion) like adding the `_cgs` suffix to the parameter in the +input file above. + +## Magnetic units + +Internally, AthenaPK (and almost all MHD codes) use +[Heaviside-Lorentz units](https://en.wikipedia.org/wiki/Heaviside%E2%80%93Lorentz_units), +where the magnetic field is transformed from $B \rightarrow B / \sqrt{4 \pi}$. +(See also the note in the +[Castro documentation](https://amrex-astro.github.io/Castro/docs/mhd.html) about this.) + +So when converting from CGS-Gaussian units to code units, it is necessary to divide +by $\sqrt{4 \pi}$ (in addition to the base dimensional factors). +This is automatically handled by the `units.code_magnetic_cgs()` factors. + + diff --git a/external/Kokkos b/external/Kokkos index 08ceff92..15dc143e 160000 --- a/external/Kokkos +++ b/external/Kokkos @@ -1 +1 @@ -Subproject commit 08ceff92bcf3a828844480bc1e6137eb74028517 +Subproject commit 15dc143e5f39949eece972a798e175c4b463d4b8 diff --git a/external/parthenon b/external/parthenon index ec61c9cb..e00017c6 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit ec61c9cb102d6a67a4dd03e7d7fa4e10c82c1032 +Subproject commit e00017c6ece547f3adf9a16b94ddf0bdf8bb4773 diff --git a/inputs/cloud.in b/inputs/cloud.in index 002febd0..e7e05f69 100644 --- a/inputs/cloud.in +++ b/inputs/cloud.in @@ -22,7 +22,7 @@ ox1_bc = outflow nx2 = 320 x2min = -400 x2max = 2100 -ix2_bc = user +ix2_bc = cloud_inflow_x2 ox2_bc = outflow nx3 = 192 diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 7bc95be7..5720346c 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -4,6 +4,7 @@ // Licensed under the BSD 3-Clause License (the "LICENSE"). //======================================================================================== +#include #include #include #include @@ -29,6 +30,7 @@ #include "diffusion/diffusion.hpp" #include "glmmhd/glmmhd.hpp" #include "hydro.hpp" +#include "interface/params.hpp" #include "outputs/outputs.hpp" #include "prolongation/custom_ops.hpp" #include "rsolvers/rsolvers.hpp" @@ -60,44 +62,7 @@ parthenon::Packages_t ProcessPackages(std::unique_ptr &pin) { // the task list is constructed (versus when the task list is being executed). // TODO(next person touching this function): If more/separate feature are required // please separate concerns. -void PreStepMeshUserWorkInLoop(Mesh *pmesh, ParameterInput *pin, SimTime &tm) { - auto hydro_pkg = pmesh->block_list[0]->packages.Get("Hydro"); - const auto num_partitions = pmesh->DefaultNumPartitions(); - - if (hydro_pkg->Param("diffint") == DiffInt::rkl2) { - auto dt_diff = std::numeric_limits::max(); - if (hydro_pkg->Param("conduction") != Conduction::none) { - for (auto i = 0; i < num_partitions; i++) { - auto &md = pmesh->mesh_data.GetOrAdd("base", i); - - dt_diff = std::min(dt_diff, EstimateConductionTimestep(md.get())); - } - } - if (hydro_pkg->Param("viscosity") != Viscosity::none) { - for (auto i = 0; i < num_partitions; i++) { - auto &md = pmesh->mesh_data.GetOrAdd("base", i); - - dt_diff = std::min(dt_diff, EstimateViscosityTimestep(md.get())); - } - } - if (hydro_pkg->Param("resistivity") != Resistivity::none) { - for (auto i = 0; i < num_partitions; i++) { - auto &md = pmesh->mesh_data.GetOrAdd("base", i); - - dt_diff = std::min(dt_diff, EstimateResistivityTimestep(md.get())); - } - } -#ifdef MPI_PARALLEL - PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &dt_diff, 1, MPI_PARTHENON_REAL, - MPI_MIN, MPI_COMM_WORLD)); -#endif - hydro_pkg->UpdateParam("dt_diff", dt_diff); - const auto max_dt_ratio = hydro_pkg->Param("rkl2_max_dt_ratio"); - if (max_dt_ratio > 0.0 && tm.dt / dt_diff > max_dt_ratio) { - tm.dt = max_dt_ratio * dt_diff; - } - } -} +void PreStepMeshUserWorkInLoop(Mesh *pmesh, ParameterInput *pin, SimTime &tm) {} template Real HydroHst(MeshData *md) { @@ -252,17 +217,23 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto glmmhd_alpha = pin->GetOrAddReal("hydro", "glmmhd_alpha", 0.1); pkg->AddParam("glmmhd_alpha", glmmhd_alpha); calc_c_h = true; - pkg->AddParam("c_h", 0.0, true); // hyperbolic divergence cleaning speed - // global minimum dx (used to calc c_h) - pkg->AddParam("mindx", std::numeric_limits::max(), true); - // hyperbolic timestep constraint - pkg->AddParam("dt_hyp", std::numeric_limits::max(), true); } else { PARTHENON_FAIL("AthenaPK hydro: Unknown fluid method."); } pkg->AddParam<>("fluid", fluid); pkg->AddParam<>("nhydro", nhydro); pkg->AddParam<>("calc_c_h", calc_c_h); + // Following params should (currently) be present independent of solver because + // they're all used in the main loop. + // TODO(pgrete) think about which approach (selective versus always is preferable) + pkg->AddParam( + "c_h", 0.0, Params::Mutability::Restart); // hyperbolic divergence cleaning speed + // global minimum dx (used to calc c_h) + pkg->AddParam("mindx", std::numeric_limits::max(), + Params::Mutability::Restart); + // hyperbolic timestep constraint + pkg->AddParam("dt_hyp", std::numeric_limits::max(), + Params::Mutability::Restart); const auto recon_str = pin->GetString("hydro", "reconstruction"); int recon_need_nghost = 3; // largest number for the choices below @@ -633,12 +604,13 @@ std::shared_ptr Initialize(ParameterInput *pin) { "Options are: none, unsplit, rkl2"); } if (diffint != DiffInt::none) { - pkg->AddParam("dt_diff", 0.0, true); // diffusive timestep constraint // As in Athena++ a cfl safety factor is also applied to the theoretical limit. // By default it is equal to the hyperbolic cfl. auto cfl_diff = pin->GetOrAddReal("diffusion", "cfl", pkg->Param("cfl")); pkg->AddParam<>("cfl_diff", cfl_diff); } + pkg->AddParam("dt_diff", std::numeric_limits::max(), + Params::Mutability::Restart); // diffusive timestep constraint pkg->AddParam<>("diffint", diffint); if (fluid == Fluid::euler) { @@ -855,9 +827,12 @@ Real EstimateTimestep(MeshData *md) { // get to package via first block in Meshdata (which exists by construction) auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); auto min_dt = std::numeric_limits::max(); + auto dt_hyp = std::numeric_limits::max(); - if (hydro_pkg->Param("calc_dt_hyp")) { - min_dt = std::min(min_dt, EstimateHyperbolicTimestep(md)); + const auto calc_dt_hyp = hydro_pkg->Param("calc_dt_hyp"); + if (calc_dt_hyp) { + dt_hyp = EstimateHyperbolicTimestep(md); + min_dt = std::min(min_dt, dt_hyp); } const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); @@ -869,17 +844,34 @@ Real EstimateTimestep(MeshData *md) { min_dt = std::min(min_dt, tabular_cooling.EstimateTimeStep(md)); } - // For RKL2 STS, the diffusive timestep is calculated separately in the driver - if (hydro_pkg->Param("diffint") == DiffInt::unsplit) { + auto dt_diff = std::numeric_limits::max(); + if (hydro_pkg->Param("diffint") != DiffInt::none) { if (hydro_pkg->Param("conduction") != Conduction::none) { - min_dt = std::min(min_dt, EstimateConductionTimestep(md)); + dt_diff = std::min(dt_diff, EstimateConductionTimestep(md)); } if (hydro_pkg->Param("viscosity") != Viscosity::none) { - min_dt = std::min(min_dt, EstimateViscosityTimestep(md)); + dt_diff = std::min(dt_diff, EstimateViscosityTimestep(md)); } if (hydro_pkg->Param("resistivity") != Resistivity::none) { - min_dt = std::min(min_dt, EstimateResistivityTimestep(md)); + dt_diff = std::min(dt_diff, EstimateResistivityTimestep(md)); + } + + // For unsplit ingegration use strict limit + if (hydro_pkg->Param("diffint") == DiffInt::unsplit) { + min_dt = std::min(min_dt, dt_diff); + // and for RKL2 integration use limit taking into account the maxium ratio + // or not constrain limit further (which is why RKL2 is there in first place) + } else if (hydro_pkg->Param("diffint") == DiffInt::rkl2) { + const auto max_dt_ratio = hydro_pkg->Param("rkl2_max_dt_ratio"); + if (max_dt_ratio > 0.0 && dt_hyp / dt_diff > max_dt_ratio) { + min_dt = std::min(min_dt, max_dt_ratio * dt_diff); + } + } else { + PARTHENON_THROW("Looks like a a new diffusion integrator was implemented without " + "taking into accout timestep contstraints yet."); } + auto dt_diff_param = hydro_pkg->Param("dt_diff"); + hydro_pkg->UpdateParam("dt_diff", std::min(dt_diff, dt_diff_param)); } if (ProblemEstimateTimestep != nullptr) { diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index 36bdce91..eeb29d3d 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -447,7 +447,10 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { // Calculate hyperbolic divergence cleaning speed // TODO(pgrete) Calculating mindx is only required after remeshing. Need to find a clean // solution for this one-off global reduction. - if (hydro_pkg->Param("calc_c_h") && (stage == 1)) { + // TODO(PG) move this to PreStepMeshUserWorkInLoop + if ((hydro_pkg->Param("calc_c_h") || + hydro_pkg->Param("diffint") != DiffInt::none) && + (stage == 1)) { // need to make sure that there's only one region in order to MPI_reduce to work TaskRegion &single_task_region = tc.AddRegion(1); auto &tl = single_task_region[0]; @@ -468,14 +471,16 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { reduce_c_h = tl.AddTask( prev_task, [](StateDescriptor *hydro_pkg) { - Real mins[2]; + Real mins[3]; mins[0] = hydro_pkg->Param("mindx"); mins[1] = hydro_pkg->Param("dt_hyp"); - PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, mins, 2, MPI_PARTHENON_REAL, + mins[2] = hydro_pkg->Param("dt_diff"); + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, mins, 3, MPI_PARTHENON_REAL, MPI_MIN, MPI_COMM_WORLD)); hydro_pkg->UpdateParam("mindx", mins[0]); hydro_pkg->UpdateParam("dt_hyp", mins[1]); + hydro_pkg->UpdateParam("dt_diff", mins[2]); return TaskStatus::complete; }, hydro_pkg.get()); @@ -657,7 +662,9 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { // Single task in single (serial) region to reset global vars used in reductions in the // first stage. - if (stage == integrator->nstages && hydro_pkg->Param("calc_c_h")) { + if (stage == integrator->nstages && + (hydro_pkg->Param("calc_c_h") || + hydro_pkg->Param("diffint") != DiffInt::none)) { TaskRegion &reset_reduction_vars_region = tc.AddRegion(1); auto &tl = reset_reduction_vars_region[0]; tl.AddTask( @@ -665,6 +672,7 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { [](StateDescriptor *hydro_pkg) { hydro_pkg->UpdateParam("mindx", std::numeric_limits::max()); hydro_pkg->UpdateParam("dt_hyp", std::numeric_limits::max()); + hydro_pkg->UpdateParam("dt_diff", std::numeric_limits::max()); return TaskStatus::complete; }, hydro_pkg.get()); diff --git a/src/main.cpp b/src/main.cpp index f239b197..23cdb825 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -63,8 +63,8 @@ int main(int argc, char *argv[]) { } else if (problem == "cloud") { pman.app_input->InitUserMeshData = cloud::InitUserMeshData; pman.app_input->ProblemGenerator = cloud::ProblemGenerator; - pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x2] = - cloud::InflowWindX2; + pman.app_input->RegisterBoundaryCondition(parthenon::BoundaryFace::inner_x2, + "cloud_inflow_x2", cloud::InflowWindX2); Hydro::ProblemCheckRefinementBlock = cloud::ProblemCheckRefinementBlock; } else if (problem == "blast") { pman.app_input->InitUserMeshData = blast::InitUserMeshData; diff --git a/src/units.hpp b/src/units.hpp index 3290e87e..fc6a2412 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -90,7 +90,8 @@ class Units { return code_energy_cgs() / kev_cgs * code_length_cgs() * code_length_cgs(); } parthenon::Real code_magnetic_cgs() const { - return sqrt(code_mass_cgs()) / sqrt(code_length_cgs()) / code_time_cgs(); + return std::sqrt(4.0 * M_PI) * sqrt(code_mass_cgs()) / sqrt(code_length_cgs()) / + code_time_cgs(); } // Physical Constants in code units