From 6b31964fa474471c671acf293d0d38a54f050249 Mon Sep 17 00:00:00 2001 From: Gurkirat Singh Date: Sun, 11 Jun 2023 19:32:00 +0530 Subject: [PATCH 1/6] update parthenon --- external/parthenon | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/parthenon b/external/parthenon index 7da5c6a2..ad23761c 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 7da5c6a29792b47cff2411390d91ef2c8f386304 +Subproject commit ad23761c10528b86a7c2ed914c6798ad14607110 From 2a09c654928f0f16b4884b17cd18cc4fb0d3ab65 Mon Sep 17 00:00:00 2001 From: Gurkirat Singh Date: Sun, 11 Jun 2023 19:37:43 +0530 Subject: [PATCH 2/6] implement RTI problem --- src/main.cpp | 2 + src/pgen/CMakeLists.txt | 1 + src/pgen/pgen.hpp | 5 ++ src/pgen/rt.cpp | 110 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 118 insertions(+) create mode 100644 src/pgen/rt.cpp diff --git a/src/main.cpp b/src/main.cpp index 1d9ef467..8ab6ed25 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -84,6 +84,8 @@ int main(int argc, char *argv[]) { pman.app_input->ProblemGenerator = rand_blast::ProblemGenerator; Hydro::ProblemInitPackageData = rand_blast::ProblemInitPackageData; Hydro::ProblemSourceFirstOrder = rand_blast::RandomBlasts; + } else if (problem == "rt") { + pman.app_input->ProblemGenerator = rt::ProblemGenerator; } else if (problem == "cluster") { pman.app_input->ProblemGenerator = cluster::ProblemGenerator; pman.app_input->MeshBlockUserWorkBeforeOutput = cluster::UserWorkBeforeOutput; diff --git a/src/pgen/CMakeLists.txt b/src/pgen/CMakeLists.txt index 443c4465..3a325196 100644 --- a/src/pgen/CMakeLists.txt +++ b/src/pgen/CMakeLists.txt @@ -19,6 +19,7 @@ target_sources(athenaPK PRIVATE linear_wave_mhd.cpp orszag_tang.cpp rand_blast.cpp + rt.cpp sod.cpp turbulence.cpp ) diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 84e36b00..3c28a873 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -121,4 +121,9 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin); void Cleanup(); } // namespace turbulence +namespace rt { + using namespace parthenon::driver::prelude; + void ProblemGenerator(MeshBlock *pm, parthenon::ParameterInput *pin); +} + #endif // PGEN_PGEN_HPP_ diff --git a/src/pgen/rt.cpp b/src/pgen/rt.cpp new file mode 100644 index 00000000..16a3433b --- /dev/null +++ b/src/pgen/rt.cpp @@ -0,0 +1,110 @@ +//! \file rt.cpp +//! \brief Rayleigh Taylor problem generator +//! Problem domain should be -.05 < x < .05; -.05 < y < .05, -.1 < z < .1, gamma=5/3 to +//! match Dimonte et al. Interface is at z=0; perturbation added to Vz. Gravity acts in +//! z-dirn. Special reflecting boundary conditions added in x3. A=1/2. Options: +//! - iprob = 1 -- Perturb V3 using single mode +//! - iprob = 2 -- Perturb V3 using multiple mode +//! - iprob = 3 -- B rotated by "angle" at interface, multimode perturbation +//! +// C headers + +// C++ headers +#include // min, max +#include // log +#include // strcmp() +#include + +// Parthenon headers +#include "mesh/mesh.hpp" +#include "parthenon/driver.hpp" +#include "parthenon/package.hpp" +#include "utils/error_checking.hpp" + +// AthenaPK headers +#include "../main.hpp" + +namespace rt { + using namespace parthenon::driver::prelude; + + void SetupSingleMode(MeshBlock *pmb, parthenon::ParameterInput *pin) { + if (pmb->pmy_mesh->ndim == 1) { + PARTHENON_THROW("This problem should be either in 2d or 3d."); + return; + } + + Real kx = 2.0*(M_PI)/(pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min); + Real ky = 2.0*(M_PI)/(pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min); + Real kz = 2.0*(M_PI)/(pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min); + + Real amp = pin->GetReal("problem/rt","amp"); + Real drat = pin->GetOrAddReal("problem/rt","drat",3.0); + Real grav_acc = pin->GetReal("hydro","const_accel_val"); + + 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); + Real p0 = 1.0/gam; + + // 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(); + + for (size_t k = kb.s; k < kb.e; k++) { + for (size_t j = jb.s; j < jb.e; j++) { + for (size_t i = ib.s; j < ib.e; i++) { + auto x1v = coords.Xc<1>(i); + auto x2v = coords.Xc<2>(j); + auto x3v = coords.Xc<3>(k); + + switch (pmb->pmy_mesh->ndim) { + case 2:{ + Real den=1.0; + if (x2v > 0.0) den *= drat; + + u(IM2,k,j,i) = (1.0 + cos(kx*x1v))*(1.0 + cos(ky*x2v))/4.0; + u(IDN,k,j,i) = den; + u(IM1,k,j,i) = 0.0; + u(IM2,k,j,i) *= (den*amp); + u(IM3,k,j,i) = 0.0; + u(IEN,k,j,i) = (p0 + grav_acc*den*x2v)/gm1 + 0.5*SQR(u(IM2,k,j,i))/den; + } + break; + case 3: { + Real den=1.0; + if (x3v > 0.0) den *= drat; + u(IM3,k,j,i) = (1.0+cos(kx*x1v))*(1.0+cos(ky*x2v))*(1.0+cos(kz*x3v))/8.0; + u(IDN,k,j,i) = den; + u(IM1,k,j,i) = 0.0; + u(IM2,k,j,i) = 0.0; + u(IM3,k,j,i) *= (den*amp); + u(IEN,k,j,i) = (p0 + grav_acc*den*x3v)/gm1 + 0.5*SQR(u(IM3,k,j,i))/den; + } + break; + } + } + } + } + } + + void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) { + auto iprob = pin->GetOrAddInteger("problem/rt", "iprob", 1); + + switch (iprob) { + case 1: + SetupSingleMode(pmb, pin); + break; + + default: + std::stringstream msg; + msg << "Problem type " << iprob << " is not supported."; + PARTHENON_THROW(msg.str()); + } + } +} // namespace rt \ No newline at end of file From 0254da3aba2973bd7347b2c650260a9d32aa2b33 Mon Sep 17 00:00:00 2001 From: Gurkirat Singh Date: Sun, 11 Jun 2023 19:38:42 +0530 Subject: [PATCH 3/6] input file for rti --- inputs/rt.in | 69 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 inputs/rt.in diff --git a/inputs/rt.in b/inputs/rt.in new file mode 100644 index 00000000..b447870c --- /dev/null +++ b/inputs/rt.in @@ -0,0 +1,69 @@ + +problem = Rayleigh-Taylor instability + + +problem_id = rt + + +iprob = 1 +amp = 0.01 +drat = 2.0 + + +refinement = adaptive +numlevel = 3 +nghost = 4 + +nx1 = 100 # Number of zones in X1-direction +x1min = -0.16666667 # minimum value of X1 +x1max = 0.16666667 # maximum value of X1 +ix1_bc = periodic # inner-X1 boundary flag +ox1_bc = periodic # outer-X1 boundary flag + +nx2 = 200 # Number of zones in X2-direction +x2min = -0.5 # minimum value of X2 +x2max = 0.5 # maximum value of X2 +ix2_bc = reflecting # inner-X2 boundary flag +ox2_bc = reflecting # outer-X2 boundary flag + +nx3 = 1 # Number of zones in X3-direction +x3min = -0.5 # minimum value of X3 +x3max = 0.5 # maximum value of X3 +ix3_bc = periodic # inner-X3 boundary flag +ox3_bc = periodic # outer-X3 boundary flag + + + +nx1 = 100 # Number of cells in each MeshBlock, X1-dir +nx2 = 200 # Number of cells in each MeshBlock, X2-dir +nx3 = 1 # Number of cells in each MeshBlock, X3-dir + + +integrator = vl2 +cfl = 0.4 +tlim = 10.0 +nlim = 100000 +perf_cycle_offset = 2 # number of inital cycles not to be included in perf calc +ncycle_out_mesh = -1000 + + +fluid = euler +eos = adiabatic +riemann = hlle +reconstruction = ppm +gamma = 1.666666666666667 # gamma = C_p/C_v +const_accel = true # add constant accleration source term +const_accel_val = -0.1 # value of constant acceleration +const_accel_dir = 2 # direction of constant acceleration +first_order_flux_correct = true +dfloor = 0.00000001 # unused, in [code units] +pfloor = 0.00000001 # unused, in [code units] + + +file_type = hdf5 +dt = 1.0 +variables = prim + + +file_type = hst +dt = 0.1 \ No newline at end of file From f5e2149024f839eba32b07103b8a7f7b8a158386 Mon Sep 17 00:00:00 2001 From: Gurkirat Singh Date: Mon, 12 Jun 2023 16:01:51 +0530 Subject: [PATCH 4/6] Update src/pgen/pgen.hpp Co-authored-by: Philipp Grete --- src/pgen/pgen.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 3c28a873..b309d239 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -123,7 +123,7 @@ void Cleanup(); namespace rt { using namespace parthenon::driver::prelude; - void ProblemGenerator(MeshBlock *pm, parthenon::ParameterInput *pin); + void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin); } #endif // PGEN_PGEN_HPP_ From 0ba1773a8bc7b7ef185a36715d3fba60befaf53e Mon Sep 17 00:00:00 2001 From: Gurkirat Singh Date: Mon, 12 Jun 2023 16:07:44 +0530 Subject: [PATCH 5/6] deep copy host vars to device --- src/pgen/rt.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pgen/rt.cpp b/src/pgen/rt.cpp index 16a3433b..7739cd11 100644 --- a/src/pgen/rt.cpp +++ b/src/pgen/rt.cpp @@ -91,6 +91,7 @@ namespace rt { } } } + u_dev.DeepCopy(u); } void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) { From 715779391b21e1927e627843d7fa6f65ee5aeb4d Mon Sep 17 00:00:00 2001 From: Gurkirat Singh Date: Mon, 12 Jun 2023 16:27:08 +0530 Subject: [PATCH 6/6] use outflow bcs for x2, for now --- inputs/rt.in | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inputs/rt.in b/inputs/rt.in index b447870c..97a137c2 100644 --- a/inputs/rt.in +++ b/inputs/rt.in @@ -23,8 +23,8 @@ ox1_bc = periodic # outer-X1 boundary flag nx2 = 200 # Number of zones in X2-direction x2min = -0.5 # minimum value of X2 x2max = 0.5 # maximum value of X2 -ix2_bc = reflecting # inner-X2 boundary flag -ox2_bc = reflecting # outer-X2 boundary flag +ix2_bc = outflow # inner-X2 boundary flag +ox2_bc = outflow # outer-X2 boundary flag nx3 = 1 # Number of zones in X3-direction x3min = -0.5 # minimum value of X3