diff --git a/CHANGELOG.md b/CHANGELOG.md index 65d55b73..152bd814 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 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) diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 00000000..018d4b68 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,36 @@ +# AthenaPK documentation + +Note that we're aware the rendering of equations in markdown on GitHub for the documentation +is currently not great. +Eventually, the docs will transition to Read the Docs (or similar). + +## Overview + +The documentation currently includes + +- [Configuring solvers in the input file](input.md) + - [An extended overview of various conventions for optically thin cooling implementations](cooling_notes.md) + - [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) +- Detailed descriptions of more complex problem generators + - [Galaxy Cluster and Cluster-like Problem Setup](cluster.md) + - [Driven turbulence](turbulence.md) + +## Tutorial + +An AthenaPK tutorial was given as part of the **Towards exascale-ready astrophysics** +workshop https://indico3-jsc.fz-juelich.de/event/169/ taking place 25-27 Sep 2024 online. + +The material is currently located at https://github.com/pgrete/athenapk_tutorial + +While the instructions for building the code are specific to the workshop environment +the tutorial itself should translate directly to other environments/systems. + +## Parthenon documenation + +Many paramters/options are directly controlled through the Parthenon framework +(both with regard to building and in the input file). + +While the [Parthenon documenation](https://parthenon-hpc-lab.github.io/parthenon) is +more geared towards developers it also contains useful information for users. \ No newline at end of file diff --git a/docs/img/turb_acc.png b/docs/img/turb_acc.png new file mode 100644 index 00000000..b828a66e Binary files /dev/null and b/docs/img/turb_acc.png differ diff --git a/docs/img/turb_evol.png b/docs/img/turb_evol.png new file mode 100644 index 00000000..de7aa7b7 Binary files /dev/null and b/docs/img/turb_evol.png differ diff --git a/docs/img/turb_spec.png b/docs/img/turb_spec.png new file mode 100644 index 00000000..e6a8c565 Binary files /dev/null and b/docs/img/turb_spec.png differ diff --git a/docs/pgen.md b/docs/pgen.md index 0b58b9a2..9a2f58e8 100644 --- a/docs/pgen.md +++ b/docs/pgen.md @@ -1,4 +1,4 @@ -# Problem generators +# Custom problem generators Different simulation setups in AthenaPK are controlled via the so-called problem generators. New problem generators can easily be added and we are happy to accept and merge contibuted problem diff --git a/docs/turbulence.md b/docs/turbulence.md new file mode 100644 index 00000000..b51b4b70 --- /dev/null +++ b/docs/turbulence.md @@ -0,0 +1,107 @@ +# Driven turbulence simulations + +The turbulence problem generator uses explicit inverse Fourier transformations (iFTs) +on each meshblock in order to reduce communication during the iFT. +Thus, it is only efficient if comparatively few modes are used (say < 100). + +Quite generally, driven turbulence simulations start from uniform initial conditions +(uniform density and pressure, some initial magnetic field configuration in case of an +MHD setup, and the fluid at rest) and reach a state of stationary, isotropic (or anisotropic +depending on the strength of the background magnetic field) turbulence after one to few +large eddy turnover times (again depending on the background magnetic field strength). +The large eddy turnover time is usually defined as `T = L/U` with `L` being the scale +of the largest eddies and `U` the root mean square Mach number in the stationary regime. + +The current implementation uses the following forcing spectrum +`(k/k_peak)^2 * (2 - (k/k_peak)^2)`. +Here, `k_peak` is the peak wavenumber of the forcing spectrum. It is related the scales of the largest eddies as +`L = 1/k_f` given that a box size of 1 is currently assumed/hardcoded. + +## Problem setup + +An example parameter file can be found in `inputs/turbulence.in`. + +A typical setup contains the following blocks in the input file: + +``` + +problem_id = turbulence + + +rho0 = 1.0 # initial mean density +p0 = 1.0 # initial mean pressure +b0 = 0.01 # initial magnetic field strength +b_config = 0 # 0 - net flux; 1 - no net flux uniform B; 2 - non net flux sin B; 4 - field loop +kpeak = 2.0 # characteristic wavenumber +corr_time = 1.0 # autocorrelation time of the OU forcing process +rseed = 20190729 # random seed of the OU forcing process +sol_weight = 1.0 # solenoidal weight of the acceleration field +accel_rms = 0.5 # root mean square value of the acceleration field +num_modes = 30 # number of wavemodes + + +k_1_0 = +2 +k_1_1 = -1 +k_1_2 = +0 +k_2_0 = +1 +... +``` + +The following parameters can be changed to control both the initial state: + +- `rho0` initial mean density +- `p0` initial mean thermal pressure +- `b0` initial mean magnetic field strength +- `b_config` + - `0`: net flux case (uniform B_x) + - `1`: no net flux case (uniform B_x with changing sign in each half of the box) + - `2`: no net flux with initial sinosoidal B_x field + - `3`: deprecated + - `4`: closed field loop/cylinder in the box (in x-y plane) located at + - `x0=0.5` (default) + - `y0=0.5` (default) + - `z0=0.5` (default) + - and radius `loop_rad=0.25` + +as well as the driving field: + +- `kpeak` peak wavenumber of the forcing spectrum. Make sure to update the wavemodes to match `kpeak`, see below. +- `corr_time` autocorrelation time of the acceleration field (in code units). +Using delta-in-time forcing, i.e., a very low value, is discouraged, see [Grete et al. 2018 ApJL](https://iopscience.iop.org/article/10.3847/2041-8213/aac0f5). +- `rseed` random seed for the OU process. Only change for new simulation, but keep unchanged for restarting simulations. +- `sol_weight` solenoidal weight of the acceleration field. `1.0` is purely solenoidal/rotational and `0.0` is purely dilatational/compressive. Any value between `0.0` and `1.0` is possible. The parameter is related to the resulting rotational power in the 3D acceleration field as +`1. - ((1-sol_weight)^2/(1-2*sol_weight+3*sol_weight^2))`, see eq (9) in [Federrath et al. 2010 A&A]( +https://doi.org/10.1051/0004-6361/200912437). +- `accel_rms` root mean square value of the acceleration (controls the "strength" of the forcing field) +- `num_modes` number of wavemodes that are specified in the `` section of the parameter file. +The modes are specified manually as an explicit inverse FT is performed and only modes set are included (all others are assumed to be 0). +This is done to make the global inverse FT possible without any +expensive communication between blocks but this becomes excessively +expensiv for large number of modes. +Typically using a few tens of modes is a good choice in practice. +In order to generate a set of modes run the `inputs/generate_fmturb_modes.py` script and replace +the corresponding parts of the parameter file with the output of the script. +Within the script, the top three variables (`k_peak`, `k_high`, and `k_low`) need to be adjusted in +order to generate a complete set (i.e., all) of wavemodes. +Important, the `k_peak` in the script should match the `k_peak` set +in the input file. +Alternatively, wavemodes can be chosen/defined manually, e.g., if not all wavemodes are desired or +only individual modes should be forced. + +## Typical results + +The results shown here are obtained from running simulations with the parameters given in the next section. + +### High level temporal evolution +![image](img/turb_evol.png) + +### Power spectra +![image](img/turb_spec.png) + +### Consistency of acceleration field +As each meshblock does a full iFT of all modes the following slices from a run with 8 meshblocks +illustrate that there is no discontinuities at the meshblock boundary. + +Plot shows x-, y-, and z-acceleration (in rows top to bottom) slices in the x-, y-, and z-direction (in columns from left to right). + +![image](img/turb_acc.png) \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index fa0c4dec..f239b197 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -81,7 +81,7 @@ int main(int argc, char *argv[]) { pman.app_input->ProblemGenerator = field_loop::ProblemGenerator; Hydro::ProblemInitPackageData = field_loop::ProblemInitPackageData; } else if (problem == "kh") { - pman.app_input->ProblemGenerator = kh::ProblemGenerator; + pman.app_input->MeshProblemGenerator = kh::ProblemGenerator; } else if (problem == "rand_blast") { pman.app_input->ProblemGenerator = rand_blast::ProblemGenerator; Hydro::ProblemInitPackageData = rand_blast::ProblemInitPackageData; diff --git a/src/pgen/kh.cpp b/src/pgen/kh.cpp index 547f2a2e..68936935 100644 --- a/src/pgen/kh.cpp +++ b/src/pgen/kh.cpp @@ -32,61 +32,32 @@ // AthenaPK headers #include "../main.hpp" +#include "utils/error_checking.hpp" namespace kh { using namespace parthenon::driver::prelude; //---------------------------------------------------------------------------------------- -//! \fn void MeshBlock::ProblemGenerator(ParameterInput *pin) +//! \fn void Mesh::ProblemGenerator(Mesh *pm, ParameterInput *pin, MeshData *md) // \brief Problem Generator for the Kelvin-Helmholtz test -void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { +void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { auto vflow = pin->GetReal("problem/kh", "vflow"); auto iprob = pin->GetInteger("problem/kh", "iprob"); + // Get pointer to first block (always exists) for common data like loop bounds + auto pmb = md->GetBlockData(0)->GetBlockPointer(); auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); auto gam = pin->GetReal("hydro", "gamma"); auto gm1 = (gam - 1.0); - // initialize conserved variables - auto &rc = pmb->meshblock_data.Get(); - auto &u_dev = rc->Get("cons").data; - auto &coords = pmb->coords; - // initializing on host - auto u = u_dev.GetHostMirrorAndCopy(); + // Initialize conserved variables + // Get a MeshBlockPack on device with all conserved variables + const auto &cons = md->PackVariables(std::vector{"cons"}); + const auto num_blocks = md->NumBlocks(); - std::mt19937 gen(pmb->gid); // Standard mersenne_twister_engine seeded with gid - std::uniform_real_distribution ran(-0.5, 0.5); - - //--- iprob=1. Uniform stream with density ratio "drat" located in region -1/4GetReal("problem/kh", "drat"); - Real amp = pin->GetReal("problem/kh", "amp"); - for (int k = kb.s; k <= kb.e; k++) { - for (int j = jb.s; j <= jb.e; j++) { - for (int i = ib.s; i <= ib.e; i++) { - u(IDN, k, j, i) = 1.0; - u(IM1, k, j, i) = vflow + amp * ran(gen); - u(IM2, k, j, i) = amp * ran(gen); - u(IM3, k, j, i) = 0.0; - if (std::abs(coords.Xc<2>(j)) < 0.25) { - u(IDN, k, j, i) = drat; - u(IM1, k, j, i) = -drat * (vflow + amp * ran(gen)); - u(IM2, k, j, i) = drat * amp * ran(gen); - } - // Pressure scaled to give a sound speed of 1 with gamma=1.4 - u(IEN, k, j, i) = - 2.5 / gm1 + - 0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i); - } - } - } - } + //--- iprob=1. This was the classic, unresolved K-H test. //--- iprob=2. Uniform density medium moving at +/-vflow seperated by a single shear // layer with tanh() profile at y=0 with a single mode perturbation, reflecting BCs at @@ -97,9 +68,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real amp = pin->GetReal("problem/kh", "amp"); Real a = 0.02; Real sigma = 0.2; - for (int k = kb.s; k <= kb.e; k++) { - for (int j = jb.s; j <= jb.e; j++) { - for (int i = ib.s; i <= ib.e; i++) { + pmb->par_for( + "KHI: iprob2", 0, num_blocks - 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) { + const auto &u = cons(b); + const auto &coords = cons.GetCoords(b); u(IDN, k, j, i) = 1.0; u(IM1, k, j, i) = vflow * std::tanh((coords.Xc<2>(j)) / a); u(IM2, k, j, i) = amp * std::cos(2.0 * M_PI * coords.Xc<1>(i)) * @@ -108,23 +81,22 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { u(IEN, k, j, i) = 1.0 / gm1 + 0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i); - } - } - } - } + }); - //--- iprob=3. Test in SR paper (Beckwith & Stone, ApJS 193, 6, 2011). Gives two - // resolved shear layers with tanh() profiles for velocity and density located at - // y = +/- 0.5, density one in middle and 0.01 elsewhere, single mode perturbation. + } else if (iprob == 3) { + //--- iprob=3. Test in SR paper (Beckwith & Stone, ApJS 193, 6, 2011). Gives two + // resolved shear layers with tanh() profiles for velocity and density located at + // y = +/- 0.5, density one in middle and 0.01 elsewhere, single mode perturbation. - if (iprob == 3) { // Read/set problem parameters Real amp = pin->GetReal("problem/kh", "amp"); Real a = 0.01; Real sigma = 0.1; - for (int k = kb.s; k <= kb.e; k++) { - for (int j = jb.s; j <= jb.e; j++) { - for (int i = ib.s; i <= ib.e; i++) { + pmb->par_for( + "KHI: iprob3", 0, num_blocks - 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) { + const auto &u = cons(b); + const auto &coords = cons.GetCoords(b); u(IDN, k, j, i) = 0.505 + 0.495 * std::tanh((std::abs(coords.Xc<2>(j)) - 0.5) / a); u(IM1, k, j, i) = vflow * std::tanh((std::abs(coords.Xc<2>(j)) - 0.5) / a); @@ -139,17 +111,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { u(IEN, k, j, i) = 1.0 / gm1 + 0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i); - } - } - } - } + }); - //--- iprob=4. "Lecoanet" test, resolved shear layers with tanh() profiles for velocity - // and density located at z1=0.5, z2=1.5 two-mode perturbation for fully periodic BCs + } else if (iprob == 4) { + //--- iprob=4. "Lecoanet" test, resolved shear layers with tanh() profiles for + // velocity + // and density located at z1=0.5, z2=1.5 two-mode perturbation for fully periodic BCs - // To promote symmetry of FP errors about midplanes, rescale z' = z - 1. ; x' = x - 0.5 - // so that domain x1 = [-0.5, 0.5] and x2 = [-1.0, 1.0] is centered about origin - if (iprob == 4) { + // To promote symmetry of FP errors about midplanes, rescale z' = z - 1. ; x' = x - + // 0.5 so that domain x1 = [-0.5, 0.5] and x2 = [-1.0, 1.0] is centered about origin // Read/set problem parameters Real amp = pin->GetReal("problem/kh", "amp"); // unstratified problem is the default @@ -164,9 +134,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real z1 = -0.5; // z1' = z1 - 1.0 Real z2 = 0.5; // z2' = z2 - 1.0 - for (int k = kb.s; k <= kb.e; k++) { - for (int j = jb.s; j <= jb.e; j++) { - for (int i = ib.s; i <= ib.e; i++) { + pmb->par_for( + "KHI: iprob4", 0, num_blocks - 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) { + const auto &u = cons(b); + const auto &coords = cons.GetCoords(b); // Lecoanet (2015) equation 8a) Real dens = 1.0 + 0.5 * drho_rho0 * (std::tanh((coords.Xc<2>(j) - z1) / a) - @@ -233,40 +205,39 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { P0 / gm1 + 0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) / u(IDN, k, j, i); - } - } - } - // copy initialized vars to device - u_dev.DeepCopy(u); - } + }); - //--- iprob=5. Uniform stream with density ratio "drat" located in region -1/4GetReal("problem/kh", "a"); Real sigma = pin->GetReal("problem/kh", "sigma"); Real drat = pin->GetReal("problem/kh", "drat"); Real amp = pin->GetReal("problem/kh", "amp"); - for (int k = kb.s; k <= kb.e; k++) { - for (int j = jb.s; j <= jb.e; j++) { - for (int i = ib.s; i <= ib.e; i++) { + pmb->par_for( + "KHI: iprob5", 0, num_blocks - 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) { + const auto &u = cons(b); + const auto &coords = cons.GetCoords(b); Real w = (std::tanh((std::abs(coords.Xc<2>(j)) - 0.25) / a) + 1.0) * 0.5; u(IDN, k, j, i) = w + (1.0 - w) * drat; - u(IM1, k, j, i) = w * vflow - (1.0 - w) * vflow * drat; + u(IM1, k, j, i) = u(IDN, k, j, i) * vflow * (w - 0.5); u(IM2, k, j, i) = - u(IDN, k, j, i) * amp * std::sin(2.0 * 2.0 * M_PI * coords.Xc<1>(i)) * + u(IDN, k, j, i) * amp * std::cos(2.0 * 2.0 * M_PI * coords.Xc<1>(i)) * std::exp(-SQR(std::abs(coords.Xc<2>(j)) - 0.25) / (sigma * sigma)); u(IM3, k, j, i) = 0.0; // Pressure scaled to give a sound speed of 1 with gamma=1.4 u(IEN, k, j, i) = 2.5 / gm1 + - 0.25 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i); - } - } - } + 0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i); + }); + } else { + PARTHENON_FAIL("Unknow iprob for KHI pgen.") } } diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index ddea80f3..877a1f13 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -84,7 +84,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg namespace kh { using namespace parthenon::driver::prelude; -void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin); +void ProblemGenerator(Mesh *pm, parthenon::ParameterInput *pin, MeshData *md); } // namespace kh namespace rand_blast { using namespace parthenon::driver::prelude;