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 diff --git a/inputs/rt.in b/inputs/rt.in new file mode 100644 index 00000000..97a137c2 --- /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 = 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 +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 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..b309d239 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 *pmb, 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..7739cd11 --- /dev/null +++ b/src/pgen/rt.cpp @@ -0,0 +1,111 @@ +//! \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; + } + } + } + } + u_dev.DeepCopy(u); + } + + 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