From 00b471169fced75ee3a31c98b47a363da21b23d9 Mon Sep 17 00:00:00 2001 From: HPC user Date: Mon, 25 Sep 2023 15:08:30 +0200 Subject: [PATCH] Updating spherical collapse routines --- src/pgen/cluster.cpp | 62 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 56 insertions(+), 6 deletions(-) diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 2882c96d..e2b4e7ca 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -691,12 +691,12 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { /************************************************************ * Initial parameters ************************************************************/ - + auto pmb = md->GetBlockData(0)->GetBlockPointer(); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); - + auto hydro_pkg = pmb->packages.Get("Hydro"); const auto fluid = hydro_pkg->Param("fluid"); auto const &cons = md->PackVariables(std::vector{"cons"}); @@ -712,6 +712,56 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + Real passive_scalar = 0.0; // Not useful here + + // Spherical collapse test with an initial overdensity + + if (spherical_collapse == true) { + + pmb->par_reduce( + "Init density field", 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, Real &lsum) { + + auto pmbb = md->GetBlockData(b)->GetBlockPointer(); // Meshblock b + + const auto gis = pmbb->loc.lx1 * pmb->block_size.nx1; + const auto gjs = pmbb->loc.lx2 * pmb->block_size.nx2; + const auto gks = pmbb->loc.lx3 * pmb->block_size.nx3; + + const auto &coords = cons.GetCoords(b); + const auto &u = cons(b); + + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); // Computing radius + const Real overdensity_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overdensity_radius", 0.01); + const Real background_density = pin->GetOrAddReal("problem/cluster/init_perturb", "background_density", 3); + const Real foreground_density = pin->GetOrAddReal("problem/cluster/init_perturb", "foreground_density", 150); + + // Setting an initial + u(IDN, k, j, i) = background_density; + + // For any cell at a radius less then overdensity radius, set the density to foreground_density + if (r < overdensity_radius){ + u(IDN, k, j, i) = foreground_density; + } + + + }, + passive_scalar); + + } + + /* -------------- Setting up a clumpy atmosphere -------------- + + 1) Extract the values of the density from an input hdf5 file using H5Easy + 2) Initiate the associated density field + 3) Optionnaly, add some overpressure ring to check behavior (overpressure_ring bool) + 4) Optionnaly, add a central overdensity + + */ + if (init_perturb_rho == true) { auto init_perturb_rho_file = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_file"); @@ -760,9 +810,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const Real z = coords.Xc<3>(k); const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); // Computing radius const Real width = 0.002; - const Real overpressure_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overpressure_radius", 0.01); + const Real overdensity_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overdensity_radius", 0.01); - if (r > 0 && r < overpressure_radius){ + if (r > 0 && r < overdensity_radius){ //u(IEN, k, j, i) *= 100; @@ -788,11 +838,11 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const Real y = coords.Xc<2>(j); const Real z = coords.Xc<3>(k); const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); // Computing radius - const Real overpressure_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overpressure_radius", 0.01); + const Real overdensity_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overdensity_radius", 0.01); u(IDN, k, j, i) = 3; - if (r > 0 && r < overpressure_radius){ + if (r > 0 && r < overdensity_radius){ u(IDN, k, j, i) *= 50;