From 413935b10c49fc6c7a2a080804d12ebecc309fd0 Mon Sep 17 00:00:00 2001 From: jtlaune Date: Tue, 10 Dec 2024 09:46:14 -0700 Subject: [PATCH 1/6] Add drag files and nu disk pgen --- src/drag/drag.cpp | 8 + src/drag/drag.hpp | 54 +++- src/pgen/nu_disk.hpp | 740 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 787 insertions(+), 15 deletions(-) create mode 100644 src/pgen/nu_disk.hpp diff --git a/src/drag/drag.cpp b/src/drag/drag.cpp index 39f8a57..e3d7840 100644 --- a/src/drag/drag.cpp +++ b/src/drag/drag.cpp @@ -42,6 +42,14 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add("x2max", pin->GetReal("parthenon/mesh", "x2max")); params.Add("x3max", pin->GetReal("parthenon/mesh", "x3max")); + params.Add("dslope",pin->GetReal("problem", "dslope")); + params.Add("tslope",pin->GetReal("problem", "tslope")); + params.Add("h0",pin->GetReal("problem", "h0")); + params.Add("r0",pin->GetReal("problem", "r0")); + params.Add("gm",pin->GetReal("gravity", "gm")); + params.Add("omf",pin->GetReal("rotating_frame", "omega")); + params.Add("nu",pin->GetReal("gas/viscosity", "nu")); + const bool do_gas = pin->GetOrAddBoolean("physics", "gas", true); const bool do_dust = pin->GetOrAddBoolean("physics", "dust", false); diff --git a/src/drag/drag.hpp b/src/drag/drag.hpp index 3f9adf8..e878305 100644 --- a/src/drag/drag.hpp +++ b/src/drag/drag.hpp @@ -14,6 +14,7 @@ #define DRAG_DRAG_HPP_ // Parthenon includes +#include #include // Artemis includes @@ -86,25 +87,24 @@ struct SelfDragParams { } SelfDragParams(std::string block_name, ParameterInput *pin) { ix[0] = pin->GetOrAddReal(block_name, "inner_x1", -Big()); - ix[1] = pin->GetOrAddReal(block_name, "inner_x2", -Big()); - ix[2] = pin->GetOrAddReal(block_name, "inner_x3", -Big()); + ix[1] = pin->GetOrAddReal(block_name, "pos_z_sh", Big()); irate[0] = pin->GetOrAddReal(block_name, "inner_x1_rate", 0.0); - irate[1] = pin->GetOrAddReal(block_name, "inner_x2_rate", 0.0); - irate[2] = pin->GetOrAddReal(block_name, "inner_x3_rate", 0.0); + irate[1] = pin->GetOrAddReal(block_name, "z_rate", 0.0); ox[0] = pin->GetOrAddReal(block_name, "outer_x1", Big()); - ox[1] = pin->GetOrAddReal(block_name, "outer_x2", Big()); - ox[2] = pin->GetOrAddReal(block_name, "outer_x3", Big()); + ox[1] = pin->GetOrAddReal(block_name, "neg_z_sh", Big()); orate[0] = pin->GetOrAddReal(block_name, "outer_x1_rate", 0.0); - orate[1] = pin->GetOrAddReal(block_name, "outer_x2_rate", 0.0); - orate[2] = pin->GetOrAddReal(block_name, "outer_x3_rate", 0.0); + orate[1] = pin->GetOrAddReal(block_name, "z_rate", 0.0); + damp_to_visc = pin->GetOrAddBoolean(block_name, "damp_to_visc", false); for (int i = 0; i < 3; i++) { PARTHENON_REQUIRE(irate[i] >= 0.0, "The damping rate in the x1 direction must be >= 0"); + if (i != 2) { PARTHENON_REQUIRE(ix[i] <= ox[i], "The damping bounds must have inner_x1 <= outer_x1"); + } } } }; @@ -185,6 +185,15 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt const int multi_d = (ndim >= 2); const int three_d = (ndim == 3); + const Real p = drag_pkg->template Param("dslope"); + const Real q = drag_pkg->template Param("tslope"); + const Real h0 = drag_pkg->template Param("h0"); + const Real r0 = drag_pkg->template Param("r0"); + const Real gm = drag_pkg->template Param("gm"); + const Real omf = drag_pkg->template Param("omf"); + const Real nu = drag_pkg->template Param("nu"); + const Real flare = 0.5 * (1.0 + q); + static auto desc = MakePackDescriptor *md, const Real time, const Real dt const auto &hx = coords.GetScaleFactors(); const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + const Real PI_2 = 1.5707963267948966; + // Compute the ramp for this cell // Ramps are quadratic, eg. the left regions is SQR( (X - ix)/(ix - xmin) ) if (do_gas) { + const Real H = xcyl[0] * h0 * std::pow(xcyl[0] / r0, flare); const Real fx1 = dt * (gasp.irate[0] * ((xv[0] < gasp.ix[0]) * SQR((xv[0] - gasp.ix[0]) / (gasp.ix[0] - x1min))) + @@ -224,6 +236,7 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt SQR((xv[2] - gasp.ix[2]) / (gasp.ix[2] - x3min))) + gasp.orate[2] * ((xv[2] > gasp.ox[2]) * SQR((xv[2] - gasp.ox[2]) / (gasp.ox[2] - x3max)))); + for (int n = 0; n < vmesh.GetSize(b, gas::cons::density()); ++n) { const Real &dens = vmesh(b, gas::cons::density(n), k, j, i); const Real vg[3] = { @@ -234,23 +247,34 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt const Real sieg = ArtemisUtils::GetSpecificInternalEnergy( vmesh, b, n, k, j, i, de_switch, dflr_gas, sieflr_gas, hx); - Real vd[3] = {0., 0., 0.}; + const Real OmKmid = std::sqrt(gm / (xcyl[0] * xcyl[0] * xcyl[0])); + const Real Omg = OmKmid * (1 + 0.5 * SQR(H / xcyl[0]) * + (p + q + 0.5 * q * SQR(xcyl[2] / H))); + const Real vp = Omg * xcyl[0]; + const Real vR = -nu * + (6 * p - 2 * q + 3 + (5 * q + 9) * SQR(xcyl[2] / H)) / + (2 * xcyl[0]); - Diffusion::DiffusionCoeff dcoeff; - const Real mu = dcoeff.Get(dp, coords, dens, sieg, eos_d); - const Real vR = -1.5 * mu / (xcyl[0] * dens); - vd[0] = ex1[0] * vR; - vd[1] = ex2[0] * vR; - vd[2] = ex3[0] * vR; + const Real vz = (-p)*xcyl[2]/xcyl[0]*vR; + + const Real vcyl[3] = {vR, vp - omf * xcyl[0], vz}; + + const Real vd[3] = { + ArtemisUtils::VDot(vcyl, ex1), + ArtemisUtils::VDot(vcyl, ex2), + ArtemisUtils::VDot(vcyl, ex3) + }; // Ep - E = 0.5 d ( vp^2 - v^2 ) // (vp-v) . (vp + v) = dv . (2v + dv) = 2 dv.v + dv.dv const Real dm1 = -fx1 * dens * (vg[0] - vd[0]) / (1.0 + fx1); const Real dm2 = -fx2 * dens * (vg[1] - vd[1]) / (1.0 + fx2); const Real dm3 = -fx3 * dens * (vg[2] - vd[2]) / (1.0 + fx3); + vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) += hx[0] * dm1; vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) += hx[1] * dm2; vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) += hx[2] * dm3; + vmesh(b, gas::cons::total_energy(n), k, j, i) += dm1 * (vg[0] + 0.5 * dm1 / dens) + dm2 * (vg[1] + 0.5 * dm2 / dens) + dm3 * (vg[2] + 0.5 * dm3 / dens); diff --git a/src/pgen/nu_disk.hpp b/src/pgen/nu_disk.hpp new file mode 100644 index 0000000..22b7dfc --- /dev/null +++ b/src/pgen/nu_disk.hpp @@ -0,0 +1,740 @@ +//======================================================================================== +// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== +// AthenaXXX astrophysical plasma code +// Copyright(C) 2020 James M. Stone and the Athena code team +// Licensed under the 3-clause BSD License (the "LICENSE") +//======================================================================================== +#ifndef PGEN_NU_DISK_HPP_ +#define PGEN_NU_DISK_HPP_ +//! \file nu_disk.hpp +//! \brief Initializes a stratified Keplerian accretion disk. Initial conditions are in +//! vertical hydrostatic equilibrium. + +// NOTE(PDM): The following is adapted from the open-source Athena++ disk.cpp +// problem generator, and adapted for Parthenon/Artemis by PDM on 10/20/23, +// and adapted for the nu disk setup by JTL on 12/9/24. + +// C/C++ headers +#include +#include +#include +#include +#include +#include +#include + +// Artemis headers +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "gravity/nbody_gravity.hpp" +#include "nbody/nbody.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/eos/eos.hpp" + +using ArtemisUtils::EOS; +using ArtemisUtils::VI; + +namespace nu_disk { +//---------------------------------------------------------------------------------------- +//! \struct DiskParams +//! \brief container for disk parameters +struct DiskParams { + Real r0, h0; + Real p, q, flare; + Real rho0, dens_min, pres_min; + Real gm, Omega0, l0; + Real omf; + Real dust_to_gas; + Real rexp; + Real rcav; + Real Gamma, gamma_gas; + Real alpha, nu0, nu_indx; + Real mdot; + Real temp_soft2; + bool const_nu; + bool do_dust; + bool nbody_temp; + bool quiet_start; +}; + +//---------------------------------------------------------------------------------------- +//! \fn Real DenProfile +//! \brief Computes density profile at cylindrical R and z +KOKKOS_INLINE_FUNCTION +Real DenProfile(struct DiskParams pgen, const Real R, const Real z) { + const Real H = R * pgen.h0 * std::pow(R / pgen.r0, pgen.flare); + return std::max(pgen.dens_min, pgen.rho0*std::pow(R/pgen.r0, pgen.p)*std::exp(-0.5*SQR(z/H))); +} + +//---------------------------------------------------------------------------------------- +//! \fn Real DenProfile +//! \brief Computes temperature profile at cylindrical R and z +KOKKOS_INLINE_FUNCTION +Real TempProfile(struct DiskParams pgen, const Real R, const Real z) { + // P = K rho^Gamma + // T = T0 (rho/rho0)^(Gamma-1) + const Real rho = DenProfile(pgen, R, z); + const Real rho0 = DenProfile(pgen, R, 0.0); + const Real H = R * pgen.h0 * std::pow(R / pgen.r0, pgen.flare); + const Real ir1 = 1.0 / std::sqrt(R * R + pgen.temp_soft2); + const Real omk2 = SQR(pgen.Omega0) * ir1 * ir1 * ir1; + const Real T0 = omk2 * H * H / pgen.Gamma; + return T0 * std::pow(rho / rho0, pgen.Gamma - 1.0); +} + +//---------------------------------------------------------------------------------------- +//! \fn Real PresProfile +//! \brief Computes pressure profile at cylindrical R and z (via dens and temp profiles) +KOKKOS_INLINE_FUNCTION +Real PresProfile(struct DiskParams pgen, EOS eos, const Real tf, const Real R, + const Real z) { + const Real df = DenProfile(pgen, R, z); + return std::max(pgen.pres_min, eos.PressureFromDensityTemperature(df, tf)); +} + +//---------------------------------------------------------------------------------------- +//! \fn Real ViscosityProfile +//! \brief Computes viscosity profile at cylindrical R and z (via dens and temp profiles) +KOKKOS_INLINE_FUNCTION +Real ViscosityProfile(struct DiskParams pgen, EOS eos, const Real R, const Real z) { + return pgen.nu0 * std::pow(R / pgen.r0, pgen.nu_indx); +} + +//---------------------------------------------------------------------------------------- +//! \fn void ComputeDiskProfile +//! \brief Initialize vertical hydrostatic and radial centrifugal equilibrium disk profile +//! at a specified index/coordinate +template +KOKKOS_INLINE_FUNCTION void +ComputeDiskProfile(const struct DiskParams pgen, const parthenon::Coordinates_t &pco, + const int k, const int j, const int i, EOS eos_d, Real &gdens, + Real >emp, Real &gvel1, Real &gvel2, Real &gvel3, const bool do_dust, + Real &ddens, Real &dvel1, Real &dvel2, Real &dvel3, + ParArray1D particles, const int npart) { + // Extract coordinates + geometry::Coords coords(pco, k, j, i); + const auto &xv = coords.GetCellCenter(); + + const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + + // compute Keplerian solution + gdens = DenProfile(pgen, xcyl[0], xcyl[2]); + const Real rt = + pgen.nbody_temp + ? -pgen.gm / Gravity::NBodyPotential(coords, xv, particles, npart) + : xcyl[0]; + gtemp = TempProfile(pgen, rt, xcyl[2]); + + // Set v_phi to centrifugal equilibrium + // vp^2/R = grad(p) + vk^2/R + if (!pgen.const_nu) { + PARTHENON_FAIL("This version only works with nu viscosity."); + } + + const Real H = xcyl[0] * pgen.h0 * std::pow(xcyl[0] / pgen.r0, pgen.flare); + const Real OmKmid = std::sqrt(pgen.gm / (xcyl[0] * xcyl[0] * xcyl[0])); + const Real Omg = OmKmid * (1 + 0.5 * SQR(H / xcyl[0]) * + (pgen.p + pgen.q + 0.5 * pgen.q * SQR(xcyl[2] / H))); + const Real vp = Omg * xcyl[0]; + const Real vr = - pgen.nu0*(6*pgen.p-2*pgen.q+3+(5*pgen.q+9)*SQR(xcyl[2]/H))/(2*xcyl[0]); + const Real vz = (-pgen.p)*xcyl[2]/xcyl[0]*vr; + + // Construct the total cylindrical velocity + const Real vcyl[3] = {vr, vp - pgen.omf * xcyl[0], vz}; + + // and convert it to the problem geometry + gvel1 = ArtemisUtils::VDot(vcyl, ex1); + gvel2 = ArtemisUtils::VDot(vcyl, ex2); + gvel3 = ArtemisUtils::VDot(vcyl, ex3); + +//---------------------------------------------------------------------------------------- +//! \fn void InitDiskParams +//! \brief Extracts disk parameters from ParameterInput. +//! NOTE(PDM): In order for our user-defined BCs to be compatible with restarts, we must +//! reset the DiskParams struct upon initialization. +inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) { + auto &artemis_pkg = pmb->packages.Get("artemis"); + Params ¶ms = artemis_pkg->AllParams(); + if (!(params.hasKey("disk_params"))) { + DiskParams disk_params; + auto &grav_pkg = pmb->packages.Get("gravity"); + auto &gas_pkg = pmb->packages.Get("gas"); + + disk_params.gm = grav_pkg->Param("gm"); + disk_params.r0 = pin->GetOrAddReal("problem", "r0", 1.0); + disk_params.Omega0 = + std::sqrt(disk_params.gm / (disk_params.r0 * disk_params.r0 * disk_params.r0)); + disk_params.rho0 = pin->GetOrAddReal("problem", "rho0", 1.0); + disk_params.p = pin->GetOrAddReal("problem", "dslope", -2.25); + disk_params.h0 = pin->GetOrAddReal("problem", "h0", 0.05); + disk_params.gamma_gas = gas_pkg->Param("adiabatic_index"); + disk_params.Gamma = + pin->GetOrAddReal("problem", "polytropic_index", disk_params.gamma_gas); + + PARTHENON_REQUIRE(disk_params.Gamma >= 1, "problem/gamma needs to be >= 1"); + + disk_params.dens_min = pin->GetOrAddReal("problem", "dens_min", 1.0e-5); + disk_params.pres_min = pin->GetOrAddReal("problem", "pres_min", 1.0e-8); + disk_params.rexp = pin->GetOrAddReal("problem", "rexp", 0.0); + disk_params.rcav = pin->GetOrAddReal("problem", "rcav", 0.0); + disk_params.l0 = pin->GetOrAddReal("problem", "l0", 0.0); + disk_params.dust_to_gas = pin->GetOrAddReal("problem", "dust_to_gas", 0.01); + disk_params.temp_soft2 = pin->GetOrAddReal("problem", "temp_soft", 0.0); + + disk_params.do_dust = params.Get("do_dust"); + + Real q = pin->GetOrAddReal("problem", "tslope", -Big()); + Real flare = pin->GetOrAddReal("problem", "flare", -Big()); + + PARTHENON_REQUIRE((flare != -Big()) || (q != -Big()), + "Set flare or tslope in "); + + if (flare == -Big()) { + flare = 0.5 * (1.0 + q); + } else if (q == -Big()) { + q = 2.0 * flare - 1.; + } else { + PARTHENON_FAIL("Set either flare or tslope in not both!"); + } + disk_params.flare = flare; + disk_params.q = q; + disk_params.alpha = 0.0; + disk_params.nu0 = 0.0; + disk_params.nu_indx = 0.0; + disk_params.mdot = 0.0; + disk_params.quiet_start = pin->GetOrAddBoolean("problem", "quiet_start", false); + + if (params.Get("do_rotating_frame")) { + auto &rf_pkg = pmb->packages.Get("rotating_frame"); + disk_params.omf = rf_pkg->Param("omega"); + } else { + disk_params.omf = 0.0; + } + if (params.Get("do_viscosity")) { + const auto vtype = pin->GetString("gas/viscosity", "type"); + if (vtype == "alpha") { + PARTHENON_FAIL("Disk pgen is only compatible with constant viscosity"); + } else if (vtype == "constant") { + disk_params.const_nu = true; + disk_params.nu0 = pin->GetReal("gas/viscosity", "nu"); + disk_params.nu_indx = 0.0; + } else { + PARTHENON_FAIL("Disk pgen is only compatible with constant viscosity"); + } + if (pin->DoesParameterExist("problem", "mdot")) { + disk_params.mdot = pin->GetReal("problem", "mdot"); + disk_params.rho0 = disk_params.mdot / (3.0 * M_PI * disk_params.nu0); + } else { + disk_params.mdot = 3.0 * M_PI * disk_params.nu0 * disk_params.rho0; + } + } + disk_params.nbody_temp = pin->GetOrAddBoolean("problem", "nbody_temp", false); + disk_params.nbody_temp = disk_params.nbody_temp && params.Get("do_nbody"); + params.Add("disk_params", disk_params); + } +} + +//---------------------------------------------------------------------------------------- +//! \fn void DiskICImpl +//! \brief Set the state vectors of cell to the initial conditions +template +KOKKOS_INLINE_FUNCTION void +DiskICImpl(V1 v, const int b, const int k, const int j, const int i, V2 pco, EOS eos_d, + DiskParams dp, ParArray1D particles, const int npart) { + // gas + Real gdens = Null(), gtemp = Null(); + Real gvel1 = Null(), gvel2 = Null(), gvel3 = Null(); + // dust + Real ddens = Null(); + Real dvel1 = Null(), dvel2 = Null(), dvel3 = Null(); + ComputeDiskProfile(dp, pco, k, j, i, eos_d, gdens, gtemp, gvel1, gvel2, gvel3, + dp.do_dust, ddens, dvel1, dvel2, dvel3, particles, npart); + + // Set state vector + v(b, gas::prim::density(0), k, j, i) = gdens; + v(b, gas::prim::velocity(0), k, j, i) = gvel1; + v(b, gas::prim::velocity(1), k, j, i) = gvel2; + v(b, gas::prim::velocity(2), k, j, i) = gvel3; + v(b, gas::prim::sie(0), k, j, i) = + eos_d.InternalEnergyFromDensityTemperature(gdens, gtemp); + if (dp.do_dust) { + for (int n = 0; n < v.GetSize(b, dust::prim::density()); ++n) { + v(b, dust::prim::density(n), k, j, i) = ddens; + v(b, dust::prim::velocity(VI(n, 0)), k, j, i) = dvel1; + v(b, dust::prim::velocity(VI(n, 1)), k, j, i) = dvel2; + v(b, dust::prim::velocity(VI(n, 2)), k, j, i) = dvel3; + } + } +} + +//---------------------------------------------------------------------------------------- +//! \fn void ProblemGenerator::Disk() +//! \brief Sets initial conditions for disk problem +template +inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { + using parthenon::MakePackDescriptor; + + // Extract parameters from packages + auto artemis_pkg = pmb->packages.Get("artemis"); + const bool do_dust = artemis_pkg->Param("do_dust"); + + auto &gas_pkg = pmb->packages.Get("gas"); + auto eos_d = gas_pkg->template Param("eos_d"); + + // Disk parameters + auto disk_params = artemis_pkg->Param("disk_params"); + + // packing and capture variables for kernel + auto &md = pmb->meshblock_data.Get(); + for (auto &var : md->GetVariableVector()) { + if (!var->IsAllocated()) pmb->AllocateSparse(var->label()); + } + static auto desc = + MakePackDescriptor( + (pmb->resolved_packages).get()); + auto v = desc.GetPack(md.get()); + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + auto &pco = pmb->coords; + auto &dp = disk_params; + + ParArray1D particles; + int npart = 0; + if (dp.nbody_temp) { + auto &nbody_pkg = pmb->packages.Get("nbody"); + particles = nbody_pkg->template Param>("particles"); + npart = static_cast(particles.size()); + } + + pmb->par_for( + "disk", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + DiskICImpl(v, 0, k, j, i, pco, eos_d, dp, particles, npart); + }); +} + +//---------------------------------------------------------------------------------------- +//! \fn void Disk::DiskBoundaryVisc() +//! \brief Sets inner or outer X1 boundary condition to the initial condition +template +void DiskBoundaryVisc(std::shared_ptr> &mbd, bool coarse) { + + PARTHENON_REQUIRE(GEOM == Coordinates::cylindrical || + GEOM == Coordinates::spherical3D || + geometry::is_axisymmetric(), + "Viscous boundary conditions only work with spherical/cylindrical " + "radial boundaries"); + auto pmb = mbd->GetBlockPointer(); + + auto artemis_pkg = pmb->packages.Get("artemis"); + auto disk_params = artemis_pkg->Param("disk_params"); + + const bool do_dust = artemis_pkg->Param("do_dust"); + + auto &gas_pkg = pmb->packages.Get("gas"); + auto eos_d = gas_pkg->template Param("eos_d"); + + const bool fine = false; + + static auto descriptors = + ArtemisUtils::GetBoundaryPackDescriptorMap(mbd); + + auto v = descriptors[coarse].GetPack(mbd.get()); + + const auto &pco = (coarse) ? pmb->pmr->GetCoarseCoords() : pmb->coords; + auto &dp = disk_params; + const auto nb = IndexRange{0, 0}; + + // Boundary index arithmetic + int is = Null(); + int ie = Null(); + const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds; + if constexpr (BDY == IndexDomain::inner_x1) { + is = bounds.GetBoundsI(IndexDomain::interior, TE::CC).s; + } else if constexpr (BDY == IndexDomain::outer_x1) { + ie = bounds.GetBoundsI(IndexDomain::interior, TE::CC).e; + } else { + PARTHENON_FAIL( + "Viscous boundary conditions only work for the inner or outer radial boundary"); + } + const int ix1 = 0; + const int ix2 = 1; + const int ix3 = 2; + + pmb->par_for_bndry( + "DiskVisc", nb, BDY, parthenon::TopologicalElement::CC, coarse, fine, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + // We are extrapolating into the ghost zone. Extrapolation is done on cylinders. + // dP/dz = - rho grad(Phi) + // rho vp^2/R = rho vk^2/R + dP/dR + // We estimate the pressure gradients as grad(P).eR and grad(P).ez + // Vertical hydrostatic balance sets rho + // Radial centrifugal balance sets vp + const int ia[3] = {k, j, (BDY == IndexDomain::inner_x1) ? is : ie}; + const int ip1[3] = {k, j, (BDY == IndexDomain::inner_x1) ? is + 1 : ie}; + const int im1[3] = {k, j, (BDY == IndexDomain::inner_x1) ? is : ie - 1}; + + // Extract coordinates at k, j, i + geometry::Coords coords(pco, k, j, i); + const auto &xv = coords.GetCellCenter(); + const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + + // Extract coordinates at ia, im, ic + geometry::Coords ca(pco, ia[0], ia[1], ia[2]); + geometry::Coords cp1(pco, ip1[0], ip1[1], ip1[2]); + geometry::Coords cm1(pco, im1[0], im1[1], im1[2]); + const auto &xva = ca.GetCellCenter(); + const auto &[xcyla, scr1, scr2, scr3] = ca.ConvertToCylWithVec(xva); + const Real eRa[3] = {scr1[0], scr2[0], scr3[0]}; + const Real epa[3] = {scr1[1], scr2[1], scr3[1]}; + const Real eza[3] = {scr1[2], scr2[2], scr3[2]}; + const auto &xvp1 = cp1.GetCellCenter(); + const auto &[xcylp1, scr1p1, scr2p1, scr3p1] = cp1.ConvertToCylWithVec(xvp1); + const Real epp1[3] = {scr1p1[1], scr2p1[1], scr3p1[1]}; + const auto &xvm1 = cm1.GetCellCenter(); + const auto &[xcylm1, scr1m1, scr2m1, scr3m1] = cm1.ConvertToCylWithVec(xvm1); + const Real epm1[3] = {scr1m1[1], scr2m1[1], scr3m1[1]}; + + // Compute cell separations (using logarithmics if necessary) + const Real xma = std::log(xv[ix1] / xva[ix1]); + const Real dx = std::log(xvp1[ix1] / xvm1[ix1]); + const Real xmadx = xma / dx; + + const Real nua = ViscosityProfile(dp, eos_d, xcyla[0], xcyla[2]); + const Real nug = ViscosityProfile(dp, eos_d, xcyl[0], xcyl[2]); + + // Extrapolate gas density and specific internal energy + Real dgsie = std::log(v(0, gas::prim::sie(0), ip1[0], ip1[1], ip1[2]) / + v(0, gas::prim::sie(0), im1[0], im1[1], im1[2])); + const Real gsieexp = std::exp(dgsie * xmadx); + const Real sieg = v(0, gas::prim::sie(0), ia[0], ia[1], ia[2]) * gsieexp; + const Real rhoa = v(0, gas::prim::density(0), ia[0], ia[1], ia[2]); + + // Extrapolate gas velocity + Real gva[3] = {v(0, gas::prim::velocity(0), ia[0], ia[1], ia[2]), + v(0, gas::prim::velocity(1), ia[0], ia[1], ia[2]), + v(0, gas::prim::velocity(2), ia[0], ia[1], ia[2])}; + Real gvp1[3] = {v(0, gas::prim::velocity(0), ip1[0], ip1[1], ip1[2]), + v(0, gas::prim::velocity(1), ip1[0], ip1[1], ip1[2]), + v(0, gas::prim::velocity(2), ip1[0], ip1[1], ip1[2])}; + Real gvm1[3] = {v(0, gas::prim::velocity(0), im1[0], im1[1], im1[2]), + v(0, gas::prim::velocity(1), im1[0], im1[1], im1[2]), + v(0, gas::prim::velocity(2), im1[0], im1[1], im1[2])}; + const Real gvp = ArtemisUtils::VDot(gva, epa) + dp.omf * xcyla[0]; + const Real gvz = ArtemisUtils::VDot(gva, eza); + const Real gvp1p = ArtemisUtils::VDot(gvp1, epp1) + dp.omf * xcylp1[0]; + const Real gvm1p = ArtemisUtils::VDot(gvm1, epm1) + dp.omf * xcylm1[0]; + const Real dgvp = std::log(gvp1p / gvm1p); + const Real vpg = gvp * std::exp(dgvp * xmadx); + Real rhog, gvR; + if constexpr (BDY == IndexDomain::inner_x1) { + rhog = rhoa * nua / nug; + gvR = -1.5 * nug / xcyl[0]; + } else if constexpr (BDY == IndexDomain::outer_x1) { + // dFnu/dl = Mdot + const Real lg = xcyl[0] * vpg; + const Real la = xcyla[0] * gvp; + rhog = (3.0 * M_PI * rhoa * nua * la + dp.mdot * (lg - la)) / + (3.0 * M_PI * nug * lg); + gvR = -dp.mdot / (2 * M_PI * xcyl[0] * rhog); + } + + const Real gvcyl[3] = {gvR, vpg - dp.omf * xcyl[0], gvz}; + const Real gvel[3] = {ArtemisUtils::VDot(gvcyl, ex1), + ArtemisUtils::VDot(gvcyl, ex2), + ArtemisUtils::VDot(gvcyl, ex3)}; + // Set extrapolated values + v(0, gas::prim::density(0), k, j, i) = rhog; + v(0, gas::prim::sie(0), k, j, i) = sieg; + v(0, gas::prim::velocity(ix1), k, j, i) = gvel[ix1]; + v(0, gas::prim::velocity(ix2), k, j, i) = gvel[ix2]; + v(0, gas::prim::velocity(ix3), k, j, i) = gvel[ix3]; + + if (do_dust) { + for (int n = 0; n < v.GetSize(0, dust::prim::density()); ++n) { + // Extrapolate dust density + Real ddrho = std::log(v(0, dust::prim::density(n), ip1[0], ip1[1], ip1[2]) / + v(0, dust::prim::density(n), im1[0], im1[1], im1[2])); + const Real drhoexp = std::exp(ddrho * xmadx); + const Real rhod = v(0, dust::prim::density(n), ia[0], ia[1], ia[2]) * drhoexp; + + // Extrapolate dust velocity + Real dva[3] = {v(0, dust::prim::velocity(VI(n, 0)), ia[0], ia[1], ia[2]), + v(0, dust::prim::velocity(VI(n, 1)), ia[0], ia[1], ia[2]), + v(0, dust::prim::velocity(VI(n, 2)), ia[0], ia[1], ia[2])}; + Real dvp1[3] = {v(0, dust::prim::velocity(VI(n, 0)), ip1[0], ip1[1], ip1[2]), + v(0, dust::prim::velocity(VI(n, 1)), ip1[0], ip1[1], ip1[2]), + v(0, dust::prim::velocity(VI(n, 2)), ip1[0], ip1[1], ip1[2])}; + Real dvm1[3] = {v(0, dust::prim::velocity(VI(n, 0)), im1[0], im1[1], im1[2]), + v(0, dust::prim::velocity(VI(n, 1)), im1[0], im1[1], im1[2]), + v(0, dust::prim::velocity(VI(n, 2)), im1[0], im1[1], im1[2])}; + const Real dvp = ArtemisUtils::VDot(dva, epa) + dp.omf * xcyla[0]; + const Real dvR = ArtemisUtils::VDot(dva, eRa); + const Real dvz = ArtemisUtils::VDot(dva, eza); + const Real dvp1p = ArtemisUtils::VDot(dvp1, epp1) + dp.omf * xcylp1[0]; + const Real dvm1p = ArtemisUtils::VDot(dvm1, epm1) + dp.omf * xcylm1[0]; + const Real ddvp = std::log(dvp1p / dvm1p); + const Real dvcyl[3] = {dvR, dvp * std::exp(dgvp * xmadx) - dp.omf * xcyl[0], + dvz}; + const Real dvel[3] = {ArtemisUtils::VDot(dvcyl, ex1), + ArtemisUtils::VDot(dvcyl, ex2), + ArtemisUtils::VDot(dvcyl, ex3)}; + + // Set extrapolated values + v(0, dust::prim::density(n), k, j, i) = rhod; + v(0, dust::prim::velocity(VI(n, ix1)), k, j, i) = dvel[ix1]; + v(0, dust::prim::velocity(VI(n, ix2)), k, j, i) = dvel[ix2]; + v(0, dust::prim::velocity(VI(n, ix3)), k, j, i) = dvel[ix3]; + } + } + }); +} + +//---------------------------------------------------------------------------------------- +//! \fn void Disk::DiskBoundaryIC() +//! \brief Sets inner or outer boundary condition to the initial condition +template +void DiskBoundaryIC(std::shared_ptr> &mbd, bool coarse) { + auto pmb = mbd->GetBlockPointer(); + + auto artemis_pkg = pmb->packages.Get("artemis"); + auto disk_params = artemis_pkg->Param("disk_params"); + + auto &gas_pkg = pmb->packages.Get("gas"); + auto eos_d = gas_pkg->template Param("eos_d"); + + static auto descriptors = + ArtemisUtils::GetBoundaryPackDescriptorMap(mbd); + + auto v = descriptors[coarse].GetPack(mbd.get()); + + const auto &pco = (coarse) ? pmb->pmr->GetCoarseCoords() : pmb->coords; + auto &dp = disk_params; + const auto nb = IndexRange{0, 0}; + const bool fine = false; + + ParArray1D particles; + int npart = 0; + if (dp.nbody_temp) { + auto &nbody_pkg = pmb->packages.Get("nbody"); + particles = nbody_pkg->template Param>("particles"); + npart = static_cast(particles.size()); + } + + pmb->par_for_bndry( + "DiskInnerX1", nb, BDY, parthenon::TopologicalElement::CC, coarse, fine, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + DiskICImpl(v, 0, k, j, i, pco, eos_d, dp, particles, npart); + }); +} + +//---------------------------------------------------------------------------------------- +//! \fn void Disk::DiskBoundaryExtrap() +//! \brief Extrapolation boundary conditions +template +void DiskBoundaryExtrap(std::shared_ptr> &mbd, bool coarse) { + const bool lnx = (GEOM != Coordinates::cartesian); + + auto pmb = mbd->GetBlockPointer(); + + // Extract artemis parameters + auto artemis_pkg = pmb->packages.Get("artemis"); + const bool do_dust = artemis_pkg->Param("do_dust"); + + // Extract gas parameters + auto &gas_pkg = pmb->packages.Get("gas"); + auto eos_d = gas_pkg->template Param("eos_d"); + + auto disk_params = artemis_pkg->Param("disk_params"); + auto &dp = disk_params; + + // Packing + static auto descriptors = + ArtemisUtils::GetBoundaryPackDescriptorMap(mbd); + + auto v = descriptors[coarse].GetPack(mbd.get()); + + const auto &pco = (coarse) ? pmb->pmr->GetCoarseCoords() : pmb->coords; + const auto nb = IndexRange{0, 0}; + const bool fine = false; + + // Boundary index arithmetic + int is = Null(), ie = Null(); + int js = Null(), je = Null(); + int ks = Null(), ke = Null(); + const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds; + if constexpr (BDY == IndexDomain::inner_x1) { + is = bounds.GetBoundsI(IndexDomain::interior, TE::CC).s; + } else if constexpr (BDY == IndexDomain::outer_x1) { + ie = bounds.GetBoundsI(IndexDomain::interior, TE::CC).e; + } else if constexpr (BDY == IndexDomain::inner_x2) { + js = bounds.GetBoundsJ(IndexDomain::interior, TE::CC).s; + } else if constexpr (BDY == IndexDomain::outer_x2) { + je = bounds.GetBoundsJ(IndexDomain::interior, TE::CC).e; + } else if constexpr (BDY == IndexDomain::inner_x3) { + ks = bounds.GetBoundsK(IndexDomain::interior, TE::CC).s; + } else if constexpr (BDY == IndexDomain::outer_x3) { + ke = bounds.GetBoundsK(IndexDomain::interior, TE::CC).e; + } + constexpr bool x1dir = + ((BDY == IndexDomain::inner_x1) || (BDY == IndexDomain::outer_x1)); + constexpr bool x2dir = + ((BDY == IndexDomain::inner_x2) || (BDY == IndexDomain::outer_x2)); + constexpr bool x3dir = + ((BDY == IndexDomain::inner_x3) || (BDY == IndexDomain::outer_x3)); + constexpr int ix1 = x1dir ? 0 : (x2dir ? 1 : 2); + constexpr int ix2 = (ix1 + 1) % 3; + constexpr int ix3 = (ix1 + 2) % 3; + + pmb->par_for_bndry( + "DiskExtrap", nb, BDY, parthenon::TopologicalElement::CC, coarse, fine, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + // We are extrapolating into the ghost zone. Extrapolation is done on cylinders. + // dP/dz = - rho grad(Phi) + // rho vp^2/R = rho vk^2/R + dP/dR + // We estimate the pressure gradients as grad(P).eR and grad(P).ez + // Vertical hydrostatic balance sets rho + // Radial centrifugal balance sets vp + const int ia[3] = {x3dir ? ((BDY == IndexDomain::inner_x3) ? ks : ke) : k, + x2dir ? ((BDY == IndexDomain::inner_x2) ? js : je) : j, + x1dir ? ((BDY == IndexDomain::inner_x1) ? is : ie) : i}; + const int ip1[3] = {x3dir ? ((BDY == IndexDomain::inner_x3) ? ks + 1 : ke) : k, + x2dir ? ((BDY == IndexDomain::inner_x2) ? js + 1 : je) : j, + x1dir ? ((BDY == IndexDomain::inner_x1) ? is + 1 : ie) : i}; + const int im1[3] = {x3dir ? ((BDY == IndexDomain::inner_x3) ? ks : ke - 1) : k, + x2dir ? ((BDY == IndexDomain::inner_x2) ? js : je - 1) : j, + x1dir ? ((BDY == IndexDomain::inner_x1) ? is : ie - 1) : i}; + + // Extract coordinates at k, j, i + geometry::Coords coords(pco, k, j, i); + const auto &xv = coords.GetCellCenter(); + const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + + // Extract coordinates at ia, im, ic + geometry::Coords ca(pco, ia[0], ia[1], ia[2]); + geometry::Coords cp1(pco, ip1[0], ip1[1], ip1[2]); + geometry::Coords cm1(pco, im1[0], im1[1], im1[2]); + const auto &xva = ca.GetCellCenter(); + const auto &[xcyla, scr1, scr2, scr3] = ca.ConvertToCylWithVec(xva); + const Real eRa[3] = {scr1[0], scr2[0], scr3[0]}; + const Real epa[3] = {scr1[1], scr2[1], scr3[1]}; + const Real eza[3] = {scr1[2], scr2[2], scr3[2]}; + const auto &xvp1 = cp1.GetCellCenter(); + const auto &[xcylp1, scr1p1, scr2p1, scr3p1] = cp1.ConvertToCylWithVec(xvp1); + const Real epp1[3] = {scr1p1[1], scr2p1[1], scr3p1[1]}; + const auto &xvm1 = cm1.GetCellCenter(); + const auto &[xcylm1, scr1m1, scr2m1, scr3m1] = cm1.ConvertToCylWithVec(xvm1); + const Real epm1[3] = {scr1m1[1], scr2m1[1], scr3m1[1]}; + + // Compute cell separations (using logarithmics if necessary) + const Real xma = (lnx) ? std::log(xv[ix1] / xva[ix1]) : xv[ix1] - xva[ix1]; + const Real dx = (lnx) ? std::log(xvp1[ix1] / xvm1[ix1]) : xvp1[ix1] - xvm1[ix1]; + const Real xmadx = xma / dx; + + // Extrapolate gas density and specific internal energy + Real dgrho = std::log(v(0, gas::prim::density(0), ip1[0], ip1[1], ip1[2]) / + v(0, gas::prim::density(0), im1[0], im1[1], im1[2])); + Real dgsie = std::log(v(0, gas::prim::sie(0), ip1[0], ip1[1], ip1[2]) / + v(0, gas::prim::sie(0), im1[0], im1[1], im1[2])); + const Real grhoexp = std::exp(dgrho * xmadx); + const Real gsieexp = std::exp(dgsie * xmadx); + const Real rhog = v(0, gas::prim::density(0), ia[0], ia[1], ia[2]) * grhoexp; + const Real sieg = v(0, gas::prim::sie(0), ia[0], ia[1], ia[2]) * gsieexp; + + // Extrapolate gas velocity + Real gva[3] = {v(0, gas::prim::velocity(0), ia[0], ia[1], ia[2]), + v(0, gas::prim::velocity(1), ia[0], ia[1], ia[2]), + v(0, gas::prim::velocity(2), ia[0], ia[1], ia[2])}; + Real gvp1[3] = {v(0, gas::prim::velocity(0), ip1[0], ip1[1], ip1[2]), + v(0, gas::prim::velocity(1), ip1[0], ip1[1], ip1[2]), + v(0, gas::prim::velocity(2), ip1[0], ip1[1], ip1[2])}; + Real gvm1[3] = {v(0, gas::prim::velocity(0), im1[0], im1[1], im1[2]), + v(0, gas::prim::velocity(1), im1[0], im1[1], im1[2]), + v(0, gas::prim::velocity(2), im1[0], im1[1], im1[2])}; + const Real gvp = ArtemisUtils::VDot(gva, epa) + dp.omf * xcyla[0]; + const Real gvR = ArtemisUtils::VDot(gva, eRa); + const Real gvz = ArtemisUtils::VDot(gva, eza); + const Real gvp1p = ArtemisUtils::VDot(gvp1, epp1) + dp.omf * xcylp1[0]; + const Real gvm1p = ArtemisUtils::VDot(gvm1, epm1) + dp.omf * xcylm1[0]; + const Real dgvp = std::log(gvp1p / gvm1p); + const Real gvcyl[3] = {gvR, gvp * std::exp(dgvp * xmadx) - dp.omf * xcyl[0], gvz}; + const Real gvel[3] = {ArtemisUtils::VDot(gvcyl, ex1), + ArtemisUtils::VDot(gvcyl, ex2), + ArtemisUtils::VDot(gvcyl, ex3)}; + + // Set extrapolated values + v(0, gas::prim::density(0), k, j, i) = rhog; + v(0, gas::prim::sie(0), k, j, i) = sieg; + v(0, gas::prim::velocity(ix1), k, j, i) = gvel[ix1]; + v(0, gas::prim::velocity(ix2), k, j, i) = gvel[ix2]; + v(0, gas::prim::velocity(ix3), k, j, i) = gvel[ix3]; + + if (do_dust) { + for (int n = 0; n < v.GetSize(0, dust::prim::density()); ++n) { + // Extrapolate dust density + Real ddrho = std::log(v(0, dust::prim::density(n), ip1[0], ip1[1], ip1[2]) / + v(0, dust::prim::density(n), im1[0], im1[1], im1[2])); + const Real drhoexp = std::exp(ddrho * xmadx); + const Real rhod = v(0, dust::prim::density(n), ia[0], ia[1], ia[2]) * drhoexp; + + // Extrapolate dust velocity + Real dva[3] = {v(0, dust::prim::velocity(VI(n, 0)), ia[0], ia[1], ia[2]), + v(0, dust::prim::velocity(VI(n, 1)), ia[0], ia[1], ia[2]), + v(0, dust::prim::velocity(VI(n, 2)), ia[0], ia[1], ia[2])}; + Real dvp1[3] = {v(0, dust::prim::velocity(VI(n, 0)), ip1[0], ip1[1], ip1[2]), + v(0, dust::prim::velocity(VI(n, 1)), ip1[0], ip1[1], ip1[2]), + v(0, dust::prim::velocity(VI(n, 2)), ip1[0], ip1[1], ip1[2])}; + Real dvm1[3] = {v(0, dust::prim::velocity(VI(n, 0)), im1[0], im1[1], im1[2]), + v(0, dust::prim::velocity(VI(n, 1)), im1[0], im1[1], im1[2]), + v(0, dust::prim::velocity(VI(n, 2)), im1[0], im1[1], im1[2])}; + const Real dvp = ArtemisUtils::VDot(dva, epa) + dp.omf * xcyla[0]; + const Real dvR = ArtemisUtils::VDot(dva, eRa); + const Real dvz = ArtemisUtils::VDot(dva, eza); + const Real dvp1p = ArtemisUtils::VDot(dvp1, epp1) + dp.omf * xcylp1[0]; + const Real dvm1p = ArtemisUtils::VDot(dvm1, epm1) + dp.omf * xcylm1[0]; + const Real ddvp = std::log(dvp1p / dvm1p); + const Real dvcyl[3] = {dvR, dvp * std::exp(dgvp * xmadx) - dp.omf * xcyl[0], + dvz}; + const Real dvel[3] = {ArtemisUtils::VDot(dvcyl, ex1), + ArtemisUtils::VDot(dvcyl, ex2), + ArtemisUtils::VDot(dvcyl, ex3)}; + + // Set extrapolated values + v(0, dust::prim::density(n), k, j, i) = rhod; + v(0, dust::prim::velocity(VI(n, ix1)), k, j, i) = dvel[ix1]; + v(0, dust::prim::velocity(VI(n, ix2)), k, j, i) = dvel[ix2]; + v(0, dust::prim::velocity(VI(n, ix3)), k, j, i) = dvel[ix3]; + } + } + }); +} + +//---------------------------------------------------------------------------------------- +//! \fn AmrTag ProblemCheckRefinementBlock() +//! \brief Refinement criterion for disk pgen +inline parthenon::AmrTag ProblemCheckRefinementBlock(MeshBlockData *mbd) { + PARTHENON_FAIL("Disk user-defined AMR criterion not yet implemented!"); + return AmrTag::same; +} + +} // namespace disk + +#endif // PGEN_DISK_HPP_ From c4862578f5184cc6636d8fe7330b0980367e0d69 Mon Sep 17 00:00:00 2001 From: jtlaune Date: Tue, 10 Dec 2024 09:46:29 -0700 Subject: [PATCH 2/6] Add pgen files for draft PR --- src/CMakeLists.txt | 1 + src/pgen/pgen.hpp | 3 +++ src/pgen/problem_modifier.hpp | 38 +++++++++++++++++++++++++++++++++++ 3 files changed, 42 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 238c5ac..ed26c34 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -64,6 +64,7 @@ set (SRC_LIST pgen/conduction.hpp pgen/constant.hpp pgen/disk.hpp + pgen/nu_disk.hpp pgen/gaussian_bump.hpp pgen/linear_wave.hpp pgen/pgen.hpp diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index abf5887..86b04ea 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -23,6 +23,7 @@ #include "conduction.hpp" #include "constant.hpp" #include "disk.hpp" +#include "nu_disk.hpp" #include "gaussian_bump.hpp" #include "linear_wave.hpp" #include "shock.hpp" @@ -48,6 +49,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { constant::ProblemGenerator(pmb, pin); } else if (name == "disk") { disk::ProblemGenerator(pmb, pin); + } else if (name == "nu_disk") { + nu_disk::ProblemGenerator(pmb, pin); } else if (name == "gaussian_bump") { gaussian_bump::ProblemGenerator(pmb, pin); } else if (name == "linear_wave") { diff --git a/src/pgen/problem_modifier.hpp b/src/pgen/problem_modifier.hpp index 1b93dee..0839d75 100644 --- a/src/pgen/problem_modifier.hpp +++ b/src/pgen/problem_modifier.hpp @@ -102,6 +102,44 @@ void ProblemModifier(parthenon::ParthenonManager *pman) { pman->app_input->RegisterBoundaryCondition(BF::outer_x1, "viscous", disk::DiskBoundaryVisc); } + } else if (artemis_problem == "nu_disk") { + pman->app_input->InitMeshBlockUserData = nu_disk::InitDiskParams; + + artemis::ProblemCheckRefinementBlock = nu_disk::ProblemCheckRefinementBlock; + + pman->app_input->RegisterBoundaryCondition(BF::inner_x1, "ic", + nu_disk::DiskBoundaryIC); + pman->app_input->RegisterBoundaryCondition(BF::outer_x1, "ic", + nu_disk::DiskBoundaryIC); + pman->app_input->RegisterBoundaryCondition(BF::inner_x2, "ic", + nu_disk::DiskBoundaryIC); + pman->app_input->RegisterBoundaryCondition(BF::outer_x2, "ic", + nu_disk::DiskBoundaryIC); + pman->app_input->RegisterBoundaryCondition(BF::inner_x3, "ic", + nu_disk::DiskBoundaryIC); + pman->app_input->RegisterBoundaryCondition(BF::outer_x3, "ic", + nu_disk::DiskBoundaryIC); + + pman->app_input->RegisterBoundaryCondition(BF::inner_x1, "extrap", + nu_disk::DiskBoundaryExtrap); + pman->app_input->RegisterBoundaryCondition(BF::outer_x1, "extrap", + nu_disk::DiskBoundaryExtrap); + pman->app_input->RegisterBoundaryCondition(BF::inner_x2, "extrap", + nu_disk::DiskBoundaryExtrap); + pman->app_input->RegisterBoundaryCondition(BF::outer_x2, "extrap", + nu_disk::DiskBoundaryExtrap); + pman->app_input->RegisterBoundaryCondition(BF::inner_x3, "extrap", + nu_disk::DiskBoundaryExtrap); + pman->app_input->RegisterBoundaryCondition(BF::outer_x3, "extrap", + nu_disk::DiskBoundaryExtrap); + + if constexpr (geometry::is_axisymmetric() || (G == Coordinates::cylindrical) || + (G == Coordinates::spherical3D)) { + pman->app_input->RegisterBoundaryCondition(BF::inner_x1, "viscous", + nu_disk::DiskBoundaryVisc); + pman->app_input->RegisterBoundaryCondition(BF::outer_x1, "viscous", + nu_disk::DiskBoundaryVisc); + } } else if (artemis_problem == "linear_wave") { pman->app_input->UserWorkAfterLoop = linear_wave::UserWorkAfterLoop; } else if (artemis_problem == "shock") { From a1c1fc3c45766988f3c86a2d87812f187b54f70c Mon Sep 17 00:00:00 2001 From: jtlaune Date: Tue, 10 Dec 2024 12:22:44 -0700 Subject: [PATCH 3/6] Fix compile errors, missing bracket --- src/pgen/nu_disk.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pgen/nu_disk.hpp b/src/pgen/nu_disk.hpp index 22b7dfc..cc0991c 100644 --- a/src/pgen/nu_disk.hpp +++ b/src/pgen/nu_disk.hpp @@ -156,6 +156,7 @@ ComputeDiskProfile(const struct DiskParams pgen, const parthenon::Coordinates_t gvel1 = ArtemisUtils::VDot(vcyl, ex1); gvel2 = ArtemisUtils::VDot(vcyl, ex2); gvel3 = ArtemisUtils::VDot(vcyl, ex3); +} //---------------------------------------------------------------------------------------- //! \fn void InitDiskParams From 99379f4fa3fcc0723face7d846078baafb44b902 Mon Sep 17 00:00:00 2001 From: jtlaune Date: Tue, 10 Dec 2024 12:23:59 -0700 Subject: [PATCH 4/6] Remove unnecessary PI_2 definition --- src/drag/drag.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/drag/drag.hpp b/src/drag/drag.hpp index e878305..bc61881 100644 --- a/src/drag/drag.hpp +++ b/src/drag/drag.hpp @@ -213,8 +213,6 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt const auto &hx = coords.GetScaleFactors(); const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); - const Real PI_2 = 1.5707963267948966; - // Compute the ramp for this cell // Ramps are quadratic, eg. the left regions is SQR( (X - ix)/(ix - xmin) ) if (do_gas) { From c1fc351aa3c1c1b313f73d07f9951a71e2eac60a Mon Sep 17 00:00:00 2001 From: jtlaune Date: Mon, 13 Jan 2025 17:39:16 -0700 Subject: [PATCH 5/6] Merge in stash --- env/bash | 60 +++++++++++---------- src/drag/drag.hpp | 59 +++++++++++++++++++- src/pgen/disk.hpp | 133 ++++++++++++---------------------------------- 3 files changed, 125 insertions(+), 127 deletions(-) diff --git a/env/bash b/env/bash index 88f8b01..b78fc57 100644 --- a/env/bash +++ b/env/bash @@ -19,35 +19,39 @@ export ARTEMIS_HOME="${exec_dir%/env}" export MAKE_PROGRAM=${MAKE_PROGRAM:-make} +export PARTITION="venado-gh" + # Identify partitions based on SLURM_JOB_PARTITION and HOSTNAME variables -PARTITION="unknown" -if [[ $HOSTNAME == ch-fe* || $SLURM_CLUSTER_NAME == "chicoma" ]]; then - if [[ $SLURM_JOB_PARTITION == *gpu* ]]; then - PARTITION="chicoma-gpu" - else - PARTITION="chicoma-cpu" - fi -elif [[ "$HOST" =~ ^ve-rfe[1-7]$ || "$HOST" =~ ^ve-fe[1-7]$ ]]; then - if [[ $SLURM_GPUS_ON_NODE > 0 || ("$HOST" =~ ^ve-rfe[1-3]$ || "$HOST" =~ ^ve-fe[1-3]$) ]]; then - PARTITION="venado-gh" - elif [[ $SLURM_GPUS_ON_NODE == 0 || ("$HOST" =~ ^ve-rfe[4-7]$ || "$HOST" =~ ^ve-fe[4-7]$) ]]; then - PARTITION="venado-gg" - fi -else # Catch-all for Darwin - if [[ $HOSTNAME == darwin-fe* ]]; then - echo "Currently on darwin frontend node; formatting is supported but not compilation" - echo "Supported partitions are" - echo " skylake-gold" - echo " volta-x86" - PARTITION="darwin-fe" - elif [[ $SLURM_JOB_PARTITION == "power9-rhel7" ]]; then - PARTITION="darwin-power9-rhel7" - elif [[ $SLURM_JOB_PARTITION == "skylake-gold" ]]; then - PARTITION="darwin-skylake-gold" - elif [[ $SLURM_JOB_PARTITION == "volta-x86" ]]; then - PARTITION="darwin-volta-x86" - fi -fi +#PARTITION="unknown" +#if [[ $HOSTNAME == ch-fe* || $SLURM_CLUSTER_NAME == "chicoma" ]]; then +# if [[ $SLURM_JOB_PARTITION == *gpu* ]]; then +# PARTITION="chicoma-gpu" +# else +# PARTITION="chicoma-cpu" +# fi +#elif [[ "$HOSTNAME" =~ ^ve-rfe[1-7]$ || "$HOSTNAME" =~ ^ve-fe[1-7]$ ]]; then +# if [[ $SLURM_GPUS_ON_NODE > 0 || ("$HOSTNAME" =~ ^ve-rfe[1-3]$ || "$HOSTNAME" =~ ^ve-fe[1-3]$) ]]; then +# PARTITION="venado-gh" +# elif [[ $SLURM_GPUS_ON_NODE == 0 || ("$HOSTNAME" =~ ^ve-rfe[4-7]$ || "$HOSTNAME" =~ ^ve-fe[4-7]$) ]]; then +# PARTITION="venado-gg" +# fi +#else # Catch-all for Darwin +# if [ -z "$SLURM_JOB_PARTITION" ]; then +# if [[ $HOSTNAME == darwin-fe* ]]; then +# echo "Do not compile on Darwin frontend nodes; use a specific partition!" +# echo "Supported partitions are" +# echo " skylake-gold" +# echo " volta-x86" +# echo "...setup FAILED" +# fi +# elif [[ $SLURM_JOB_PARTITION == "power9-rhel7" ]]; then +# PARTITION="darwin-power9-rhel7" +# elif [[ $SLURM_JOB_PARTITION == "skylake-gold" ]]; then +# PARTITION="darwin-skylake-gold" +# elif [[ $SLURM_JOB_PARTITION == "volta-x86" ]]; then +# PARTITION="darwin-volta-x86" +# fi +#fi # Absolute path for conda environment can be too long function shorten_prompt { diff --git a/src/drag/drag.hpp b/src/drag/drag.hpp index bc61881..71489d5 100644 --- a/src/drag/drag.hpp +++ b/src/drag/drag.hpp @@ -213,6 +213,8 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt const auto &hx = coords.GetScaleFactors(); const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + const Real PI_2 = 1.5707963267948966; + // Compute the ramp for this cell // Ramps are quadratic, eg. the left regions is SQR( (X - ix)/(ix - xmin) ) if (do_gas) { @@ -221,7 +223,11 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt dt * (gasp.irate[0] * ((xv[0] < gasp.ix[0]) * SQR((xv[0] - gasp.ix[0]) / (gasp.ix[0] - x1min))) + gasp.orate[0] * ((xv[0] > gasp.ox[0]) * - SQR((xv[0] - gasp.ox[0]) / (gasp.ox[0] - x1max)))); + SQR((xv[0] - gasp.ox[0]) / (gasp.ox[0] - x1max))));// + + //gasp.irate[1] * ((xv[0] >= gasp.ix[0]) * (xcyl[2] > gasp.ix[1]*H) * // pos z + // SQR((xcyl[2] - gasp.ix[1]*H) / (gasp.ix[0]*H - xv[0]*std::cos(x2min)))) + + //gasp.orate[1] * ((xv[0] <= gasp.ox[0]) * (xcyl[2] < -gasp.ox[1]*H) * // neg z + // SQR((xcyl[2] + gasp.ox[1]*H) / (-gasp.ox[0]*H - xv[0]*std::cos(x2max))))); const Real fx2 = multi_d * dt * (gasp.irate[1] * ((xv[1] < gasp.ix[1]) * @@ -262,6 +268,40 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt ArtemisUtils::VDot(vcyl, ex2), ArtemisUtils::VDot(vcyl, ex3) }; + //Real vd[3] = { + // vR, + // 0., + // vp - omf * xcyl[0], + //}; + + //if (i==2 && j==2 && k==2) { + // std::cout << "(" << i << ", " << j << ", " << k << ")" << std::endl; + // std::cout << i << j << k << std::endl; + // std::cout << "drag vd:" << std::fixed << std::setprecision(12) << vd[0] + // << ", " << vd[1] << ", " << vd[2] << std::endl; + // std::cout << "drag fx1:" << fx1 << std::endl; + //} + //if (i==2 && j==32 && k==2) { + // std::cout << "(" << i << ", " << j << ", " << k << ")" << std::endl; + // std::cout << "drag vd:" << std::fixed << std::setprecision(12) << vd[0] + // << ", " << vd[1] << ", " << vd[2] << std::endl; + // std::cout << "drag fx1:" << fx1 << std::endl; + //} + //if (i==64 && j==2 && k==2) { + // std::cout << "(" << i << ", " << j << ", " << k << ")" << std::endl; + // std::cout << "drag vd:" << std::fixed << std::setprecision(12) << vd[0] + // << ", " << vd[1] << ", " << vd[2] << std::endl; + // std::cout << "drag fx1:" << fx1 << std::endl; + //} + //if (i==64 && j==32 && k==2) { + // std::cout << "(" << i << ", " << j << ", " << k << ")" << std::endl; + // std::cout << "drag vd:" << std::fixed << std::setprecision(12) << vd[0] + // << ", " << vd[1] << ", " << vd[2] << std::endl; + // std::cout << "drag fx1:" << fx1 << std::endl; + //} + // debugging code + + // std::cout << "drag:" << i << ", " << j << ", " << k << std::endl; // Ep - E = 0.5 d ( vp^2 - v^2 ) // (vp-v) . (vp + v) = dv . (2v + dv) = 2 dv.v + dv.dv @@ -273,6 +313,23 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) += hx[1] * dm2; vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) += hx[2] * dm3; + //if (i==2 && j==2 && k==2) { + // std::cout << std::fixed << std::setprecision(12) << "drag nu: " << nu << std::endl; + // std::cout << "gas vel after damping:" + // << std::fixed << std::setprecision(12) + // << vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i)/dens << ", " + // << vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i)/dens << ", " + // << vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i)/dens << std::endl; + // std::cout << "drag cc:" << std::fixed << std::setprecision(12) << xv[0] + // << ", " << xv[1] << ", " << xv[2] << std::endl; + // std::cout << "drag vd:" << std::fixed << std::setprecision(12) << vd[0] + // << ", " << vd[1] << ", " << vd[2] << std::endl; + // std::cout << "drag vcyl:" << std::fixed << std::setprecision(12) << vcyl[0] + // << ", " << vcyl[1] << ", " << vcyl[2] << std::endl; + // std::cout << "drag ex1:" << std::fixed << std::setprecision(12)<< ex1[0] + // << ", " << ex1[1] << ", " << ex1[2] << std::endl; + //} + vmesh(b, gas::cons::total_energy(n), k, j, i) += dm1 * (vg[0] + 0.5 * dm1 / dens) + dm2 * (vg[1] + 0.5 * dm2 / dens) + dm3 * (vg[2] + 0.5 * dm3 / dens); diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp index 7c34ea3..c08e5fe 100644 --- a/src/pgen/disk.hpp +++ b/src/pgen/disk.hpp @@ -60,6 +60,7 @@ struct DiskParams { Real alpha, nu0, nu_indx; Real mdot; Real temp_soft2; + bool const_nu; bool do_dust; bool nbody_temp; bool quiet_start; @@ -70,23 +71,8 @@ struct DiskParams { //! \brief Computes density profile at cylindrical R and z KOKKOS_INLINE_FUNCTION Real DenProfile(struct DiskParams pgen, const Real R, const Real z) { - const Real r = std::sqrt(R * R + z * z); - const Real h = pgen.h0 * std::pow(R / pgen.r0, pgen.flare); - const Real sig0 = pgen.rho0; // / (std::sqrt(2.0 * M_PI) * pgen.h0 * pgen.r0); - const Real exp_fac = (pgen.rexp == 0.) ? 1. : std::exp(-SQR(R / pgen.rexp)); - const Real dmid = - (sig0 * std::pow(R / pgen.r0, pgen.p)) * - (1. - pgen.l0 * std::sqrt(pgen.r0 / R)) * // correction for an inner binary - (pgen.dens_min / pgen.rho0 + // The inner cavity - (1. - pgen.dens_min / pgen.rho0) * std::exp(-std::pow(pgen.rcav / R, 12.0))) * - exp_fac; // the outer cutoff - const Real sint = (r == 0.0) ? 1.0 : R / r; // TODO(ADM): should it be 1? - const Real efac = (1. - sint) / (h * h); - if (pgen.Gamma == 1.) return std::max(pgen.dens_min, dmid * std::exp(-efac)); - // sint <= h^2/(gamma-1), efac*(g-1) = 1 - eps - const Real pfac = 1. - (pgen.Gamma - 1) * efac; - return std::max(pgen.dens_min, - dmid * std::pow(pfac + Fuzz(), 1. / (pgen.Gamma - 1))); + const Real H = R * pgen.h0 * std::pow(R / pgen.r0, pgen.flare); + return std::max(pgen.dens_min, pgen.rho0*std::pow(R/pgen.r0, pgen.p)*std::exp(-0.5*SQR(z/H))); } //---------------------------------------------------------------------------------------- @@ -148,98 +134,47 @@ ComputeDiskProfile(const struct DiskParams pgen, const parthenon::Coordinates_t : xcyl[0]; gtemp = TempProfile(pgen, rt, xcyl[2]); - // Construct grad(P) for this zone - const std::array fx1m{coords.bnds.x1[0], xv[1], xv[2]}; - const std::array fx1p{coords.bnds.x1[1], xv[1], xv[2]}; - const std::array fx2m{xv[0], coords.bnds.x2[0], xv[2]}; - const std::array fx2p{xv[0], coords.bnds.x2[1], xv[2]}; - const std::array fx3m{xv[0], xv[1], coords.bnds.x3[0]}; - const std::array fx3p{xv[0], xv[1], coords.bnds.x3[1]}; - Real pgrad[3] = {Null()}; - Real pfm = Null(), pfp = Null(); - - // X1 Faces - auto xf = coords.ConvertToCyl(fx1m); - Real rtm = pgen.nbody_temp ? -pgen.gm / Gravity::NBodyPotential(coords, fx1m, - particles, npart) - : xf[0]; - Real tfm = TempProfile(pgen, rtm, xf[2]); - pfm = PresProfile(pgen, eos_d, tfm, xf[0], xf[2]); - xf = coords.ConvertToCyl(fx1p); - Real rtp = pgen.nbody_temp ? -pgen.gm / Gravity::NBodyPotential(coords, fx1p, - particles, npart) - : xf[0]; - Real tfp = TempProfile(pgen, rtp, xf[2]); - pfp = - (pfm = pgen.pres_min) ? pgen.pres_min : PresProfile(pgen, eos_d, tfp, xf[0], xf[2]); - pfm = (pfp == pgen.pres_min) ? pgen.pres_min : pfm; - pgrad[0] = (pfp - pfm) / coords.GetCellWidthX1(); - - // X2 Faces - xf = coords.ConvertToCyl(fx2m); - rtm = pgen.nbody_temp - ? -pgen.gm / Gravity::NBodyPotential(coords, fx2m, particles, npart) - : xf[0]; - tfm = TempProfile(pgen, rtm, xf[2]); - pfm = PresProfile(pgen, eos_d, tfm, xf[0], xf[2]); - xf = coords.ConvertToCyl(fx2p); - rtp = pgen.nbody_temp - ? -pgen.gm / Gravity::NBodyPotential(coords, fx2p, particles, npart) - : xf[0]; - tfp = TempProfile(pgen, rtp, xf[2]); - pfp = - (pfm = pgen.pres_min) ? pgen.pres_min : PresProfile(pgen, eos_d, tfp, xf[0], xf[2]); - pfm = (pfp == pgen.pres_min) ? pgen.pres_min : pfm; - pgrad[1] = (pfp - pfm) / coords.GetCellWidthX2(); - - // X3 Faces - xf = coords.ConvertToCyl(fx3m); - rtm = pgen.nbody_temp - ? -pgen.gm / Gravity::NBodyPotential(coords, fx3m, particles, npart) - : xf[0]; - tfm = TempProfile(pgen, rtm, xf[2]); - pfm = PresProfile(pgen, eos_d, tfm, xf[0], xf[2]); - xf = coords.ConvertToCyl(fx3p); - rtp = pgen.nbody_temp - ? -pgen.gm / Gravity::NBodyPotential(coords, fx3p, particles, npart) - : xf[0]; - tfp = TempProfile(pgen, rtp, xf[2]); - pfp = - (pfm = pgen.pres_min) ? pgen.pres_min : PresProfile(pgen, eos_d, tfp, xf[0], xf[2]); - pfm = (pfp == pgen.pres_min) ? pgen.pres_min : pfm; - pgrad[2] = (pfp - pfm) / coords.GetCellWidthX3(); - - // Convert to the cylindrical radial gradient - const Real eR[3] = {ex1[0], ex2[0], ex3[0]}; - const Real dpdr = ArtemisUtils::VDot(pgrad, eR); - // Set v_phi to centrifugal equilibrium // vp^2/R = grad(p) + vk^2/R - const Real r = pgen.nbody_temp ? rt : std::sqrt(SQR(xcyl[0]) + SQR(xcyl[2])); - const Real omk2 = pgen.gm / (r * r * r); - const Real vk2 = omk2 * SQR(xcyl[0]); - const Real vp = std::sqrt(vk2 + dpdr * xcyl[0] / gdens); - const Real nu = ViscosityProfile(pgen, eos_d, rt, xcyl[2]); - const Real vr = pgen.quiet_start ? 0.0 : -1.5 * nu / xcyl[0]; + if (!pgen.const_nu) { + PARTHENON_FAIL("jtlaune: This version only works with nu viscosity."); + } + + const Real H = xcyl[0] * pgen.h0 * std::pow(xcyl[0] / pgen.r0, pgen.flare); + const Real OmKmid = std::sqrt(pgen.gm / (xcyl[0] * xcyl[0] * xcyl[0])); + const Real Omg = OmKmid * (1 + 0.5 * SQR(H / xcyl[0]) * + (pgen.p + pgen.q + 0.5 * pgen.q * SQR(xcyl[2] / H))); + const Real vp = Omg * xcyl[0]; + const Real vr = - pgen.nu0*(6*pgen.p-2*pgen.q+3+(5*pgen.q+9)*SQR(xcyl[2]/H))/(2*xcyl[0]); + + const Real vz = (-pgen.p)*xcyl[2]/xcyl[0]*vr; // Construct the total cylindrical velocity - const Real vcyl[3] = {vr, vp - pgen.omf * xcyl[0], 0.0}; + const Real vcyl[3] = {vr, vp - pgen.omf * xcyl[0], vz}; // and convert it to the problem geometry gvel1 = ArtemisUtils::VDot(vcyl, ex1); gvel2 = ArtemisUtils::VDot(vcyl, ex2); gvel3 = ArtemisUtils::VDot(vcyl, ex3); - if (!(do_dust)) return; + // debugging code + + //if (i == 2 && j== 2&& k==2){ + //std::cout << "init:" << gvel1 << ", " << gvel2 << ", " << gvel3 << std::endl; + //} + //if (i==2 && j==2 && k==2){ + // std::cout << "init cc:" << std::fixed << std::setprecision(12)<< xv[0] << ", " << xv[1] << ", " << xv[2] << std::endl; + // std::cout << "init v:" << std::fixed << std::setprecision(12)<< gvel1 << ", " << gvel2 << ", " << gvel3 << std::endl; + // std::cout << "init ex1:" << std::fixed << std::setprecision(12)<< ex1[0] << ", " << ex1[1] << ", " << ex1[2] << std::endl; + // std::cout << "init vcyl:" << std::fixed << std::setprecision(12) << vcyl[0] + // << ", " << vcyl[1] << ", " << vcyl[2] << std::endl; + //} + + //if (!(do_dust)) return; - // Dust is just Keplerian - ddens = pgen.dust_to_gas * gdens; - const Real vkep[3] = {0.0, std::sqrt(vk2) - pgen.omf * xcyl[0], 0.0}; - dvel1 = ArtemisUtils::VDot(vkep, ex1); - dvel2 = ArtemisUtils::VDot(vkep, ex2); - dvel3 = ArtemisUtils::VDot(vkep, ex3); + // PARTHENON_FAIL("jtlaune: This version only works without dust."); - return; + //return; } //---------------------------------------------------------------------------------------- @@ -308,11 +243,13 @@ inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) { if (params.Get("do_viscosity")) { const auto vtype = pin->GetString("gas/viscosity", "type"); if (vtype == "alpha") { + disk_params.const_nu = false; disk_params.alpha = pin->GetReal("gas/viscosity", "alpha"); disk_params.nu0 = disk_params.alpha * disk_params.gamma_gas * SQR(disk_params.h0 * disk_params.r0 * disk_params.Omega0); disk_params.nu_indx = 1.5 + disk_params.q; - } else if (vtype == "powerlaw") { + } else if (vtype == "constant") { + disk_params.const_nu = true; disk_params.nu0 = pin->GetReal("gas/viscosity", "nu"); disk_params.nu_indx = pin->GetOrAddReal("gas/viscosity", "r_exp", 0.0); } else { From f3ba1c6b5287ac8183ca8754ea1d72c2eaa7e758 Mon Sep 17 00:00:00 2001 From: jtlaune Date: Tue, 14 Jan 2025 08:24:30 -0700 Subject: [PATCH 6/6] Fix viscosity specification in disk.hpp --- src/pgen/disk.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp index c08e5fe..3fc8ee6 100644 --- a/src/pgen/disk.hpp +++ b/src/pgen/disk.hpp @@ -248,7 +248,7 @@ inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) { disk_params.nu0 = disk_params.alpha * disk_params.gamma_gas * SQR(disk_params.h0 * disk_params.r0 * disk_params.Omega0); disk_params.nu_indx = 1.5 + disk_params.q; - } else if (vtype == "constant") { + } else if (vtype == "powerlaw") { disk_params.const_nu = true; disk_params.nu0 = pin->GetReal("gas/viscosity", "nu"); disk_params.nu_indx = pin->GetOrAddReal("gas/viscosity", "r_exp", 0.0);