diff --git a/doc/src/physics.rst b/doc/src/physics.rst index 96fad328..75e1882f 100644 --- a/doc/src/physics.rst +++ b/doc/src/physics.rst @@ -81,6 +81,63 @@ An example that triggers refinement when the gas density is larger than ``10.0`` deref_thr = 0.02 +MHD +""" + +Ideal magnetohydrodynamics is enabled by setting ``mhd = true`` in the ```` block. +When MHD is active, |code| evolves a face-centered magnetic field using constrained transport (CT). +This keeps the magnetic update tied to the divergence-preserving CT representation. + +At minimum, an MHD run should enable both gas and MHD: + +:: + + + gas = true + mhd = true + +If there are multiple gas species defined, only the first species couples to the magnetic field, i.e., its total energy contains the magnetic energy and only its velocity induces an electric fied. +Magnetic fluxes are computed through the gas Riemann solver selected in ````. +The available solver choices are: + +* ``llf`` +* ``hlle`` +* ``hlld`` + +The ``hlld`` solver is the most specialized MHD option currently available. +An example input block that enables MHD is, + +:: + + + gas = true + mhd = true + + + riemann = hlle + reconstruct = plm + + + gamma = 1.4 + +In dimensional runs, magnetic quantities use the same unit-conversion system as the rest of the code. +When ``artemis/physical_units = cgs`` and ``artemis/unit_conversion = base``, the base-unit specification now includes electrical current: + +:: + + + physical_units = cgs + unit_conversion = base + length = 1.0 + time = 1.0 + mass = 1.0 + temperature = 1.0 + current = 1.0 + +The ``current`` parameter is used together with length, time, and mass to derive the magnetic-field conversion and the code-unit value of :math:`\mu_0`. +In ``scalefree`` runs, these conversions reduce to unity. + + Gasses support several microphysics models including `Cooling`_, `Viscosity`_, and `Conduction`_. These are controlled by adding additional nodes under the ```` node, e.g., ````. We describe each of these below. diff --git a/inputs/orszag_tang/orszag_tang.in b/inputs/orszag_tang/orszag_tang.in new file mode 100644 index 00000000..83f903f7 --- /dev/null +++ b/inputs/orszag_tang/orszag_tang.in @@ -0,0 +1,84 @@ +# ======================================================================================== +# (C) (or copyright) 2023-2026. 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. +# ======================================================================================== + + +problem = orszag_tang # name of the pgen +coordinates = cartesian # coordinate system + + +problem_id = orszag_tang # problem ID: basename of output filenames + + +variables = gas.prim.density, & + gas.prim.pressure, & + gas.prim.velocity, & + field.cell.B, & + field.cell.divB +file_type = hdf5 # HDF5 data dump +dt = 0.05 # time increment between outputs + + +nlim = -1 # cycle limit +tlim = 0.5 # time limit +integrator = rk2 # time integration algorithm +ncycle_out = 1 # interval for stdout summary info + + +nghost = 2 +refinement = none +numlevel = 1 + +nx1 = 64 # Number of zones in X1-direction +x1min = 0.0 # minimum value of X1 +x1max = 1.0 # maximum value of X1 +ix1_bc = periodic # inner-X1 boundary flag +ox1_bc = periodic # outer-X1 boundary flag + +nx2 = 64 # Number of zones in X2-direction +x2min = 0.0 # minimum value of X2 +x2max = 1.0 # maximum value of X2 +ix2_bc = periodic # inner-X2 boundary flag +ox2_bc = periodic # 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 = 64 # Number of cells in each MeshBlock, X1-dir +nx2 = 64 # Number of cells in each MeshBlock, X2-dir +nx3 = 1 # Number of cells in each MeshBlock, X3-dir + + +gas = true +mhd = true +dust = false + + +cfl = 0.4 +reconstruct = plm +riemann = llf +update_fluxes = true +dfloor = 1.0e-10 +siefloor = 1.0e-10 + + +gamma = 1.66666666667 + + +rho0 = 0.22104853207207686 +p0 = 0.13262911924324611 +v0 = 1.0 +b0 = 0.28209479177387814 diff --git a/inputs/shock/brio_wu.in b/inputs/shock/brio_wu.in new file mode 100644 index 00000000..fd5aa589 --- /dev/null +++ b/inputs/shock/brio_wu.in @@ -0,0 +1,92 @@ +# ======================================================================================== +# (C) (or copyright) 2026. 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. +# ======================================================================================== + + +problem = shock # name of the pgen +coordinates = cartesian # coordinate system + + +problem_id = brio_wu # problem ID: basename of output filenames + + +variables = gas.prim.density, & + gas.prim.velocity, & + gas.prim.pressure, & + field.cell.B, & + field.cell.divB +file_type = hdf5 # HDF5 data dump +dt = 0.02 # time increment between outputs + + +nlim = -1 # cycle limit +tlim = 0.1 # time limit +integrator = rk2 # time integration algorithm +ncycle_out = 10 + + +nghost = 2 +refinement = none +numlevel = 1 + +nx1 = 512 # Number of zones in X1-direction +x1min = -0.5 # minimum value of X1 +x1max = 0.5 # maximum value of X1 +ix1_bc = ic # inner-X1 boundary flag +ox1_bc = ic # outer-X1 boundary flag + +nx2 = 1 # Number of zones in X2-direction +x2min = 0.0 # minimum value of X2 +x2max = 1.0 # maximum value of X2 +ix2_bc = periodic # inner-X2 boundary flag +ox2_bc = periodic # outer-X2 boundary flag + +nx3 = 1 # Number of zones in X3-direction +x3min = 0.0 # minimum value of X3 +x3max = 1.0 # maximum value of X3 +ix3_bc = periodic # inner-X3 boundary flag +ox3_bc = periodic # outer-X3 boundary flag + + +nx1 = 512 # Number of cells in each MeshBlock, X1-dir +nx2 = 1 # Number of cells in each MeshBlock, X2-dir +nx3 = 1 # Number of cells in each MeshBlock, X3-dir + + +gas = true +mhd = true +dust = false + + +cfl = 0.4 +reconstruct = plm +riemann = llf +update_fluxes = true +dfloor = 1.0e-10 +siefloor = 1.0e-10 + + +gamma = 2.0 + + +rhol = 1.0 +rhor = 0.125 +vxl = 0.0 +vxr = 0.0 +pl = 1.0 +pr = 0.1 +bx = 0.75 +byl = 1.0 +byr = -1.0 +bzl = 0.0 +bzr = 0.0 +xdisc = 0.0 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 51386b94..d821f2c4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -48,6 +48,9 @@ set (SRC_LIST gravity/point_mass.cpp gravity/uniform.cpp + mhd/mhd.cpp + mhd/mhd.hpp + nbody/nbody.cpp nbody/nbody.hpp nbody/nbody_advance.cpp diff --git a/src/artemis.cpp b/src/artemis.cpp index 7139c290..79bf706a 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -21,6 +21,7 @@ #include "gas/gas.hpp" #include "geometry/geometry.hpp" #include "gravity/gravity.hpp" +#include "mhd/mhd.hpp" #include "nbody/nbody.hpp" #include "radiation/moments/moments.hpp" #include "radiation/radiation.hpp" @@ -102,6 +103,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { const bool do_radiation = pin->GetOrAddBoolean("physics", "radiation", false); const bool do_coagulation = pin->GetOrAddBoolean("physics", "coagulation", false); const bool do_raytrace = pin->GetOrAddBoolean("physics", "raytrace", false); + const bool do_mhd = pin->GetOrAddBoolean("physics", "mhd", false); const bool do_orbital_advection = pin->GetOrAddBoolean("physics", "orbital_advection", false); @@ -144,6 +146,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { artemis->AddParam("do_moment", do_moment); artemis->AddParam("do_orbital_advection", do_orbital_advection); artemis->AddParam("do_raytrace", do_raytrace); + artemis->AddParam("do_mhd", do_mhd); artemis->AddParam("update_fluxes", update_fluxes); // Set coordinate system @@ -152,6 +155,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { Coordinates coords = geometry::CoordSelect(sys, ndim); artemis->AddParam("coords", coords); artemis->AddParam("coord_sys", sys); + artemis->AddParam("ndim", ndim); geometry::CoordParams cpars(pin.get()); artemis->AddParam("coord_params", cpars); @@ -202,6 +206,12 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { } } if (do_raytrace) packages.Add(RT::Initialize(pin.get(), units, constants)); + if (do_mhd) packages.Add(MHD::Initialize(pin.get(), units, constants, packages)); + + if (do_mhd && pin->GetOrAddBoolean("mhd", "monitor_divb", false)) { + artemis->PreStepDiagnosticsMesh = ArtemisUtils::PreStepDiagnosticsRemeshDivB; + artemis->PostStepDiagnosticsMesh = ArtemisUtils::PostStepDiagnosticsRemeshDivB; + } // Assign geometry-specific FillDerived functions if (do_gas || do_dust) { diff --git a/src/artemis.hpp b/src/artemis.hpp index 20e13e93..cc6f2e49 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -96,6 +96,17 @@ PAR_VAR(grav, phi); PAR_VAR(grav, rhs); } // namespace grav +namespace field { +namespace face { +PAR_VAR(field.face, B); +} // namespace face +namespace cell { +PAR_VAR(field.cell, B); +PAR_VAR(field.cell, energy); +PAR_VAR(field.cell, divB); +} // namespace cell +} // namespace field + namespace geom { PAR_VAR(geom, x1v); PAR_VAR(geom, x2v); @@ -112,9 +123,15 @@ PAR_VAR(geom, hx3f2); PAR_VAR(geom, hx1f3); PAR_VAR(geom, hx2f3); PAR_VAR(geom, hx3f3); +PAR_VAR(geom, hx1e1); +PAR_VAR(geom, hx2e2); +PAR_VAR(geom, hx3e3); PAR_VAR(geom, dx1); PAR_VAR(geom, dx2); PAR_VAR(geom, dx3); +PAR_VAR(geom, dl1); +PAR_VAR(geom, dl2); +PAR_VAR(geom, dl3); PAR_VAR(geom, vol); PAR_VAR(geom, ax1); PAR_VAR(geom, ax2); @@ -159,7 +176,7 @@ enum class Coordinates { }; // ...Riemann solvers -enum class RSolver { hllc_general, hlle, llf, hllc_gamma, null }; +enum class RSolver { hllc_general, hlle, llf, hllc_gamma, hlld, null }; // ... Upwinding (left vs right state) enum class Upwind { l, r, null }; // ...Reconstruction algorithms diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index a4af03db..75905119 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -75,6 +75,7 @@ ArtemisDriver::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in do_moment = artemis_pkg->template Param("do_moment"); do_coagulation = artemis_pkg->template Param("do_coagulation"); do_raytrace = artemis_pkg->template Param("do_raytrace"); + do_mhd = artemis_pkg->template Param("do_mhd"); // Update fluxes option--gas fields are needed for radiation temperature updates but for // rad-only test problems turn off advection @@ -273,9 +274,14 @@ TaskCollection ArtemisDriver::StepTasks() { diff_flx = vflx | tflx; } + TaskID edge_emf = none; + if (do_mhd) { + edge_emf = tl.AddTask(gas_flx, MHD::AssembleEdgeEMF, u0.get()); + } + // Communicate and set fluxes auto send_flx = - tl.AddTask(gas_flx | dust_flx | diff_flx, + tl.AddTask(gas_flx | dust_flx | diff_flx | edge_emf, parthenon::SendBoundBufs, u0); auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0); auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0); @@ -284,10 +290,16 @@ TaskCollection ArtemisDriver::StepTasks() { auto update = tl.AddTask(gas_flx | dust_flx | set_flx, ArtemisUtils::ApplyUpdate, u0.get(), u1.get(), g0, g1, bdt); + auto update_mhd = gas_flx | set_flx; + if (do_mhd) { + update_mhd = tl.AddTask(edge_emf | set_flx, ArtemisUtils::ApplyFaceUpdate, + u0.get(), u1.get(), g0, g1, bdt); + } // Apply "coordinate source terms" - TaskID gas_coord_src = update, dust_coord_src = update; - if (do_gas) gas_coord_src = tl.AddTask(update, Gas::FluxSource, u0.get(), bdt); + TaskID gas_coord_src = update | update_mhd, dust_coord_src = update; + if (do_gas) + gas_coord_src = tl.AddTask(update | update_mhd, Gas::FluxSource, u0.get(), bdt); if (do_dust) dust_coord_src = tl.AddTask(update, Dust::FluxSource, u0.get(), bdt); // Apply (gas) diffusion sources diff --git a/src/artemis_driver.hpp b/src/artemis_driver.hpp index 0ac79ce4..bf5fcef5 100644 --- a/src/artemis_driver.hpp +++ b/src/artemis_driver.hpp @@ -57,7 +57,7 @@ class ArtemisDriver : public EvolutionDriver { bool do_gravity, do_self_gravity, do_nbody, do_rotating_frame; bool do_orbital_advection; bool do_cooling, do_drag, do_viscosity, do_conduction, do_diffusion; - bool do_coagulation, do_raytrace; + bool do_coagulation, do_raytrace, do_mhd; bool update_fluxes; const bool is_restart; }; diff --git a/src/artemis_params.yaml b/src/artemis_params.yaml index b15ca028..9c01dbd0 100644 --- a/src/artemis_params.yaml +++ b/src/artemis_params.yaml @@ -39,7 +39,7 @@ artemis: _description: "How to provide unit conversions between code and physical units" base: _type: opt - _description: "Provide base unit conversions (length, time, mass, temperature)" + _description: "Provide base unit conversions (length, time, mass, temperature, current)" ppd: _type: opt _description: "AU, Year/(2 pi), solar mass units for protoplanetary disks" @@ -59,6 +59,10 @@ artemis: _type: "Real" _default: "1.0" _description: "Physical units value of temperature equal to 1 code unit" + current: + _type: "Real" + _default: "1.0" + _description: "Physical units value of electrical current equal to 1 code unit; used to derive magnetic-field and magnetic-permeability conversions" radial_spacing: _type: string _default: "uniform" @@ -70,6 +74,10 @@ physics: _type: "bool" _default: "true" _description: "Turn on the gas fluid" + mhd: + _type: "bool" + _default: "false" + _description: "Turn on ideal MHD and constrained-transport magnetic fields" dust: _type: "bool" _default: "false" diff --git a/src/derived/fill_derived.cpp b/src/derived/fill_derived.cpp index b65d13e3..a363dfa9 100644 --- a/src/derived/fill_derived.cpp +++ b/src/derived/fill_derived.cpp @@ -15,6 +15,7 @@ #include "fill_derived.hpp" #include "artemis.hpp" #include "geometry/geometry.hpp" +#include "mhd/mhd.hpp" #include "radiation/moments/moments.hpp" #include "utils/artemis_utils.hpp" #include "utils/eos/eos.hpp" @@ -23,6 +24,7 @@ using ArtemisUtils::EOS; using ArtemisUtils::VI; namespace ArtemisDerived { + //---------------------------------------------------------------------------------------- //! \fn TaskStatus ArtemisDerived::SetAuxillaryFields(MeshData *md) //! \brief Sets auxillary fields over IndexDomain::interior after an integration stage @@ -37,6 +39,7 @@ TaskStatus SetAuxillaryFields(MeshData *md) { // Return immediately if not evolving gas hydrodynamics auto &artemis_pkg = pm->packages.Get("artemis"); const bool do_gas = artemis_pkg->template Param("do_gas"); + const bool do_mhd = artemis_pkg->template Param("do_mhd"); if (!(do_gas)) return TaskStatus::complete; // Extract gas parameters @@ -47,7 +50,8 @@ TaskStatus SetAuxillaryFields(MeshData *md) { // Packing and indexing static auto desc = MakePackDescriptor(resolved_pkgs.get()); + gas::cons::internal_energy, field::cell::energy>( + resolved_pkgs.get()); auto vmesh = desc.GetPack(md); static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); @@ -78,7 +82,10 @@ TaskStatus SetAuxillaryFields(MeshData *md) { u_d = (dfloor)*u_d + (!dfloor) * dflr_gas; // Compute SIE via dual energy formalism and apply floor - Real sie = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx); + const Real emag = + (do_mhd && (n == 0)) ? vmesh(b, field::cell::energy(), k, j, i) : 0.0; + Real sie = + ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx, emag); const Real efloor = (sie > sieflr_gas); sie = (efloor)*sie + (!efloor) * sieflr_gas; @@ -105,6 +112,9 @@ void ConsToPrim(MeshData *md) { const bool do_gas = artemis_pkg->template Param("do_gas"); const bool do_dust = artemis_pkg->template Param("do_dust"); const bool do_rad = artemis_pkg->template Param("do_moment"); + const bool do_mhd = artemis_pkg->template Param("do_mhd"); + const Real mu0_code = + do_mhd ? pm->packages.Get("mhd")->template Param("mu0_code") : 1.0; // Extract gas parameters Real dflr_gas = Null(), sieflr_gas = Null(); @@ -134,15 +144,25 @@ void ConsToPrim(MeshData *md) { gas::cons::density, gas::cons::momentum, gas::cons::internal_energy, gas::prim::density, gas::prim::velocity, gas::prim::sie, dust::cons::density, dust::cons::momentum, dust::prim::density, dust::prim::velocity, rad::cons::energy, - rad::cons::flux, rad::prim::energy, rad::prim::flux>(resolved_pkgs.get()); + rad::cons::flux, rad::prim::energy, rad::prim::flux, field::face::B, field::cell::B, + field::cell::energy, field::cell::divB>(resolved_pkgs.get()); auto vmesh = desc.GetPack(md); - static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); + static auto desc_g = + MakePackDescriptor( + resolved_pkgs.get()); auto vg = desc_g.GetPack(md); const int nblocks = md->NumBlocks(); IndexRange ib = md->GetBoundsI(IndexDomain::interior); IndexRange jb = md->GetBoundsJ(IndexDomain::interior); IndexRange kb = md->GetBoundsK(IndexDomain::interior); + IndexRange ibe = md->GetBoundsI(IndexDomain::entire); + IndexRange jbe = md->GetBoundsJ(IndexDomain::entire); + IndexRange kbe = md->GetBoundsK(IndexDomain::entire); + + const int ndim = pm->ndim; + const int multid = ndim >= 2; + const int threed = ndim == 3; parthenon::par_for( DEFAULT_LOOP_PATTERN, "ConsToPrim", parthenon::DevExecSpace(), 0, @@ -156,21 +176,25 @@ void ConsToPrim(MeshData *md) { if (do_gas) { for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { // Set primitive density and apply floor - const Real u_d = vmesh(b, gas::cons::density(n), k, j, i); + const Real u_d = vmesh(b, TE::CC, gas::cons::density(n), k, j, i); const bool dfloor = (u_d > dflr_gas); - Real &w_d = vmesh(b, gas::prim::density(n), k, j, i); + Real &w_d = vmesh(b, TE::CC, gas::prim::density(n), k, j, i); w_d = (dfloor)*u_d + (!dfloor) * dflr_gas; // Set primitive velocity - Real &vel1 = vmesh(b, gas::prim::velocity(VI(n, 0)), k, j, i); - Real &vel2 = vmesh(b, gas::prim::velocity(VI(n, 1)), k, j, i); - Real &vel3 = vmesh(b, gas::prim::velocity(VI(n, 2)), k, j, i); - vel1 = vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) / (w_d * hx[0]); - vel2 = vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) / (w_d * hx[1]); - vel3 = vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) / (w_d * hx[2]); + Real &vel1 = vmesh(b, TE::CC, gas::prim::velocity(VI(n, 0)), k, j, i); + Real &vel2 = vmesh(b, TE::CC, gas::prim::velocity(VI(n, 1)), k, j, i); + Real &vel3 = vmesh(b, TE::CC, gas::prim::velocity(VI(n, 2)), k, j, i); + vel1 = + vmesh(b, TE::CC, gas::cons::momentum(VI(n, 0)), k, j, i) / (w_d * hx[0]); + vel2 = + vmesh(b, TE::CC, gas::cons::momentum(VI(n, 1)), k, j, i) / (w_d * hx[1]); + vel3 = + vmesh(b, TE::CC, gas::cons::momentum(VI(n, 2)), k, j, i) / (w_d * hx[2]); // Set primitive specific internal energy and apply floor - const Real w_s = vmesh(b, gas::cons::internal_energy(n), k, j, i) / w_d; + const Real w_s = + vmesh(b, TE::CC, gas::cons::internal_energy(n), k, j, i) / w_d; const bool siefloor = (w_s > sieflr_gas); vmesh(b, gas::prim::sie(n), k, j, i) = (siefloor)*w_s + (!siefloor) * sieflr_gas; @@ -180,18 +204,21 @@ void ConsToPrim(MeshData *md) { if (do_dust) { for (int n = 0; n < vmesh.GetSize(b, dust::prim::density()); ++n) { // Set primitive density - const Real u_d = vmesh(b, dust::cons::density(n), k, j, i); + const Real u_d = vmesh(b, TE::CC, dust::cons::density(n), k, j, i); const bool dfloor = (u_d > dflr_dust); - Real &w_d = vmesh(b, dust::prim::density(n), k, j, i); + Real &w_d = vmesh(b, TE::CC, dust::prim::density(n), k, j, i); w_d = (dfloor)*u_d + (u_d <= dflr_dust) * dflr_dust; // Set primitive velocity - Real &vel1 = vmesh(b, dust::prim::velocity(VI(n, 0)), k, j, i); - Real &vel2 = vmesh(b, dust::prim::velocity(VI(n, 1)), k, j, i); - Real &vel3 = vmesh(b, dust::prim::velocity(VI(n, 2)), k, j, i); - vel1 = vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) / (w_d * hx[0]); - vel2 = vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) / (w_d * hx[1]); - vel3 = vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) / (w_d * hx[2]); + Real &vel1 = vmesh(b, TE::CC, dust::prim::velocity(VI(n, 0)), k, j, i); + Real &vel2 = vmesh(b, TE::CC, dust::prim::velocity(VI(n, 1)), k, j, i); + Real &vel3 = vmesh(b, TE::CC, dust::prim::velocity(VI(n, 2)), k, j, i); + vel1 = + vmesh(b, TE::CC, dust::cons::momentum(VI(n, 0)), k, j, i) / (w_d * hx[0]); + vel2 = + vmesh(b, TE::CC, dust::cons::momentum(VI(n, 1)), k, j, i) / (w_d * hx[1]); + vel3 = + vmesh(b, TE::CC, dust::cons::momentum(VI(n, 2)), k, j, i) / (w_d * hx[2]); } } @@ -206,19 +233,105 @@ void ConsToPrim(MeshData *md) { // Set primitive radiation flux const Real cer = c * w_er; const std::array conv = {cer * hx[0], cer * hx[1], cer * hx[2]}; - const Real hfx1 = vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) / conv[0]; - const Real hfx2 = vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) / conv[1]; - const Real hfx3 = vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) / conv[2]; + const Real hfx1 = + vmesh(b, TE::CC, rad::cons::flux(VI(n, 0)), k, j, i) / conv[0]; + const Real hfx2 = + vmesh(b, TE::CC, rad::cons::flux(VI(n, 1)), k, j, i) / conv[1]; + const Real hfx3 = + vmesh(b, TE::CC, rad::cons::flux(VI(n, 2)), k, j, i) / conv[2]; const auto fx = Moments::NormalizeFlux(hfx1, hfx2, hfx3); - vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) = fx[0] * conv[0]; - vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) = fx[1] * conv[1]; - vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) = fx[2] * conv[2]; - vmesh(b, rad::prim::flux(VI(n, 0)), k, j, i) = fx[0]; - vmesh(b, rad::prim::flux(VI(n, 1)), k, j, i) = fx[1]; - vmesh(b, rad::prim::flux(VI(n, 2)), k, j, i) = fx[2]; + vmesh(b, TE::CC, rad::cons::flux(VI(n, 0)), k, j, i) = fx[0] * conv[0]; + vmesh(b, TE::CC, rad::cons::flux(VI(n, 1)), k, j, i) = fx[1] * conv[1]; + vmesh(b, TE::CC, rad::cons::flux(VI(n, 2)), k, j, i) = fx[2] * conv[2]; + vmesh(b, TE::CC, rad::prim::flux(VI(n, 0)), k, j, i) = fx[0]; + vmesh(b, TE::CC, rad::prim::flux(VI(n, 1)), k, j, i) = fx[1]; + vmesh(b, TE::CC, rad::prim::flux(VI(n, 2)), k, j, i) = fx[2]; } } + if (do_mhd) { + const Real vol = coords.GetVolume(vg, b, k, j, i); + const auto ax1 = coords.GetFaceAreaX1(vg, b, k, j, i); + const auto ax2 = coords.GetFaceAreaX2(vg, b, k, j, i); + const auto ax3 = coords.GetFaceAreaX3(vg, b, k, j, i); + const auto xv = coords.GetCellCenter(vg, b, k, j, i); + const auto &bnds = coords.GetBounds(); + vmesh(b, TE::CC, field::cell::divB(), k, j, i) = + ((ax1[1] * vmesh(b, TE::F1, field::face::B(), k, j, i + 1) - + ax1[0] * vmesh(b, TE::F1, field::face::B(), k, j, i)) + + (ax2[1] * vmesh(b, TE::F2, field::face::B(), k, j + multid, i) - + ax2[0] * vmesh(b, TE::F2, field::face::B(), k, j, i)) + + (ax3[1] * vmesh(b, TE::F3, field::face::B(), k + threed, j, i) - + ax3[0] * vmesh(b, TE::F3, field::face::B(), k, j, i))) / + vol; + vmesh(b, TE::CC, field::cell::B(0), k, j, i) = + ((bnds.x1[1] - xv[0]) * vmesh(b, TE::F1, field::face::B(), k, j, i) + + (xv[0] - bnds.x1[0]) * vmesh(b, TE::F1, field::face::B(), k, j, i + 1)) / + (bnds.x1[1] - bnds.x1[0]); + vmesh(b, TE::CC, field::cell::B(1), k, j, i) = + multid + ? (((bnds.x2[1] - xv[1]) * vmesh(b, TE::F2, field::face::B(), k, j, i) + + (xv[1] - bnds.x2[0]) * + vmesh(b, TE::F2, field::face::B(), k, j + multid, i)) / + (bnds.x2[1] - bnds.x2[0])) + : vmesh(b, TE::F2, field::face::B(), k, j, i); + vmesh(b, TE::CC, field::cell::B(2), k, j, i) = + threed + ? (((bnds.x3[1] - xv[2]) * vmesh(b, TE::F3, field::face::B(), k, j, i) + + (xv[2] - bnds.x3[0]) * + vmesh(b, TE::F3, field::face::B(), k + threed, j, i)) / + (bnds.x3[1] - bnds.x3[0])) + : vmesh(b, TE::F3, field::face::B(), k, j, i); + vmesh(b, TE::CC, field::cell::energy(), k, j, i) = MHD::MagneticEnergyDensity( + vmesh(b, TE::CC, field::cell::B(0), k, j, i), + vmesh(b, TE::CC, field::cell::B(1), k, j, i), + vmesh(b, TE::CC, field::cell::B(2), k, j, i), mu0_code); + } }); + + if (do_mhd) { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ConsToPrim::MHDGhost", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kbe.s, kbe.e, jbe.s, jbe.e, ibe.s, ibe.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const Real vol = coords.GetVolume(vg, b, k, j, i); + const auto ax1 = coords.GetFaceAreaX1(vg, b, k, j, i); + const auto ax2 = coords.GetFaceAreaX2(vg, b, k, j, i); + const auto ax3 = coords.GetFaceAreaX3(vg, b, k, j, i); + const auto xv = coords.GetCellCenter(vg, b, k, j, i); + const auto &bnds = coords.GetBounds(); + vmesh(b, TE::CC, field::cell::divB(), k, j, i) = + ((ax1[1] * vmesh(b, TE::F1, field::face::B(), k, j, i + 1) - + ax1[0] * vmesh(b, TE::F1, field::face::B(), k, j, i)) + + (ax2[1] * vmesh(b, TE::F2, field::face::B(), k, j + multid, i) - + ax2[0] * vmesh(b, TE::F2, field::face::B(), k, j, i)) + + (ax3[1] * vmesh(b, TE::F3, field::face::B(), k + threed, j, i) - + ax3[0] * vmesh(b, TE::F3, field::face::B(), k, j, i))) / + vol; + vmesh(b, TE::CC, field::cell::B(0), k, j, i) = + ((bnds.x1[1] - xv[0]) * vmesh(b, TE::F1, field::face::B(), k, j, i) + + (xv[0] - bnds.x1[0]) * vmesh(b, TE::F1, field::face::B(), k, j, i + 1)) / + (bnds.x1[1] - bnds.x1[0] + Fuzz()); + vmesh(b, TE::CC, field::cell::B(1), k, j, i) = + multid + ? (((bnds.x2[1] - xv[1]) * vmesh(b, TE::F2, field::face::B(), k, j, i) + + (xv[1] - bnds.x2[0]) * + vmesh(b, TE::F2, field::face::B(), k, j + multid, i)) / + (bnds.x2[1] - bnds.x2[0] + Fuzz())) + : vmesh(b, TE::F2, field::face::B(), k, j, i); + vmesh(b, TE::CC, field::cell::B(2), k, j, i) = + threed + ? (((bnds.x3[1] - xv[2]) * vmesh(b, TE::F3, field::face::B(), k, j, i) + + (xv[2] - bnds.x3[0]) * + vmesh(b, TE::F3, field::face::B(), k + threed, j, i)) / + (bnds.x3[1] - bnds.x3[0] + Fuzz())) + : vmesh(b, TE::F3, field::face::B(), k, j, i); + vmesh(b, TE::CC, field::cell::energy(), k, j, i) = MHD::MagneticEnergyDensity( + vmesh(b, TE::CC, field::cell::B(0), k, j, i), + vmesh(b, TE::CC, field::cell::B(1), k, j, i), + vmesh(b, TE::CC, field::cell::B(2), k, j, i), mu0_code); + }); + } } //---------------------------------------------------------------------------------------- @@ -228,6 +341,7 @@ template void PrimToCons(T *md) { PARTHENON_INSTRUMENT using parthenon::MakePackDescriptor; + using TE = parthenon::TopologicalElement; auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; @@ -235,7 +349,10 @@ void PrimToCons(T *md) { auto &artemis_pkg = pm->packages.Get("artemis"); const bool do_gas = artemis_pkg->template Param("do_gas"); const bool do_dust = artemis_pkg->template Param("do_dust"); + const bool do_mhd = artemis_pkg->template Param("do_mhd"); const bool do_rad = artemis_pkg->template Param("do_moment"); + const Real mu0_code = + do_mhd ? pm->packages.Get("mhd")->template Param("mu0_code") : 1.0; // Extract gas parameters Real dflr_gas = Null(); @@ -272,7 +389,8 @@ void PrimToCons(T *md) { gas::prim::bmod, gas::prim::temperature, dust::cons::density, dust::cons::momentum, dust::prim::density, dust::prim::velocity, rad::cons::energy, rad::cons::flux, rad::prim::energy, - rad::prim::flux, rad::prim::pressure>(resolved_pkgs.get()); + rad::prim::flux, rad::prim::pressure, field::face::B>( + resolved_pkgs.get()); auto vmesh = desc.GetPack(md); static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); @@ -280,6 +398,9 @@ void PrimToCons(T *md) { IndexRange ibe = md->GetBoundsI(IndexDomain::entire); IndexRange jbe = md->GetBoundsJ(IndexDomain::entire); IndexRange kbe = md->GetBoundsK(IndexDomain::entire); + const int ndim = md->GetMeshPointer()->ndim; + const int multid = ndim >= 2; + const int threed = ndim == 3; parthenon::par_for( DEFAULT_LOOP_PATTERN, "PrimToCons", parthenon::DevExecSpace(), 0, @@ -326,8 +447,20 @@ void PrimToCons(T *md) { // Sync conserved total energy const Real ke = 0.5 * w_d * (SQR(vel1) + SQR(vel2) + SQR(vel3)); + Real me = 0.0; + if (do_mhd && (n == 0)) { + const Real bx = 0.5 * (vmesh(b, TE::F1, field::face::B(), k, j, i) + + vmesh(b, TE::F1, field::face::B(), k, j, i + 1)); + const Real by = + 0.5 * (vmesh(b, TE::F2, field::face::B(), k, j, i) + + vmesh(b, TE::F2, field::face::B(), k, j + multid, i)); + const Real bz = + 0.5 * (vmesh(b, TE::F3, field::face::B(), k, j, i) + + vmesh(b, TE::F3, field::face::B(), k + threed, j, i)); + me = MHD::MagneticEnergyDensity(bx, by, bz, mu0_code); + } Real &u_e = vmesh(b, gas::cons::total_energy(n), k, j, i); - u_e = u_u + ke; + u_e = u_u + ke + me; } } @@ -392,6 +525,7 @@ template void PostInitialization(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_INSTRUMENT auto &md = pmb->meshblock_data.Get(); + MHD::SetCellCenteredMagneticFields, GEOM>(md.get()); PrimToCons, GEOM>(md.get()); } diff --git a/src/drag/drag.hpp b/src/drag/drag.hpp index 4358cdd8..1b86d8d4 100644 --- a/src/drag/drag.hpp +++ b/src/drag/drag.hpp @@ -178,6 +178,7 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt auto &artemis_pkg = pm->packages.Get("artemis"); const bool do_gas = artemis_pkg->template Param("do_gas"); const bool do_dust = artemis_pkg->template Param("do_dust"); + const bool do_mhd = artemis_pkg->template Param("do_mhd"); // Extract gas parameters Real de_switch = Null(); @@ -208,8 +209,8 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt // Packing and indexing static auto desc = MakePackDescriptor(resolved_pkgs.get()); + gas::cons::internal_energy, field::cell::energy, + dust::cons::momentum, dust::cons::density>(resolved_pkgs.get()); auto vmesh = desc.GetPack(md); static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); @@ -263,7 +264,10 @@ TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt dens = (dfloor)*dens + (!dfloor) * dflr_gas; // Compute SIE via dual energy formalism and apply floor - Real sieg = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx); + const Real emag = + (do_mhd && (n == 0)) ? vmesh(b, field::cell::energy(), k, j, i) : 0.0; + Real sieg = + ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx, emag); const Real efloor = (sieg > sieflr_gas); sieg = (efloor)*sieg + (!efloor) * sieflr_gas; @@ -355,6 +359,7 @@ TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real const int multi_d = (ndim >= 2); const int three_d = (ndim == 3); auto &artemis_pkg = pm->packages.Get("artemis"); + const bool do_mhd = artemis_pkg->template Param("do_mhd"); // Extract gas package and params auto &gas_pkg = pm->packages.Get("gas"); @@ -380,8 +385,8 @@ TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real // Packing and indexing static auto desc = MakePackDescriptor(resolved_pkgs.get()); + gas::cons::internal_energy, field::cell::energy, + dust::cons::momentum, dust::cons::density>(resolved_pkgs.get()); auto vmesh = desc.GetPack(md); static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); @@ -446,7 +451,9 @@ TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real dg = (dfloor)*dg + (!dfloor) * dflr_gas; // Compute SIE via dual energy formalism and apply floor - Real sieg = ArtemisUtils::DualEnergySIE(vmesh, b, 0, k, j, i, de_switch, hx); + const Real emag = do_mhd ? vmesh(b, field::cell::energy(), k, j, i) : 0.0; + Real sieg = + ArtemisUtils::DualEnergySIE(vmesh, b, 0, k, j, i, de_switch, hx, emag); const Real efloor = (sieg > sieflr_gas); sieg = (efloor)*sieg + (!efloor) * sieflr_gas; diff --git a/src/gas/cooling/beta_cooling.cpp b/src/gas/cooling/beta_cooling.cpp index b1aa6e36..27d21411 100644 --- a/src/gas/cooling/beta_cooling.cpp +++ b/src/gas/cooling/beta_cooling.cpp @@ -46,6 +46,7 @@ TaskStatus BetaCooling(MeshData *md, const Real time, const Real dt) { auto &resolved_pkgs = pm->resolved_packages; // Extract gas package and params + const bool do_mhd = pm->packages.Get("artemis")->template Param("do_mhd"); auto &gas_pkg = pm->packages.Get("gas"); const auto &eos_d = gas_pkg->template Param("eos_d"); const auto de_switch = gas_pkg->template Param("de_switch"); @@ -66,8 +67,8 @@ TaskStatus BetaCooling(MeshData *md, const Real time, const Real dt) { // Packing and indexing static auto desc = MakePackDescriptor( - resolved_pkgs.get()); + gas::cons::internal_energy, gas::cons::density, + field::cell::energy>(resolved_pkgs.get()); auto vmesh = desc.GetPack(md); static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); @@ -126,7 +127,10 @@ TaskStatus BetaCooling(MeshData *md, const Real time, const Real dt) { dens = (dfloor)*dens + (!dfloor) * dflr_gas; // Compute SIE via dual energy formalism and apply floor - Real sie = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx); + const Real emag = + (do_mhd && (n == 0)) ? vmesh(b, field::cell::energy(), k, j, i) : 0.0; + Real sie = + ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx, emag); const Real efloor = (sie > sieflr_gas); sie = (efloor)*sie + (!efloor) * sieflr_gas; diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index f643b472..59a7c40c 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -202,6 +202,12 @@ std::shared_ptr Initialize(ParameterInput *pin, riemann_solver = RSolver::hlle; } else if (riemann.compare("llf") == 0) { riemann_solver = RSolver::llf; + } else if (riemann.compare("hlld") == 0) { + PARTHENON_REQUIRE(eos_type == "ideal", + "The hlld Riemann solver requires the ideal eos."); + PARTHENON_REQUIRE(pin->GetOrAddBoolean("physics", "mhd", false), + "The hlld Riemann solver requires MHD to be enabled."); + riemann_solver = RSolver::hlld; } else { PARTHENON_FAIL("Riemann solver (gas) not recognized."); } @@ -544,9 +550,14 @@ Real EstimateTimestepMesh(MeshData *md) { auto ¶ms = gas_pkg->AllParams(); auto eos_d = params.template Get("eos_d"); + auto &artemis_pkg = pm->packages.Get("artemis"); + const bool do_mhd = artemis_pkg->template Param("do_mhd"); + const Real mu0_code = + do_mhd ? pm->packages.Get("mhd")->template Param("mu0_code") : 1.0; + static auto desc = MakePackDescriptor(resolved_pkgs.get()); + gas::prim::bmod, field::cell::B>(resolved_pkgs.get()); auto vmesh = desc.GetPack(md); static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); @@ -570,11 +581,23 @@ Real EstimateTimestepMesh(MeshData *md) { for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { const Real &dens = vmesh(b, gas::prim::density(n), k, j, i); const Real &bulk = vmesh(b, gas::prim::bmod(n), k, j, i); - const Real cs = std::sqrt(bulk / dens); + + // NOTE(AMD): This is an upper bound on the fast magnetosonic speed. + // It is much cheaper to compute since we would need to compute + // the fast magnetosonic speed at each face of the zone to do it more + // accurately. + Real b2 = 0.0; + if (do_mhd) { + b2 += SQR(vmesh(b, TE::CC, field::cell::B(0), k, j, i)) + + SQR(vmesh(b, TE::CC, field::cell::B(1), k, j, i)) + + SQR(vmesh(b, TE::CC, field::cell::B(2), k, j, i)); + } + const Real ss = std::sqrt((bulk + b2 / mu0_code) / dens); + Real denom = 0.0; for (int d = 0; d < ndim; d++) { denom += - (std::abs(vmesh(b, gas::prim::velocity(VI(n, d)), k, j, i)) + cs) / dx[d]; + (std::abs(vmesh(b, gas::prim::velocity(VI(n, d)), k, j, i)) + ss) / dx[d]; } ldt = std::min(ldt, 1.0 / denom); } @@ -630,12 +653,14 @@ TaskStatus CalculateFluxes(MeshData *md, const bool pcm) { static auto desc_prim = parthenon::MakePackDescriptor( + gas::prim::pressure, gas::prim::sie, gas::prim::bmod, + field::cell::B, field::cell::energy>( resolved_pkgs.get(), {}, {parthenon::PDOpt::WithFluxes}); static auto desc_flux = parthenon::MakePackDescriptor( - resolved_pkgs.get(), {}, {parthenon::PDOpt::WithFluxes}); + gas::cons::total_energy, gas::cons::internal_energy, + field::face::B>(resolved_pkgs.get(), {}, + {parthenon::PDOpt::WithFluxes}); static auto desc_face = parthenon::MakePackDescriptor(resolved_pkgs.get()); static auto desc_g = @@ -663,7 +688,8 @@ TaskStatus FluxSource(MeshData *md, const Real dt) { static auto desc_prim = parthenon::MakePackDescriptor(resolved_pkgs.get(), {}, + gas::prim::pressure, field::cell::B, + field::cell::energy>(resolved_pkgs.get(), {}, {parthenon::PDOpt::WithFluxes}); static auto desc_cons = parthenon::MakePackDescriptor( diff --git a/src/gas/params.yaml b/src/gas/params.yaml index b2a579e9..0bb1712a 100644 --- a/src/gas/params.yaml +++ b/src/gas/params.yaml @@ -76,6 +76,9 @@ gas: hlle: _type: opt _description: "HLLE" + hlld: + _type: opt + _description: "HLLD for ideal-MHD problems with the ideal EOS." hllc-general: _type: opt _description: "HLLC for a general EOS." diff --git a/src/geometry/axisymmetric.hpp b/src/geometry/axisymmetric.hpp index ed9d3b79..007b69b0 100644 --- a/src/geometry/axisymmetric.hpp +++ b/src/geometry/axisymmetric.hpp @@ -46,7 +46,7 @@ constexpr bool is_x1dep() { std::is_same_v || std::is_same_v || std::is_same_v) || std::is_same_v || std::is_same_v || - std::is_same_v; + std::is_same_v || std::is_same_v; } template constexpr bool is_x2dep() { diff --git a/src/geometry/cylindrical.hpp b/src/geometry/cylindrical.hpp index fdc763ec..65121542 100644 --- a/src/geometry/cylindrical.hpp +++ b/src/geometry/cylindrical.hpp @@ -43,7 +43,7 @@ constexpr bool is_x1dep() { std::is_same_v || std::is_same_v || std::is_same_v) || std::is_same_v || std::is_same_v || - std::is_same_v; + std::is_same_v || std::is_same_v; } template constexpr bool is_x2dep() { diff --git a/src/geometry/geometry.cpp b/src/geometry/geometry.cpp index bf896980..8693565f 100644 --- a/src/geometry/geometry.cpp +++ b/src/geometry/geometry.cpp @@ -50,6 +50,9 @@ void EnrollFields(StateDescriptor *pkg, CoordParams &cpars) { ADD_FIELD(geom::hx1f3); ADD_FIELD(geom::hx2f3); ADD_FIELD(geom::hx3f3); + ADD_FIELD(geom::hx1e1); + ADD_FIELD(geom::hx2e2); + ADD_FIELD(geom::hx3e3); ADD_FIELD(geom::dx1); ADD_FIELD(geom::dx2); ADD_FIELD(geom::dx3); @@ -121,8 +124,8 @@ void InitBlockGeom(MeshBlock *pmb, ParameterInput *pin) { geom::dh2dx1, geom::dh3dx1, geom::dh1dx2, geom::dh2dx2, geom::dh3dx2, geom::dh1dx3, geom::dh2dx3, geom::dh3dx3, geom::rfw1m, geom::rfw1p, geom::rfw2m, geom::rfw2p, geom::rfw3m, geom::rfw3p, geom::hx1f1, geom::hx1f2, geom::hx1f3, geom::hx2f1, - geom::hx2f2, geom::hx2f3, geom::hx3f1, geom::hx3f2, geom::hx3f3>( - (pm->resolved_packages).get()); + geom::hx2f2, geom::hx2f3, geom::hx3f1, geom::hx3f2, geom::hx3f3, geom::hx1e1, + geom::hx2e2, geom::hx3e3>((pm->resolved_packages).get()); auto vg = desc_g.GetPack(md.get()); IndexRange ib = md->GetBoundsI(IndexDomain::entire); IndexRange jb = md->GetBoundsJ(IndexDomain::entire); @@ -280,6 +283,79 @@ void InitBlockGeom(MeshBlock *pmb, ParameterInput *pin) { } } + { + auto xe = coords.EdgeCenX1(CellFace::lower, CellFace::lower); + Real &hx1e = + vg(b, geom::hx1e1())(coords.template index(k, j, i)); + Kokkos::atomic_store(&hx1e, coords.hx1(xe[0], xe[1], xe[2])); + if (x2end) { + xe = coords.EdgeCenX1(CellFace::upper, CellFace::lower); + Real &hx1e = + vg(b, geom::hx1e1())(coords.template index(k, j + 1, i)); + Kokkos::atomic_store(&hx1e, coords.hx1(xe[0], xe[1], xe[2])); + } + if (x3end) { + xe = coords.EdgeCenX1(CellFace::lower, CellFace::upper); + Real &hx1e = + vg(b, geom::hx1e1())(coords.template index(k + 1, j, i)); + Kokkos::atomic_store(&hx1e, coords.hx1(xe[0], xe[1], xe[2])); + } + if (x2end && x3end) { + xe = coords.EdgeCenX1(CellFace::upper, CellFace::upper); + Real &hx1e = vg(b, geom::hx1e1())( + coords.template index(k + 1, j + 1, i)); + Kokkos::atomic_store(&hx1e, coords.hx1(xe[0], xe[1], xe[2])); + } + } + { + auto xe = coords.EdgeCenX2(CellFace::lower, CellFace::lower); + Real &hx2e = + vg(b, geom::hx2e2())(coords.template index(k, j, i)); + Kokkos::atomic_store(&hx2e, coords.hx2(xe[0], xe[1], xe[2])); + if (x1end) { + xe = coords.EdgeCenX2(CellFace::upper, CellFace::lower); + Real &hx2e = + vg(b, geom::hx2e2())(coords.template index(k, j, i + 1)); + Kokkos::atomic_store(&hx2e, coords.hx2(xe[0], xe[1], xe[2])); + } + if (x3end) { + xe = coords.EdgeCenX2(CellFace::lower, CellFace::upper); + Real &hx2e = + vg(b, geom::hx2e2())(coords.template index(k + 1, j, i)); + Kokkos::atomic_store(&hx2e, coords.hx2(xe[0], xe[1], xe[2])); + } + if (x1end && x3end) { + xe = coords.EdgeCenX2(CellFace::upper, CellFace::upper); + Real &hx2e = vg(b, geom::hx2e2())( + coords.template index(k + 1, j, i + 1)); + Kokkos::atomic_store(&hx2e, coords.hx2(xe[0], xe[1], xe[2])); + } + } + { + auto xe = coords.EdgeCenX3(CellFace::lower, CellFace::lower); + Real &hx3e = + vg(b, geom::hx3e3())(coords.template index(k, j, i)); + Kokkos::atomic_store(&hx3e, coords.hx3(xe[0], xe[1], xe[2])); + if (x1end) { + xe = coords.EdgeCenX3(CellFace::upper, CellFace::lower); + Real &hx3e = + vg(b, geom::hx3e3())(coords.template index(k, j, i + 1)); + Kokkos::atomic_store(&hx3e, coords.hx3(xe[0], xe[1], xe[2])); + } + if (x2end) { + xe = coords.EdgeCenX3(CellFace::lower, CellFace::upper); + Real &hx3e = + vg(b, geom::hx3e3())(coords.template index(k, j + 1, i)); + Kokkos::atomic_store(&hx3e, coords.hx3(xe[0], xe[1], xe[2])); + } + if (x1end && x2end) { + xe = coords.EdgeCenX3(CellFace::upper, CellFace::upper); + Real &hx3e = vg(b, geom::hx3e3())( + coords.template index(k, j + 1, i + 1)); + Kokkos::atomic_store(&hx3e, coords.hx3(xe[0], xe[1], xe[2])); + } + } + // connection coeffs { auto dh = coords.GetConnX1(); diff --git a/src/geometry/geometry.hpp b/src/geometry/geometry.hpp index 6c3911b2..aee48f44 100644 --- a/src/geometry/geometry.hpp +++ b/src/geometry/geometry.hpp @@ -173,14 +173,20 @@ constexpr bool staggered_field() { (std::is_same_v && DIR == 1) || (std::is_same_v && DIR == 1) || (std::is_same_v && DIR == 1) || + (std::is_same_v && DIR == 1) || + (std::is_same_v && DIR == 1) || (std::is_same_v && DIR == 2) || (std::is_same_v && DIR == 2) || (std::is_same_v && DIR == 2) || (std::is_same_v && DIR == 2) || + (std::is_same_v && DIR == 2) || + (std::is_same_v && DIR == 2) || (std::is_same_v && DIR == 3) || (std::is_same_v && DIR == 3) || (std::is_same_v && DIR == 3) || - (std::is_same_v && DIR == 3); + (std::is_same_v && DIR == 3) || + (std::is_same_v && DIR == 3) || + (std::is_same_v && DIR == 3); } // NOTE(@amd) @@ -319,6 +325,22 @@ class CoordsBase { bnds.x3[static_cast(f)]}; } + KOKKOS_INLINE_FUNCTION std::array EdgeCenX1(const CellFace f2, + const CellFace f3) const { + return {static_cast(this)->x1v(), bnds.x2[static_cast(f2)], + bnds.x3[static_cast(f3)]}; + } + KOKKOS_INLINE_FUNCTION std::array EdgeCenX2(const CellFace f1, + const CellFace f3) const { + return {bnds.x1[static_cast(f1)], static_cast(this)->x2v(), + bnds.x3[static_cast(f3)]}; + } + KOKKOS_INLINE_FUNCTION std::array EdgeCenX3(const CellFace f1, + const CellFace f2) const { + return {bnds.x1[static_cast(f1)], bnds.x2[static_cast(f2)], + static_cast(this)->x3v()}; + } + template KOKKOS_INLINE_FUNCTION Real GetVolume(const V1 &vg, const int b, const int k, const int j, const int i) const { @@ -465,7 +487,9 @@ class CoordsBase { return static_cast(this)->template shape_(); } - KOKKOS_INLINE_FUNCTION BBox &GetBounds() const { + KOKKOS_INLINE_FUNCTION BBox &GetBounds() { return static_cast(this)->bnds; } + + KOKKOS_INLINE_FUNCTION const BBox &GetBounds() const { return static_cast(this)->bnds; } @@ -519,6 +543,58 @@ class CoordsBase { vg(b, geom::dx3())(index(k, j, i))}; } + template + KOKKOS_INLINE_FUNCTION Real GetEdgeScaleFactor() const { + PARTHENON_REQUIRE(DIR > 0 && DIR <= 3, "Invalid edge direction!"); + if constexpr (DIR == 1) { + const auto xe = EdgeCenX1(CellFace::lower, CellFace::lower); + return static_cast(this)->hx1(xe[0], xe[1], xe[2]); + } else if constexpr (DIR == 2) { + const auto xe = EdgeCenX2(CellFace::lower, CellFace::lower); + return static_cast(this)->hx2(xe[0], xe[1], xe[2]); + } + const auto xe = EdgeCenX3(CellFace::lower, CellFace::lower); + return static_cast(this)->hx3(xe[0], xe[1], xe[2]); + } + + template + KOKKOS_INLINE_FUNCTION Real GetEdgeScaleFactor(const V1 &vg, const int b, const int k, + const int j, const int i) const { + PARTHENON_REQUIRE(DIR > 0 && DIR <= 3, "Invalid edge direction!"); + if constexpr (CoordsTrait::value == Coordinates::cartesian) { + return 1.0; + } + if constexpr (DIR == 1) { + return vg(b, geom::hx1e1())(index(k, j, i)); + } else if constexpr (DIR == 2) { + return vg(b, geom::hx2e2())(index(k, j, i)); + } + return vg(b, geom::hx3e3())(index(k, j, i)); + } + template + KOKKOS_INLINE_FUNCTION Real GetEdgeLengthX1(const V1 &vg, const int b, const int k, + const int j, const int i) const { + return GetEdgeScaleFactor(vg, b, k, j, i) * (bnds.x1[1] - bnds.x1[0]); + } + KOKKOS_INLINE_FUNCTION Real GetEdgeLengthX1() const { + return GetEdgeScaleFactor() * (bnds.x1[1] - bnds.x1[0]); + } + template + KOKKOS_INLINE_FUNCTION Real GetEdgeLengthX2(const V1 &vg, const int b, const int k, + const int j, const int i) const { + return GetEdgeScaleFactor(vg, b, k, j, i) * (bnds.x2[1] - bnds.x2[0]); + } + KOKKOS_INLINE_FUNCTION Real GetEdgeLengthX2() const { + return GetEdgeScaleFactor() * (bnds.x2[1] - bnds.x2[0]); + } + template + KOKKOS_INLINE_FUNCTION Real GetEdgeLengthX3(const V1 &vg, const int b, const int k, + const int j, const int i) const { + return GetEdgeScaleFactor(vg, b, k, j, i) * (bnds.x3[1] - bnds.x3[0]); + } + KOKKOS_INLINE_FUNCTION Real GetEdgeLengthX3() const { + return GetEdgeScaleFactor() * (bnds.x3[1] - bnds.x3[0]); + } KOKKOS_INLINE_FUNCTION std::array GetCellCenter() const { // Get the cell centroid return {static_cast(this)->x1v(), static_cast(this)->x2v(), diff --git a/src/geometry/spherical.hpp b/src/geometry/spherical.hpp index b391a296..93daf27d 100644 --- a/src/geometry/spherical.hpp +++ b/src/geometry/spherical.hpp @@ -43,6 +43,7 @@ constexpr bool is_x1dep() { std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || + std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v; @@ -54,8 +55,9 @@ constexpr bool is_x2dep() { std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || - std::is_same_v || std::is_same_v || - std::is_same_v || std::is_same_v; + std::is_same_v || std::is_same_v || + std::is_same_v || std::is_same_v || + std::is_same_v; } template constexpr bool is_x3dep() { diff --git a/src/mhd/mhd.cpp b/src/mhd/mhd.cpp new file mode 100644 index 00000000..7b36b091 --- /dev/null +++ b/src/mhd/mhd.cpp @@ -0,0 +1,109 @@ +//======================================================================================== +// (C) (or copyright) 2026. 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. +//======================================================================================== + +// C++ headers +#include + +// Artemis includes +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "mhd.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/history.hpp" +#include "utils/refinement/amr_criteria.hpp" +#include "utils/units.hpp" + +using ArtemisUtils::VI; + +namespace MHD { +//---------------------------------------------------------------------------------------- +//! \fn StateDescriptor MHD::Initialize +//! \brief Adds intialization function for gas hydrodynamics package +std::shared_ptr Initialize(ParameterInput *pin, + ArtemisUtils::Units &units, + ArtemisUtils::Constants &constants, + Packages_t &packages) { + + auto mhd = std::make_shared("mhd"); + Params ¶ms = mhd->AllParams(); + + params.Add("mu0", constants.GetMu0Physical()); + params.Add("mu0_code", constants.GetMu0Code()); + + const int ndim = ProblemDimension(pin); + std::string sys = pin->GetOrAddString("artemis", "coordinates", "cartesian"); + Coordinates coords = geometry::CoordSelect(sys, ndim); + const bool log = + pin->GetOrAddString("artemis", "radial_spacing", "uniform") == "logarithmic"; + + Metadata m = Metadata({Metadata::Face, Metadata::Conserved, Metadata::Independent, + Metadata::WithFluxes, Metadata::FillGhost}); + ArtemisUtils::EnrollArtemisFaceRefinementOps(m, coords, log); + mhd->AddField(m); + + m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, Metadata::OneCopy, + Metadata::FillGhost, Metadata::WithFluxes}, + std::vector({3})); + mhd->AddField(m); + + m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, Metadata::OneCopy, + Metadata::FillGhost, Metadata::WithFluxes}); + mhd->AddField(m); + + m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}); + mhd->AddField(m); + return mhd; +} + +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus ArtemisUtils::AssembleEdgeEMF +//! \brief Runtime dispatch for geometry-aware edge EMF assembly. +TaskStatus AssembleEdgeEMF(MeshData *md) { + PARTHENON_INSTRUMENT + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + const auto &artemis_pkg = pm->packages.Get("artemis"); + const auto do_mhd = artemis_pkg->template Param("do_mhd"); + if (!do_mhd) return TaskStatus::complete; + + static auto desc = MakePackDescriptor( + resolved_pkgs.get(), {}, {parthenon::PDOpt::WithFluxes}); + static auto desc_g = + MakePackDescriptor(resolved_pkgs.get()); + const auto v = desc.GetPack(md); + const auto vg = desc_g.GetPack(md); + + const auto sys = artemis_pkg->template Param("coords"); + const auto &cpars = artemis_pkg->template Param("coord_params"); + typedef Coordinates G; + if (sys == G::cartesian) { + return AssembleEdgeEMFImpl(md, v, vg, cpars); + } else if (sys == G::spherical3D) { + return AssembleEdgeEMFImpl(md, v, vg, cpars); + } else if (sys == G::spherical1D) { + return AssembleEdgeEMFImpl(md, v, vg, cpars); + } else if (sys == G::spherical2D) { + return AssembleEdgeEMFImpl(md, v, vg, cpars); + } else if (sys == G::cylindrical) { + return AssembleEdgeEMFImpl(md, v, vg, cpars); + } else if (sys == G::axisymmetric) { + return AssembleEdgeEMFImpl(md, v, vg, cpars); + } else { + PARTHENON_FAIL("Coordinate type not recognized!"); + } +} + +} // namespace MHD \ No newline at end of file diff --git a/src/mhd/mhd.hpp b/src/mhd/mhd.hpp new file mode 100644 index 00000000..c5d5b91e --- /dev/null +++ b/src/mhd/mhd.hpp @@ -0,0 +1,268 @@ +//======================================================================================== +// (C) (or copyright) 2026. 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. +//======================================================================================== +#ifndef MHD_MHD_HPP_ +#define MHD_MHD_HPP_ + +#include "artemis.hpp" +#include "utils/units.hpp" + +namespace MHD { + +KOKKOS_FORCEINLINE_FUNCTION Real MagneticEnergyDensity(const Real bx, const Real by, + const Real bz, const Real mu0) { + return 0.5 * (SQR(bx) + SQR(by) + SQR(bz)) / mu0; +} + +KOKKOS_FORCEINLINE_FUNCTION Real FastMagnetosonicSpeed(const Real bulk, + const Real density, const Real b2, + const Real bx, const Real mu0) { + const Real wave_sum = (bulk + b2 / mu0) / density; + const Real discriminant = + std::max(0.0, SQR(wave_sum) - 4.0 * bulk * SQR(bx) / (mu0 * SQR(density))); + return std::sqrt(0.5 * (wave_sum + std::sqrt(discriminant))); +} + +KOKKOS_FORCEINLINE_FUNCTION Real FastMagnetosonicSpeed(const Real bulk, + const Real density, const Real bx, + const Real by, const Real bz, + const Real mu0) { + // NOTE(AMD): bx is the component normal to the interface + return FastMagnetosonicSpeed(bulk, density, SQR(bx) + SQR(by) + SQR(bz), bx, mu0); +} + +std::shared_ptr Initialize(ParameterInput *pin, + ArtemisUtils::Units &units, + ArtemisUtils::Constants &constants, + Packages_t &packages); +TaskStatus AssembleEdgeEMF(MeshData *md); + +template +void SetCellCenteredMagneticFields(T *md) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + + auto &artemis_pkg = pm->packages.Get("artemis"); + const bool do_mhd = artemis_pkg->template Param("do_mhd"); + if (!do_mhd) return; + + const auto &cpars = artemis_pkg->template Param("coord_params"); + + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + if (vmesh.GetMaxNumberOfVars() == 0) return; + + static auto desc_g = + MakePackDescriptor(resolved_pkgs.get()); + auto vg = desc_g.GetPack(md); + + const int ndim = md->GetMeshPointer()->ndim; + const int multid = ndim >= 2; + const int threed = ndim == 3; + IndexRange ibe = md->GetBoundsI(IndexDomain::entire); + IndexRange jbe = md->GetBoundsJ(IndexDomain::entire); + IndexRange kbe = md->GetBoundsK(IndexDomain::entire); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "SetCellCenteredMagneticFields", parthenon::DevExecSpace(), 0, + vmesh.GetNBlocks() - 1, kbe.s, kbe.e, jbe.s, jbe.e, ibe.s, ibe.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + + const auto xv = coords.GetCellCenter(vg, b, k, j, i); + const auto &bnds = coords.GetBounds(); + + vmesh(b, TE::CC, field::cell::B(0), k, j, i) = + ((bnds.x1[1] - xv[0]) * vmesh(b, TE::F1, field::face::B(), k, j, i) + + (xv[0] - bnds.x1[0]) * vmesh(b, TE::F1, field::face::B(), k, j, i + 1)) / + (bnds.x1[1] - bnds.x1[0]); + vmesh(b, TE::CC, field::cell::B(1), k, j, i) = + multid + ? (((bnds.x2[1] - xv[1]) * vmesh(b, TE::F2, field::face::B(), k, j, i) + + (xv[1] - bnds.x2[0]) * + vmesh(b, TE::F2, field::face::B(), k, j + multid, i)) / + (bnds.x2[1] - bnds.x2[0])) + : vmesh(b, TE::F2, field::face::B(), k, j, i); + vmesh(b, TE::CC, field::cell::B(2), k, j, i) = + threed + ? (((bnds.x3[1] - xv[2]) * vmesh(b, TE::F3, field::face::B(), k, j, i) + + (xv[2] - bnds.x3[0]) * + vmesh(b, TE::F3, field::face::B(), k + threed, j, i)) / + (bnds.x3[1] - bnds.x3[0])) + : vmesh(b, TE::F3, field::face::B(), k, j, i); + }); +} + +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus ArtemisUtils::AssembleEdgeEMFImpl +//! \brief Assemble unique edge EMFs from raw face induction fluxes on field::cell::B. +template +inline TaskStatus AssembleEdgeEMFImpl(MeshData *md, PACK v, GEO vg, + const geometry::CoordParams &cpars) { + PARTHENON_INSTRUMENT + auto pm = md->GetParentPointer(); + + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + const bool multi_d = (pm->ndim > 1); + const bool three_d = (pm->ndim > 2); + + // + // Fx(By) = -Ez and Fy(Bx) = Ez + // + // E_i-1j--Fy_i-1j+1--E_ij+1---Fy_ij+1--E_i+1j+1 + // | | | + // | | | + // Fx_i-1j Fx_ij Fx_i+1j + // | | | + // | | | + // E_i-1j---Fy_i-1j---E_ij------Fy_ij---E_i+1j + // | | | + // | | | + // Fx_i-1j-1 Fx_ij-1 Fx_i+1j-1 + // | | | + // | | | + // E_i-1j-1--Fy_i-1j--E_ij-1-----Fy_ij---E_i+1j-1 + // + if (multi_d) { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "AssembleEdgeEMF::E3", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e + 1, ib.s, ib.e + 1, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vg.GetCoordinates(b), k, j, i); + const auto hx1 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const auto hx1_jm = + coords.template GetScaleFactorsFace(vg, b, k, j - 1, i); + const auto hx2 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const auto hx2_im = + coords.template GetScaleFactorsFace(vg, b, k, j, i - 1); + const Real h3e = coords.template GetEdgeScaleFactor(vg, b, k, j, i); + Real &emf = v.flux(b, TE::E3, field::face::B(), k, j, i); + emf = h3e * 0.25 * + (-v.flux(b, X1DIR, field::cell::B(1), k, j, i) / hx1[1] - + v.flux(b, X1DIR, field::cell::B(1), k, j - 1, i) / hx1_jm[1] + + v.flux(b, X2DIR, field::cell::B(0), k, j, i) / hx2[0] + + v.flux(b, X2DIR, field::cell::B(0), k, j, i - 1) / hx2_im[0]); + }); + } else { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "AssembleEdgeEMF::E3", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e + 1, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vg.GetCoordinates(b), k, j, i); + const auto hx1 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const Real h3e = coords.template GetEdgeScaleFactor(vg, b, k, j, i); + Real &emf = v.flux(b, TE::E3, field::face::B(), k, j, i); + emf = -h3e * v.flux(b, X1DIR, field::cell::B(1), k, j, i) / hx1[1]; + }); + } + + if (three_d) { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "AssembleEdgeEMF::E2", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e + 1, jb.s, jb.e, ib.s, ib.e + 1, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vg.GetCoordinates(b), k, j, i); + const auto hx1 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const auto hx1_km = + coords.template GetScaleFactorsFace(vg, b, k - 1, j, i); + const auto hx3 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const auto hx3_im = + coords.template GetScaleFactorsFace(vg, b, k, j, i - 1); + const Real h2e = coords.template GetEdgeScaleFactor(vg, b, k, j, i); + Real &emf = v.flux(b, TE::E2, field::face::B(), k, j, i); + emf = h2e * 0.25 * + (v.flux(b, X1DIR, field::cell::B(2), k, j, i) / hx1[2] + + v.flux(b, X1DIR, field::cell::B(2), k - 1, j, i) / hx1_km[2] - + v.flux(b, X3DIR, field::cell::B(0), k, j, i) / hx3[0] - + v.flux(b, X3DIR, field::cell::B(0), k, j, i - 1) / hx3_im[0]); + }); + + if (multi_d) { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "AssembleEdgeEMF::E1", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e + 1, jb.s, jb.e + 1, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vg.GetCoordinates(b), k, j, i); + const auto hx2 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const auto hx2_km = + coords.template GetScaleFactorsFace(vg, b, k - 1, j, i); + const auto hx3 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const auto hx3_jm = + coords.template GetScaleFactorsFace(vg, b, k, j - 1, i); + const Real h1e = coords.template GetEdgeScaleFactor(vg, b, k, j, i); + Real &emf = v.flux(b, TE::E1, field::face::B(), k, j, i); + emf = h1e * 0.25 * + (-v.flux(b, X2DIR, field::cell::B(2), k, j, i) / hx2[2] - + v.flux(b, X2DIR, field::cell::B(2), k - 1, j, i) / hx2_km[2] + + v.flux(b, X3DIR, field::cell::B(1), k, j, i) / hx3[1] + + v.flux(b, X3DIR, field::cell::B(1), k, j - 1, i) / hx3_jm[1]); + }); + } + } else { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "AssembleEdgeEMF::E2", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e + 1, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vg.GetCoordinates(b), k, j, i); + const auto hx1 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const Real h2e = coords.template GetEdgeScaleFactor(vg, b, k, j, i); + Real &emf = v.flux(b, TE::E2, field::face::B(), k, j, i); + emf = h2e * v.flux(b, X1DIR, field::cell::B(2), k, j, i) / hx1[2]; + }); + + if (multi_d) { + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "AssembleEdgeEMF::E1", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e + 1, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vg.GetCoordinates(b), k, j, i); + const auto hx2 = coords.template GetScaleFactorsFace(vg, b, k, j, i); + const Real h1e = coords.template GetEdgeScaleFactor(vg, b, k, j, i); + Real &emf = v.flux(b, TE::E1, field::face::B(), k, j, i); + emf = -h1e * v.flux(b, X2DIR, field::cell::B(2), k, j, i) / hx2[2]; + }); + } + } + + return TaskStatus::complete; +} + +//---------------------------------------------------------------------------------------- +//! \fn void ArtemisUtils::ScaleMHDFlux +//! \brief Scales raw induction fluxes by scale factors associated with relevant coord sys +template +KOKKOS_INLINE_FUNCTION void ScaleMHDFlux(parthenon::team_mbr_t const &member, + const geometry::CoordParams &cpars, const int b, + const int k, const int j, const int il, + const int iu, const V4 &vg, const V3 &p) { + if constexpr (G == Coordinates::cartesian) return; + PARTHENON_REQUIRE(DIR > 0 && DIR <= 3, "Invalid flux direction!"); + + parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, il, iu, [&](const int i) { + geometry::Coords coords(cpars, p.GetCoordinates(b), k, j, i); + const auto &hx = coords.template GetScaleFactorsFace(vg, b, k, j, i); + p.flux(b, DIR, field::cell::B(0), k, j, i) *= hx[0]; + p.flux(b, DIR, field::cell::B(1), k, j, i) *= hx[1]; + p.flux(b, DIR, field::cell::B(2), k, j, i) *= hx[2]; + }); + + return; +} + +} // namespace MHD + +#endif // MHD_MHD_HPP_ diff --git a/src/pgen/blast.hpp b/src/pgen/blast.hpp index bba557c8..edd3b722 100644 --- a/src/pgen/blast.hpp +++ b/src/pgen/blast.hpp @@ -129,6 +129,7 @@ template inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_INSTRUMENT using parthenon::MakePackDescriptor; + using TE = parthenon::TopologicalElement; // Extract blast parameters blast_params.rinit = pin->GetOrAddReal("problem", "radius", 1.0); @@ -139,6 +140,11 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { blast_params.x0[0] = pin->GetOrAddReal("problem", "x1", 0.0); blast_params.x0[1] = pin->GetOrAddReal("problem", "x2", 0.0); blast_params.x0[2] = pin->GetOrAddReal("problem", "x3", 0.0); + + const Real bx1 = pin->GetOrAddReal("problem", "bx1", 0.0); + const Real bx2 = pin->GetOrAddReal("problem", "bx2", 0.0); + const Real bx3 = pin->GetOrAddReal("problem", "bx3", 0.0); + blast_params.samples = pin->GetOrAddInteger("problem", "samples", -1); std::string sym = pin->GetOrAddString("problem", "symmetry", "spherical"); if (sym == "spherical") { @@ -152,6 +158,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // Extract parameters from packages auto artemis_pkg = pmb->packages.Get("artemis"); const bool do_dust = artemis_pkg->Param("do_dust"); + const bool do_mhd = artemis_pkg->Param("do_mhd"); // TODO(PDM): Replace the below with a call to singularity-eos auto gas_pkg = pmb->packages.Get("gas"); const auto &eos = gas_pkg->Param("eos_d"); @@ -163,7 +170,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } static auto desc = MakePackDescriptor( + dust::prim::density, dust::prim::velocity, field::face::B>( (pmb->resolved_packages).get()); auto v = desc.GetPack(md.get()); static auto desc_g = MakePackDescriptor( @@ -176,6 +183,10 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto pars = blast_params; const auto &cpars = artemis_pkg->template Param("coord_params"); + const int ndim = pmb->pmy_mesh->ndim; + const int multid = ndim >= 2; + const int threed = ndim == 3; + // setup uniform ambient medium with spherical over-pressured region pmb->par_for( "blast", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -214,10 +225,8 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real vol = (pars.samples > 0) ? compute_overlap_cyl(coords.bnds, pars.rinit, pars.samples) - : ((SQR(xcart[0]) + SQR(xcart[1]) + SQR(xcart[2]) < - pars.rinit * pars.rinit) - ? total_vol - : 0.0); + : ((SQR(xcart[0]) + SQR(xcart[1]) < pars.rinit * pars.rinit) ? total_vol + : 0.0); internal_energy = e0 * (1.0 - vol / total_vol) + pars.internal_energy * vol / total_vol / (M_PI * pars.rinit * pars.rinit); @@ -231,11 +240,23 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real vx3 = 0.0; // compute cell-centered conserved variables - v(0, gas::prim::density(), k, j, i) = den; - v(0, gas::prim::velocity(0), k, j, i) = vx1; - v(0, gas::prim::velocity(1), k, j, i) = vx2; - v(0, gas::prim::velocity(2), k, j, i) = vx3; - v(0, gas::prim::sie(), k, j, i) = internal_energy / den; + v(0, TE::CC, gas::prim::density(), k, j, i) = den; + v(0, TE::CC, gas::prim::velocity(0), k, j, i) = vx1; + v(0, TE::CC, gas::prim::velocity(1), k, j, i) = vx2; + v(0, TE::CC, gas::prim::velocity(2), k, j, i) = vx3; + v(0, TE::CC, gas::prim::sie(), k, j, i) = internal_energy / den; + if (do_mhd) { + v(0, TE::F1, field::face::B(), k, j, i) = bx1; + if (i == ib.e) v(0, TE::F1, field::face::B(), k, j, ib.e + 1) = bx1; + if (multid || bx2 != 0.0) { + v(0, TE::F2, field::face::B(), k, j, i) = bx2; + if (j == jb.e) v(0, TE::F2, field::face::B(), k, jb.e + multid, i) = bx2; + } + if (threed || bx3 != 0.0) { + v(0, TE::F3, field::face::B(), k, j, i) = bx3; + if (k == kb.e) v(0, TE::F3, field::face::B(), kb.e + threed, j, i) = bx3; + } + } }); } diff --git a/src/pgen/orszag_tang.hpp b/src/pgen/orszag_tang.hpp new file mode 100644 index 00000000..703c8396 --- /dev/null +++ b/src/pgen/orszag_tang.hpp @@ -0,0 +1,153 @@ +//======================================================================================== +// (C) (or copyright) 2023-2026. 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. +//======================================================================================== +#ifndef PGEN_ORSZAG_TANG_HPP_ +#define PGEN_ORSZAG_TANG_HPP_ +//! \file orszag_tang.hpp +//! \brief Orszag-Tang vortex problem generator. + +// C++ headers +#include + +// Artemis headers +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/eos/eos.hpp" + +using ArtemisUtils::EOS; + +namespace { + +struct OrszagTangParams { + Real rho0; + Real p0; + Real v0; + Real b0; + Real x1min; + Real x1max; + Real x2min; + Real x2max; +}; + +} // namespace + +namespace orszag_tang { + +static OrszagTangParams otv; + +//---------------------------------------------------------------------------------------- +//! \fn void ProblemGenerator::OrszagTang() +//! \brief Sets initial conditions for the Orszag-Tang vortex. +template +inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + using TE = parthenon::TopologicalElement; + + auto artemis_pkg = pmb->packages.Get("artemis"); + const bool do_gas = artemis_pkg->Param("do_gas"); + const bool do_dust = artemis_pkg->Param("do_dust"); + const bool do_mhd = artemis_pkg->Param("do_mhd"); + const int ndim = pmb->pmy_mesh->ndim; + + PARTHENON_REQUIRE(GEOM == Coordinates::cartesian, + "orszag_tang pgen requires Cartesian geometry!"); + PARTHENON_REQUIRE(ndim >= 2, + "orszag_tang pgen requires at least 2 spatial dimensions!"); + PARTHENON_REQUIRE(do_gas, "orszag_tang pgen requires gas hydrodynamics!"); + PARTHENON_REQUIRE(!(do_dust), "orszag_tang pgen does not support dust!"); + PARTHENON_REQUIRE(do_mhd, "orszag_tang pgen requires MHD!"); + + auto gas_pkg = pmb->packages.Get("gas"); + PARTHENON_REQUIRE(gas_pkg->Param("nspecies") == 1, + "orszag_tang pgen requires a single gas species."); + const auto &eos = gas_pkg->Param("eos_h"); + + otv.rho0 = pin->GetOrAddReal("problem", "rho0", 25.0 / (36.0 * M_PI)); + otv.p0 = pin->GetOrAddReal("problem", "p0", 5.0 / (12.0 * M_PI)); + otv.v0 = pin->GetOrAddReal("problem", "v0", 1.0); + otv.b0 = pin->GetOrAddReal("problem", "b0", 1.0 / std::sqrt(4.0 * M_PI)); + otv.x1min = pmb->pmy_mesh->mesh_size.xmin(X1DIR); + otv.x1max = pmb->pmy_mesh->mesh_size.xmax(X1DIR); + otv.x2min = pmb->pmy_mesh->mesh_size.xmin(X2DIR); + otv.x2max = pmb->pmy_mesh->mesh_size.xmax(X2DIR); + + 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); + IndexRange ib1 = pmb->cellbounds.GetBoundsI(IndexDomain::entire, TE::F1); + IndexRange jb1 = pmb->cellbounds.GetBoundsJ(IndexDomain::entire, TE::F1); + IndexRange kb1 = pmb->cellbounds.GetBoundsK(IndexDomain::entire, TE::F1); + IndexRange ib2 = pmb->cellbounds.GetBoundsI(IndexDomain::entire, TE::F2); + IndexRange jb2 = pmb->cellbounds.GetBoundsJ(IndexDomain::entire, TE::F2); + IndexRange kb2 = pmb->cellbounds.GetBoundsK(IndexDomain::entire, TE::F2); + IndexRange ib3 = pmb->cellbounds.GetBoundsI(IndexDomain::entire, TE::F3); + IndexRange jb3 = pmb->cellbounds.GetBoundsJ(IndexDomain::entire, TE::F3); + IndexRange kb3 = pmb->cellbounds.GetBoundsK(IndexDomain::entire, TE::F3); + auto &pco = pmb->coords; + auto pars = otv; + + const Real lx = pars.x1max - pars.x1min; + const Real ly = pars.x2max - pars.x2min; + const Real sie0 = ArtemisUtils::EofPR(eos, pars.p0, pars.rho0); + const bool three_d = (ndim > 2); + + pmb->par_for( + "orszag_tang_cc", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + const Real x = (pco.template Xc(i) - pars.x1min) / lx; + const Real y = (pco.template Xc(j) - pars.x2min) / ly; + v(0, TE::CC, gas::prim::density(), k, j, i) = pars.rho0; + v(0, TE::CC, gas::prim::velocity(0), k, j, i) = + -pars.v0 * std::sin(2.0 * M_PI * y); + v(0, TE::CC, gas::prim::velocity(1), k, j, i) = + pars.v0 * std::sin(2.0 * M_PI * x); + v(0, TE::CC, gas::prim::velocity(2), k, j, i) = 0.0; + v(0, TE::CC, gas::prim::sie(), k, j, i) = sie0; + }); + + pmb->par_for( + "orszag_tang_b1", kb1.s, kb1.e, jb1.s, jb1.e, ib1.s, ib1.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + const Real y = (pco.template Xc(j) - pars.x2min) / ly; + v(0, TE::F1, field::face::B(), k, j, i) = -pars.b0 * std::sin(2.0 * M_PI * y); + }); + + pmb->par_for( + "orszag_tang_b2", kb2.s, kb2.e, jb2.s, jb2.e, ib2.s, ib2.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + const Real x = (pco.template Xc(i) - pars.x1min) / lx; + v(0, TE::F2, field::face::B(), k, j, i) = pars.b0 * std::sin(4.0 * M_PI * x); + }); + + if (three_d) { + pmb->par_for( + "orszag_tang_b3", kb3.s, kb3.e, jb3.s, jb3.e, ib3.s, ib3.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + v(0, TE::F3, field::face::B(), k, j, i) = 0.0; + }); + } +} + +} // namespace orszag_tang + +#endif // PGEN_ORSZAG_TANG_HPP_ diff --git a/src/pgen/params.yaml b/src/pgen/params.yaml index 8b919e2d..0783cae1 100644 --- a/src/pgen/params.yaml +++ b/src/pgen/params.yaml @@ -219,10 +219,25 @@ problem(linear_wave): _type: Real _default: 1.0 _description: "Sets new time limit to the specified number of periods" - along_x(1,2,3): - _type: bool - _default: "false" - _description: "Wavevctor along x(1,2,3) axis" + +problem(orszag_tang): + _description: " block parameters when artemis/problem=orszag_tang" + rho0: + _type: Real + _default: 25.0 / (36.0 * pi) + _description: "Uniform background density" + p0: + _type: Real + _default: 5.0 / (12.0 * pi) + _description: "Uniform background gas pressure" + v0: + _type: Real + _default: 1.0 + _description: "Velocity amplitude for the vortex pattern" + b0: + _type: Real + _default: 1.0 / sqrt(4.0 * pi) + _description: "Magnetic-field amplitude for the standard Orszag-Tang pattern" problem(shock): _description: " block parameters when artemis/problem=shock" @@ -237,7 +252,11 @@ problem(shock): tl: _type: Real _default: 0.6 - _description: "Temperature (left)." + _description: "Temperature (left). Used when pl < 0." + pl: + _type: Real + _default: -1.0 + _description: "Pressure (left). If >= 0, overrides tl." rhor: _type: Real _default: 2.285714 @@ -249,7 +268,31 @@ problem(shock): tr: _type: Real _default: 1.246875 - _description: "Temperature (right)." + _description: "Temperature (right). Used when pr < 0." + pr: + _type: Real + _default: -1.0 + _description: "Pressure (right). If >= 0, overrides tr." + bx: + _type: Real + _default: 0.0 + _description: "Uniform normal magnetic field Bx." + byl: + _type: Real + _default: 0.0 + _description: "Left transverse magnetic field By." + byr: + _type: Real + _default: 0.0 + _description: "Right transverse magnetic field By." + bzl: + _type: Real + _default: 0.0 + _description: "Left transverse magnetic field Bz." + bzr: + _type: Real + _default: 0.0 + _description: "Right transverse magnetic field Bz." xdisc: _type: Real _default: 0.0005 diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 9ae84c8f..68caac8d 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -32,6 +32,7 @@ #include "kh.hpp" #include "linear_wave.hpp" #include "lw.hpp" +#include "orszag_tang.hpp" #include "polytrope.hpp" #include "rt.hpp" #include "shock.hpp" @@ -70,6 +71,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { linear_wave::ProblemGenerator(pmb, pin); } else if (name == "lw") { lw::ProblemGenerator(pmb, pin); + } else if (name == "orszag_tang") { + orszag_tang::ProblemGenerator(pmb, pin); } else if (name == "polytrope") { polytrope::ProblemGenerator(pmb, pin); } else if (name == "kh") { diff --git a/src/pgen/shock.hpp b/src/pgen/shock.hpp index 9be64caf..b76d0f90 100644 --- a/src/pgen/shock.hpp +++ b/src/pgen/shock.hpp @@ -40,8 +40,11 @@ using ArtemisUtils::EOS; namespace shock { struct ShockParams { - Real rhol, vxl, tl; - Real rhor, vxr, tr; + Real rhol, vxl, tl, pl; + Real rhor, vxr, tr, pr; + Real bx; + Real byl, byr; + Real bzl, bzr; Real xdisc; }; @@ -56,9 +59,16 @@ inline void InitShockParams(MeshBlock *pmb, ParameterInput *pin) { shock_params.rhol = pin->GetOrAddReal("problem", "rhol", 5.69); shock_params.vxl = pin->GetOrAddReal("problem", "vxl", 5.19e7); shock_params.tl = pin->GetOrAddReal("problem", "tl", 2.18e6); + shock_params.pl = pin->GetOrAddReal("problem", "pl", Null()); shock_params.rhor = pin->GetOrAddReal("problem", "rhor", 17.1); shock_params.vxr = pin->GetOrAddReal("problem", "vxr", 1.73e7); shock_params.tr = pin->GetOrAddReal("problem", "tr", 7.98e6); + shock_params.pr = pin->GetOrAddReal("problem", "pr", Null()); + shock_params.bx = pin->GetOrAddReal("problem", "bx", 0.0); + shock_params.byl = pin->GetOrAddReal("problem", "byl", 0.0); + shock_params.byr = pin->GetOrAddReal("problem", "byr", 0.0); + shock_params.bzl = pin->GetOrAddReal("problem", "bzl", 0.0); + shock_params.bzr = pin->GetOrAddReal("problem", "bzr", 0.0); shock_params.xdisc = pin->GetOrAddReal("problem", "xdisc", 0.0005); params.Add("shock_params", shock_params); } @@ -71,6 +81,7 @@ template inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_INSTRUMENT using parthenon::MakePackDescriptor; + using TE = parthenon::TopologicalElement; // Extract parameters from packages auto artemis_pkg = pmb->packages.Get("artemis"); @@ -78,6 +89,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const bool do_dust = artemis_pkg->Param("do_dust"); const bool do_imc = artemis_pkg->Param("do_imc"); const bool do_moment = artemis_pkg->Param("do_moment"); + const bool do_mhd = artemis_pkg->Param("do_mhd"); PARTHENON_REQUIRE(do_gas, "The shock problem requires gas hydrodynamics!"); PARTHENON_REQUIRE(!(do_dust), "The shock problem does not permit dust hydrodynamics!"); auto eos_d = pmb->packages.Get("gas")->Param("eos_d"); @@ -94,7 +106,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } static auto desc = MakePackDescriptor( + rad::prim::energy, rad::prim::flux, field::face::B>( (pmb->resolved_packages).get()); auto v = desc.GetPack(md.get()); static auto desc_g = MakePackDescriptor((pmb->resolved_packages).get()); @@ -102,6 +114,15 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + IndexRange ib1 = pmb->cellbounds.GetBoundsI(IndexDomain::entire, TE::F1); + IndexRange jb1 = pmb->cellbounds.GetBoundsJ(IndexDomain::entire, TE::F1); + IndexRange kb1 = pmb->cellbounds.GetBoundsK(IndexDomain::entire, TE::F1); + IndexRange ib2 = pmb->cellbounds.GetBoundsI(IndexDomain::entire, TE::F2); + IndexRange jb2 = pmb->cellbounds.GetBoundsJ(IndexDomain::entire, TE::F2); + IndexRange kb2 = pmb->cellbounds.GetBoundsK(IndexDomain::entire, TE::F2); + IndexRange ib3 = pmb->cellbounds.GetBoundsI(IndexDomain::entire, TE::F3); + IndexRange jb3 = pmb->cellbounds.GetBoundsJ(IndexDomain::entire, TE::F3); + IndexRange kb3 = pmb->cellbounds.GetBoundsK(IndexDomain::entire, TE::F3); auto &pco = pmb->coords; // Shock parameters @@ -117,13 +138,17 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const bool upwind = (xi <= shkp.xdisc); const Real rho = upwind ? shkp.rhol : shkp.rhor; const Real vx = upwind ? shkp.vxl : shkp.vxr; - const Real T = upwind ? shkp.tl : shkp.tr; + const Real pres = upwind ? shkp.pl : shkp.pr; + const Real sie = (pres != Null()) + ? ArtemisUtils::EofPR(eos_d, pres, rho) + : eos_d.InternalEnergyFromDensityTemperature( + rho, upwind ? shkp.tl : shkp.tr); + const Real T = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); v(0, gas::prim::density(0), k, j, i) = rho; v(0, gas::prim::velocity(0), k, j, i) = vx; v(0, gas::prim::velocity(1), k, j, i) = 0.0; v(0, gas::prim::velocity(2), k, j, i) = 0.0; - v(0, gas::prim::sie(0), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(rho, T); + v(0, gas::prim::sie(0), k, j, i) = sie; if (do_moment) { v(0, rad::prim::energy(0), k, j, i) = ar * SQR(SQR(T)); v(0, rad::prim::flux(0), k, j, i) = 0.0; @@ -132,6 +157,30 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } }); + if (do_mhd) { + pmb->par_for( + "shock_b1", kb1.s, kb1.e, jb1.s, jb1.e, ib1.s, ib1.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + v(0, TE::F1, field::face::B(), k, j, i) = shkp.bx; + }); + pmb->par_for( + "shock_b2", kb2.s, kb2.e, jb2.s, jb2.e, ib2.s, ib2.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + geometry::Coords coords(cpars, pco, k, j, i); + const auto xi = coords.x1v(); + const bool upwind = (xi <= shkp.xdisc); + v(0, TE::F2, field::face::B(), k, j, i) = upwind ? shkp.byl : shkp.byr; + }); + pmb->par_for( + "shock_b3", kb3.s, kb3.e, jb3.s, jb3.e, ib3.s, ib3.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + geometry::Coords coords(cpars, pco, k, j, i); + const auto xi = coords.x1v(); + const bool upwind = (xi <= shkp.xdisc); + v(0, TE::F3, field::face::B(), k, j, i) = upwind ? shkp.bzl : shkp.bzr; + }); + } + if (do_imc) jaybenne::InitializeRadiation(md.get(), true); } @@ -146,6 +195,7 @@ inline void ShockInnerX1(std::shared_ptr> &mbd, bool coarse) auto artemis_pkg = pmb->packages.Get("artemis"); const bool do_moment = artemis_pkg->Param("do_moment"); + const bool do_mhd = artemis_pkg->Param("do_mhd"); auto shkp = artemis_pkg->Param("shock_params"); auto eos_d = pmb->packages.Get("gas")->Param("eos_d"); Real ar = Null(); @@ -157,25 +207,37 @@ inline void ShockInnerX1(std::shared_ptr> &mbd, bool coarse) static auto descriptors = ArtemisUtils::GetBoundaryPackDescriptorMap(mbd); + rad::prim::flux, field::cell::B>(mbd); auto v = descriptors[coarse].GetPack(mbd.get()); if (v.GetMaxNumberOfVars() == 0) return; + static auto descriptors_b = + ArtemisUtils::GetBoundaryPackDescriptorMap(mbd); + auto vb = descriptors_b[coarse].GetPack(mbd.get()); + const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds; pmb->par_for_bndry( - "ShockInnerX1", nb, IndexDomain::inner_x1, parthenon::TopologicalElement::CC, - coarse, false, + "ShockInnerX1", nb, IndexDomain::inner_x1, TE::CC, coarse, false, KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + const Real sie = + (shkp.pl != Null()) + ? ArtemisUtils::EofPR(eos_d, shkp.pl, shkp.rhol) + : eos_d.InternalEnergyFromDensityTemperature(shkp.rhol, shkp.tl); + const Real T = eos_d.TemperatureFromDensityInternalEnergy(shkp.rhol, sie); for (int n = 0; n < v.GetSize(0, gas::prim::density()); ++n) { v(0, gas::prim::density(n), k, j, i) = shkp.rhol; v(0, gas::prim::velocity(VI(n, 0)), k, j, i) = shkp.vxl; v(0, gas::prim::velocity(VI(n, 1)), k, j, i) = 0.0; v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = 0.0; - v(0, gas::prim::sie(n), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(shkp.rhol, shkp.tl); + v(0, gas::prim::sie(n), k, j, i) = sie; + } + if (do_mhd) { + v(0, TE::CC, field::cell::B(0), k, j, i) = shkp.bx; + v(0, TE::CC, field::cell::B(1), k, j, i) = shkp.byl; + v(0, TE::CC, field::cell::B(2), k, j, i) = shkp.bzl; } if (do_moment) { for (int n = 0; n < v.GetSize(0, rad::prim::energy()); ++n) { - v(0, rad::prim::energy(n), k, j, i) = ar * SQR(SQR(shkp.tl)); + v(0, rad::prim::energy(n), k, j, i) = ar * SQR(SQR(T)); v(0, rad::prim::flux(VI(n, 0)), k, j, i) = 0.0; v(0, rad::prim::flux(VI(n, 1)), k, j, i) = 0.0; v(0, rad::prim::flux(VI(n, 2)), k, j, i) = 0.0; @@ -183,6 +245,24 @@ inline void ShockInnerX1(std::shared_ptr> &mbd, bool coarse) } }); + if (do_mhd && vb.GetMaxNumberOfVars() > 0) { + pmb->par_for_bndry( + "ShockInnerX1::B1", nb, IndexDomain::inner_x1, TE::F1, coarse, false, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + vb(0, TE::F1, field::face::B(), k, j, i) = shkp.bx; + }); + pmb->par_for_bndry( + "ShockInnerX1::B2", nb, IndexDomain::inner_x1, TE::F2, coarse, false, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + vb(0, TE::F2, field::face::B(), k, j, i) = shkp.byl; + }); + pmb->par_for_bndry( + "ShockInnerX1::B3", nb, IndexDomain::inner_x1, TE::F3, coarse, false, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + vb(0, TE::F3, field::face::B(), k, j, i) = shkp.bzl; + }); + } + return; } @@ -197,6 +277,7 @@ inline void ShockOuterX1(std::shared_ptr> &mbd, bool coarse) auto artemis_pkg = pmb->packages.Get("artemis"); const bool do_moment = artemis_pkg->Param("do_moment"); + const bool do_mhd = artemis_pkg->Param("do_mhd"); auto shkp = artemis_pkg->Param("shock_params"); auto eos_d = pmb->packages.Get("gas")->Param("eos_d"); Real ar = Null(); @@ -208,25 +289,36 @@ inline void ShockOuterX1(std::shared_ptr> &mbd, bool coarse) static auto descriptors = ArtemisUtils::GetBoundaryPackDescriptorMap(mbd); + rad::prim::flux, field::cell::B>(mbd); auto v = descriptors[coarse].GetPack(mbd.get()); if (v.GetMaxNumberOfVars() == 0) return; + static auto descriptors_b = + ArtemisUtils::GetBoundaryPackDescriptorMap(mbd); + auto vb = descriptors_b[coarse].GetPack(mbd.get()); pmb->par_for_bndry( - "ShockOuterX1", nb, IndexDomain::outer_x1, parthenon::TopologicalElement::CC, - coarse, false, + "ShockOuterX1", nb, IndexDomain::outer_x1, TE::CC, coarse, false, KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + const Real sie = + (shkp.pr != Null()) + ? ArtemisUtils::EofPR(eos_d, shkp.pr, shkp.rhor) + : eos_d.InternalEnergyFromDensityTemperature(shkp.rhor, shkp.tr); + const Real T = eos_d.TemperatureFromDensityInternalEnergy(shkp.rhor, sie); for (int n = 0; n < v.GetSize(0, gas::prim::density()); ++n) { v(0, gas::prim::density(n), k, j, i) = shkp.rhor; v(0, gas::prim::velocity(VI(n, 0)), k, j, i) = shkp.vxr; v(0, gas::prim::velocity(VI(n, 1)), k, j, i) = 0.0; v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = 0.0; - v(0, gas::prim::sie(n), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(shkp.rhor, shkp.tr); + v(0, gas::prim::sie(n), k, j, i) = sie; + } + if (do_mhd) { + v(0, TE::CC, field::cell::B(0), k, j, i) = shkp.bx; + v(0, TE::CC, field::cell::B(1), k, j, i) = shkp.byr; + v(0, TE::CC, field::cell::B(2), k, j, i) = shkp.bzr; } if (do_moment) { for (int n = 0; n < v.GetSize(0, rad::prim::energy()); ++n) { - v(0, rad::prim::energy(n), k, j, i) = ar * SQR(SQR(shkp.tr)); + v(0, rad::prim::energy(n), k, j, i) = ar * SQR(SQR(T)); v(0, rad::prim::flux(VI(n, 0)), k, j, i) = 0.0; v(0, rad::prim::flux(VI(n, 1)), k, j, i) = 0.0; v(0, rad::prim::flux(VI(n, 2)), k, j, i) = 0.0; @@ -234,6 +326,24 @@ inline void ShockOuterX1(std::shared_ptr> &mbd, bool coarse) } }); + if (do_mhd && vb.GetMaxNumberOfVars() > 0) { + pmb->par_for_bndry( + "ShockOuterX1::B1", nb, IndexDomain::outer_x1, TE::F1, coarse, false, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + vb(0, TE::F1, field::face::B(), k, j, i) = shkp.bx; + }); + pmb->par_for_bndry( + "ShockOuterX1::B2", nb, IndexDomain::outer_x1, TE::F2, coarse, false, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + vb(0, TE::F2, field::face::B(), k, j, i) = shkp.byr; + }); + pmb->par_for_bndry( + "ShockOuterX1::B3", nb, IndexDomain::outer_x1, TE::F3, coarse, false, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + vb(0, TE::F3, field::face::B(), k, j, i) = shkp.bzr; + }); + } + return; } diff --git a/src/utils/artemis_utils.cpp b/src/utils/artemis_utils.cpp index 5b014c38..5ca52544 100644 --- a/src/utils/artemis_utils.cpp +++ b/src/utils/artemis_utils.cpp @@ -13,11 +13,146 @@ // C++ headers #include "artemis_utils.hpp" +#include "geometry/geometry.hpp" #include "nbody/nbody_utils.hpp" #include "units.hpp" namespace ArtemisUtils { +namespace { + +template +Real ComputeLocalMaxAbsFaceDivBImpl(const BlockList_t &blocks, + const geometry::CoordParams &cpars, const int ndim) { + using TE = TopologicalElement; + constexpr int f1 = static_cast(TE::F1) % 3; + constexpr int f2 = static_cast(TE::F2) % 3; + constexpr int f3 = static_cast(TE::F3) % 3; + constexpr Real tiny = 1.0e-300; + + Real max_divb = 0.0; + for (const auto &pmb : blocks) { + auto &md = pmb->meshblock_data.Get(); + if (!md->HasVariable("field.face.B")) continue; + + auto var = md->GetVarPtr("field.face.B"); + if (!var->IsAllocated()) continue; + + auto b = var->data; + auto pco = pmb->coords; + const auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + const auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + const auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + + Real block_max_divb = 0.0; + if (ndim == 3) { + parthenon::par_reduce( + parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), kb.s, + kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i, Real &lmax_divb) { + geometry::Coords coords(cpars, pco, k, j, i); + const Real vol = coords.Volume(); + const auto ax1 = coords.GetFaceAreaX1(); + const auto ax2 = coords.GetFaceAreaX2(); + const auto ax3 = coords.GetFaceAreaX3(); + const Real divb = ((ax1[1] * b(f1, 0, 0, 0, k, j, i + 1) - + ax1[0] * b(f1, 0, 0, 0, k, j, i)) + + (ax2[1] * b(f2, 0, 0, 0, k, j + 1, i) - + ax2[0] * b(f2, 0, 0, 0, k, j, i)) + + (ax3[1] * b(f3, 0, 0, 0, k + 1, j, i) - + ax3[0] * b(f3, 0, 0, 0, k, j, i))) / + (vol + tiny); + const Real abs_divb = fabs(divb); + if (abs_divb > lmax_divb) lmax_divb = abs_divb; + }, + Kokkos::Max(block_max_divb)); + } else if (ndim == 2) { + const int k = kb.s; + parthenon::par_reduce( + parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), jb.s, + jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int j, const int i, Real &lmax_divb) { + geometry::Coords coords(cpars, pco, k, j, i); + const Real vol = coords.Volume(); + const auto ax1 = coords.GetFaceAreaX1(); + const auto ax2 = coords.GetFaceAreaX2(); + const Real divb = ((ax1[1] * b(f1, 0, 0, 0, k, j, i + 1) - + ax1[0] * b(f1, 0, 0, 0, k, j, i)) + + (ax2[1] * b(f2, 0, 0, 0, k, j + 1, i) - + ax2[0] * b(f2, 0, 0, 0, k, j, i))) / + (vol + tiny); + const Real abs_divb = fabs(divb); + if (abs_divb > lmax_divb) lmax_divb = abs_divb; + }, + Kokkos::Max(block_max_divb)); + } else { + const int k = kb.s; + const int j = jb.s; + parthenon::par_reduce( + parthenon::loop_pattern_flatrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), + ib.s, ib.e, + KOKKOS_LAMBDA(const int i, Real &lmax_divb) { + geometry::Coords coords(cpars, pco, k, j, i); + const Real vol = coords.Volume(); + const auto ax1 = coords.GetFaceAreaX1(); + const Real divb = (ax1[1] * b(f1, 0, 0, 0, k, j, i + 1) - + ax1[0] * b(f1, 0, 0, 0, k, j, i)) / + (vol + tiny); + const Real abs_divb = fabs(divb); + if (abs_divb > lmax_divb) lmax_divb = abs_divb; + }, + Kokkos::Max(block_max_divb)); + } + + max_divb = std::max(max_divb, block_max_divb); + } + + return max_divb; +} + +void PrintMeshMaxAbsDivB(Mesh *pm, const char *label) { + auto artemis = pm->packages.Get("artemis"); + if (!artemis->Param("do_mhd")) return; + + const auto coords = artemis->Param("coords"); + const auto &cpars = artemis->Param("coord_params"); + const int ndim = artemis->Param("ndim"); + + Real max_divb = 0.0; + if (coords == Coordinates::cartesian) { + max_divb = ComputeLocalMaxAbsFaceDivBImpl(pm->block_list, + cpars, ndim); + } else if (coords == Coordinates::spherical1D) { + max_divb = ComputeLocalMaxAbsFaceDivBImpl(pm->block_list, + cpars, ndim); + } else if (coords == Coordinates::spherical2D) { + max_divb = ComputeLocalMaxAbsFaceDivBImpl(pm->block_list, + cpars, ndim); + } else if (coords == Coordinates::spherical3D) { + max_divb = ComputeLocalMaxAbsFaceDivBImpl(pm->block_list, + cpars, ndim); + } else if (coords == Coordinates::cylindrical) { + max_divb = ComputeLocalMaxAbsFaceDivBImpl(pm->block_list, + cpars, ndim); + } else if (coords == Coordinates::axisymmetric) { + max_divb = ComputeLocalMaxAbsFaceDivBImpl(pm->block_list, + cpars, ndim); + } else { + PARTHENON_FAIL("Invalid Artemis coordinate system for divB diagnostic."); + } + +#ifdef MPI_PARALLEL + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &max_divb, 1, MPI_PARTHENON_REAL, + MPI_MAX, MPI_COMM_WORLD)); +#endif + if (parthenon::Globals::my_rank == 0 && max_divb > 0.0) { + std::cout << "AMR remesh: max |divB| [" << label << "] = " << std::scientific + << max_divb << std::defaultfloat << std::endl; + } +} + +} // namespace + //---------------------------------------------------------------------------------------- //! \fn void ArtemisUtils::PrintArtemisConfiguration //! \brief @@ -92,6 +227,17 @@ void PrintArtemisConfiguration(Packages_t &packages) { } } +void PreStepDiagnosticsRemeshDivB(SimTime const &simtime, MeshData *rc) { + auto pm = rc->GetMeshPointer(); + if ((simtime.ncycle > 0) && pm->modified) { + PrintMeshMaxAbsDivB(pm, "post-remesh"); + } +} + +void PostStepDiagnosticsRemeshDivB(SimTime const &, MeshData *rc) { + PrintMeshMaxAbsDivB(rc->GetMeshPointer(), "pre-remesh"); +} + //---------------------------------------------------------------------------------------- //! \fn void ArtemisUtils::EnrollArtemisRefinementOps //! \brief Registers custom prolongation and restriction operators on provided Metadata @@ -222,6 +368,158 @@ void EnrollArtemisRefinementOps(parthenon::Metadata &m, Coordinates coords, } } +//---------------------------------------------------------------------------------------- +//! \fn void ArtemisUtils::EnrollArtemisFaceRefinementOps +//! \brief Registers custom face-centered prolongation and restriction operators. +void EnrollArtemisFaceRefinementOps(parthenon::Metadata &m, Coordinates coords, + const bool log, const bool use_minmod_slope) { + typedef Coordinates G; + + if (coords == G::cartesian) { + if (use_minmod_slope) { + m.RegisterRefinementOps, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } else if (coords == G::spherical1D) { + if (log) { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } else { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } + } else if (coords == G::spherical2D) { + if (log) { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } else { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } + } else if (coords == G::spherical3D) { + if (log) { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } else { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } + } else if (coords == G::cylindrical) { + if (log) { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } else { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } + } else if (coords == G::axisymmetric) { + if (log) { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } else { + if (use_minmod_slope) { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } else { + m.RegisterRefinementOps< + ArtemisUtils::ProlongateShared, + ArtemisUtils::RestrictAverage, + ArtemisUtils::ProlongateTothAndRoe>(); + } + } + } else { + PARTHENON_FAIL("Invalid artemis/coordinate system!"); + } +} + //---------------------------------------------------------------------------------------- //! \fn std::vector> NBody::loadtxt //! \brief diff --git a/src/utils/artemis_utils.hpp b/src/utils/artemis_utils.hpp index c81a797b..5885fc9a 100644 --- a/src/utils/artemis_utils.hpp +++ b/src/utils/artemis_utils.hpp @@ -46,14 +46,17 @@ KOKKOS_FORCEINLINE_FUNCTION Real VDot(const V1 &a, const V2 &b) { //---------------------------------------------------------------------------------------- //! \fn Real ArtemisUtils::DualEnergySIE(vmesh, const int b, const int n, const int k, //! const int j, const int i, const Real de_switch, -//! const Real hx[3]) +//! const Real hx[3], const Real emag) //! \brief Returns appropriate specific internal energy variable based on de_switch +//! If `de_switch <= 0`, the dual-energy switch is disabled and the auxiliary internal +//! energy is used. //! NOTE(@pdmullen): Floors should be handled outside this function call template KOKKOS_FORCEINLINE_FUNCTION Real DualEnergySIE(T &vmesh, const int b, const int n, const int k, const int j, const int i, const Real de_switch, - const std::array &hx) { + const std::array &hx, + const Real emag = 0.0) { // Extract state vector const Real invd = 1.0 / vmesh(b, gas::cons::density(n), k, j, i); const Real &rv1 = vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) / hx[0]; @@ -64,9 +67,10 @@ KOKKOS_FORCEINLINE_FUNCTION Real DualEnergySIE(T &vmesh, const int b, const int const Real ke = 0.5 * invd * (SQR(rv1) + SQR(rv2) + SQR(rv3)); // Calculate conserved representation of internal energy - const Real ut_sie = invd * (u_e - ke); - const bool dual_switch = (ut_sie > invd * de_switch * u_e); - return (dual_switch)*ut_sie + (!dual_switch) * invd * u_u; + const Real ut_sie = invd * (u_e - ke - emag); + const Real u_hyd = u_e - emag; + const bool use_total = (de_switch > 0.0) && (ut_sie > invd * de_switch * u_hyd); + return (use_total)*ut_sie + (!use_total) * invd * u_u; } //---------------------------------------------------------------------------------------- @@ -173,8 +177,12 @@ struct SumMyArray { //! Defined in artemis_utils.cpp //! NOTE(@pdmullen): We should likely move everything above to implementation file too... void PrintArtemisConfiguration(Packages_t &packages); +void PreStepDiagnosticsRemeshDivB(SimTime const &simtime, MeshData *rc); +void PostStepDiagnosticsRemeshDivB(SimTime const &simtime, MeshData *rc); void EnrollArtemisRefinementOps(parthenon::Metadata &m, Coordinates coords, const bool log, const bool use_minmod_slope = true); +void EnrollArtemisFaceRefinementOps(parthenon::Metadata &m, Coordinates coords, + const bool log, const bool use_minmod_slope = true); std::vector> loadtxt(std::string fname); // 4D outer parallel loop using Kokkos Teams diff --git a/src/utils/fluxes/fluid_fluxes.hpp b/src/utils/fluxes/fluid_fluxes.hpp index 8d35c205..cc79ab48 100644 --- a/src/utils/fluxes/fluid_fluxes.hpp +++ b/src/utils/fluxes/fluid_fluxes.hpp @@ -16,6 +16,7 @@ // Artemis includes #include "artemis.hpp" #include "geometry/geometry.hpp" +#include "mhd/mhd.hpp" #include "rotating_frame/rotating_frame.hpp" #include "utils/artemis_utils.hpp" #include "utils/fluxes/reconstruction/reconstruction.hpp" @@ -64,6 +65,7 @@ ScaleMomentumFlux(parthenon::team_mbr_t const &member, const geometry::CoordPara q.flux(b, DIR, IVZ, k, j, i) *= hx[2]; }); } + return; } @@ -92,8 +94,11 @@ TaskStatus CalculateFluxesImpl(MeshData *md, PKG &pkg, PRIM vp, FLUX vflx, eos = pkg->template Param("eos_d"); } - const auto &cpars = - pm->packages.Get("artemis")->template Param("coord_params"); + const auto &artemis_pkg = pm->packages.Get("artemis"); + const auto &cpars = artemis_pkg->template Param("coord_params"); + const auto do_mhd = artemis_pkg->template Param("do_mhd"); + const Real mu0_code = + do_mhd ? pm->packages.Get("mhd")->template Param("mu0_code") : 1.0; // Speed of light (and reduced), if used Real chat = Null(); @@ -124,6 +129,14 @@ TaskStatus CalculateFluxesImpl(MeshData *md, PKG &pkg, PRIM vp, FLUX vflx, // X1-Flux int il = ib.s, iu = ib.e + 1; int jl = jb.s, ju = jb.e, kl = kb.s, ku = kb.e; + if constexpr (F == Fluid::gas) { + if (do_mhd) { + jl -= multi_d; + ju += multi_d; + kl -= three_d; + ku += three_d; + } + } parthenon::par_for_outer( DEFAULT_OUTER_LOOP_PATTERN, "CalculateFluxes::X1-Flux", DevExecSpace(), scr_size, scr_level, 0, md->NumBlocks() - 1, kl, ku, jl, ju, @@ -136,22 +149,35 @@ TaskStatus CalculateFluxesImpl(MeshData *md, PKG &pkg, PRIM vp, FLUX vflx, recon(mbr, cpars, b, k, j, il - 1, iu, vp, vg, wl, wr); mbr.team_barrier(); - post_recon(eos, dfloor, siefloor, mbr, X1DIR, b, k, j, il - 1, iu, vp, wl, wr); + post_recon(eos, dfloor, siefloor, do_mhd, mbr, X1DIR, b, k, j, il - 1, iu, vp, + vflx, wl, wr); mbr.team_barrier(); // Compute fluxes over[is, ie + 1] RiemannSolver riemann; - riemann(eos, c, chat, mbr, b, k, j, il, iu, X1DIR, wl, wr, vp, vflx, vface); + riemann(eos, c, chat, mu0_code, do_mhd, mbr, b, k, j, il, iu, X1DIR, wl, wr, vp, + vflx, vface); mbr.team_barrier(); // Scale X1-momentum flux by appropriate scale factor for coord system ScaleMomentumFlux(mbr, cpars, b, k, j, il, iu, vg, vflx); + if constexpr (F == Fluid::gas) { + if (do_mhd) MHD::ScaleMHDFlux(mbr, cpars, b, k, j, il, iu, vg, vp); + } }); // X2-Flux if (multi_d) { jl = jb.s - 1, ju = jb.e + 1; il = ib.s, iu = ib.e, kl = kb.s, ku = kb.e; + if constexpr (F == Fluid::gas) { + if (do_mhd) { + il -= 1; + iu += 1; + kl -= three_d; + ku += three_d; + } + } scr_size = ScratchPad2D::shmem_size(nvars, ncells1) * 3; parthenon::par_for_outer( DEFAULT_OUTER_LOOP_PATTERN, "CalculateFluxes::X2-Flux", DevExecSpace(), scr_size, @@ -176,18 +202,23 @@ TaskStatus CalculateFluxesImpl(MeshData *md, PKG &pkg, PRIM vp, FLUX vflx, recon(mbr, cpars, b, k, j, il, iu, vp, vg, wl_jp1, wr); mbr.team_barrier(); - post_recon(eos, dfloor, siefloor, mbr, X2DIR, b, k, j, il, iu, vp, wl_jp1, - wr); + post_recon(eos, dfloor, siefloor, do_mhd, mbr, X2DIR, b, k, j, il, iu, vp, + vflx, wl_jp1, wr); mbr.team_barrier(); if (j > jl) { // compute fluxes over [js,je+1] RiemannSolver riemann; - riemann(eos, c, chat, mbr, b, k, j, il, iu, X2DIR, wl, wr, vp, vflx, vface); + riemann(eos, c, chat, mu0_code, do_mhd, mbr, b, k, j, il, iu, X2DIR, wl, wr, + vp, vflx, vface); mbr.team_barrier(); // Scale X2-momentum flux by appropriate scale factor for coord system ScaleMomentumFlux(mbr, cpars, b, k, j, il, iu, vg, vflx); + if constexpr (F == Fluid::gas) { + if (do_mhd) + MHD::ScaleMHDFlux(mbr, cpars, b, k, j, il, iu, vg, vp); + } } } }); @@ -197,6 +228,14 @@ TaskStatus CalculateFluxesImpl(MeshData *md, PKG &pkg, PRIM vp, FLUX vflx, if (three_d) { kl = kb.s - 1, ku = kb.e + 1; il = ib.s, iu = ib.e, jl = jb.s, ju = jb.e; + if constexpr (F == Fluid::gas) { + if (do_mhd) { + il -= 1; + iu += 1; + jl -= 1; + ju += 1; + } + } scr_size = ScratchPad2D::shmem_size(nvars, ncells1) * 3; parthenon::par_for_outer( DEFAULT_OUTER_LOOP_PATTERN, "Hydro::X3-Flux", DevExecSpace(), scr_size, scr_level, @@ -221,18 +260,23 @@ TaskStatus CalculateFluxesImpl(MeshData *md, PKG &pkg, PRIM vp, FLUX vflx, recon(mbr, cpars, b, k, j, il, iu, vp, vg, wl_kp1, wr); mbr.team_barrier(); - post_recon(eos, dfloor, siefloor, mbr, X3DIR, b, k, j, il, iu, vp, wl_kp1, - wr); + post_recon(eos, dfloor, siefloor, do_mhd, mbr, X3DIR, b, k, j, il, iu, vp, + vflx, wl_kp1, wr); mbr.team_barrier(); // compute fluxes over [ks,ke+1] if (k > kl) { RiemannSolver riemann; - riemann(eos, c, chat, mbr, b, k, j, il, iu, X3DIR, wl, wr, vp, vflx, vface); + riemann(eos, c, chat, mu0_code, do_mhd, mbr, b, k, j, il, iu, X3DIR, wl, wr, + vp, vflx, vface); mbr.team_barrier(); // Scale X3-momentum flux by appropriate scale factor for coord system ScaleMomentumFlux(mbr, cpars, b, k, j, il, iu, vg, vflx); + if constexpr (F == Fluid::gas) { + if (do_mhd) + MHD::ScaleMHDFlux(mbr, cpars, b, k, j, il, iu, vg, vp); + } } } }); @@ -270,10 +314,14 @@ TaskStatus FluxSourceImpl(MeshData *md, PKG &pkg, PRIM vp, CONS vcons, FAC if constexpr (F == Fluid::radiation) { hcchat = 0.5 * pkg->template Param("c") * pkg->template Param("chat"); } - const auto &cpars = md->GetParentPointer() - ->packages.Get("artemis") - ->template Param("coord_params"); + const auto &artemis_pkg = md->GetParentPointer()->packages.Get("artemis"); + const auto &cpars = artemis_pkg->template Param("coord_params"); + const auto do_mhd = artemis_pkg->template Param("do_mhd"); + const Real mu0_ = + do_mhd + ? md->GetParentPointer()->packages.Get("mhd")->template Param("mu0_code") + : 1.0; // Apply flux sources parthenon::par_for( DEFAULT_LOOP_PATTERN, "FluxSourceTerms", parthenon::DevExecSpace(), 0, @@ -342,6 +390,7 @@ TaskStatus FluxSourceImpl(MeshData *md, PKG &pkg, PRIM vp, CONS vcons, FAC // Add the "flux source terms" for (int n = 0; n < nspecies; ++n) { + const bool mhd = (F == Fluid::gas) && do_mhd && (n == 0); const int IMX = VI(n, 0); const int IMY = VI(n, 1); const int IMZ = VI(n, 2); @@ -352,7 +401,7 @@ TaskStatus FluxSourceImpl(MeshData *md, PKG &pkg, PRIM vp, CONS vcons, FAC const int IEG = nspecies * 3 + n; // may not be used // Pressure gradient force for gas and radiation - if constexpr (F == Fluid::gas || F == Fluid::radiation) { + if constexpr (F == Fluid::radiation) { // Pressure gradient force vc_(b, IMX, k, j, i) += dtdx[0] * (vp_.flux(b, d1, IPR, k, j, i) - vp_.flux(b, d1, IPR, k, j, i + 1)); @@ -360,10 +409,26 @@ TaskStatus FluxSourceImpl(MeshData *md, PKG &pkg, PRIM vp, CONS vcons, FAC vp_.flux(b, d2, IPR, k, j + multi_d, i)); vc_(b, IMZ, k, j, i) += dtdx[2] * (vp_.flux(b, d3, IPR, k, j, i) - vp_.flux(b, d3, IPR, k + three_d, j, i)); - } + } else if constexpr (F == Fluid::gas) { + // Pressure gradient force + Real Pl = vp_.flux(b, d1, IPR, k, j, i) + + (mhd ? vp_.flux(b, d1, field::cell::energy(), k, j, i) : 0.0); + Real Pr = vp_.flux(b, d1, IPR, k, j, i + 1) + + (mhd ? vp_.flux(b, d1, field::cell::energy(), k, j, i + 1) : 0.0); + vc_(b, IMX, k, j, i) += dtdx[0] * (Pl - Pr); + + Pl = vp_.flux(b, d2, IPR, k, j, i) + + (mhd ? vp_.flux(b, d2, field::cell::energy(), k, j, i) : 0.0); + Pr = vp_.flux(b, d2, IPR, k, j + multi_d, i) + + (mhd ? vp_.flux(b, d2, field::cell::energy(), k, j + multi_d, i) : 0.0); + vc_(b, IMY, k, j, i) += dtdx[1] * (Pl - Pr); + + Pl = vp_.flux(b, d3, IPR, k, j, i) + + (mhd ? vp_.flux(b, d3, field::cell::energy(), k, j, i) : 0.0); + Pr = vp_.flux(b, d3, IPR, k + three_d, j, i) + + (mhd ? vp_.flux(b, d3, field::cell::energy(), k + three_d, j, i) : 0.0); + vc_(b, IMZ, k, j, i) += dtdx[2] * (Pl - Pr); - // pdV source for gas internal energy equation - if constexpr (F == Fluid::gas) { // pdV source term const auto &ax1 = coords.GetFaceAreaX1(vg, b, k, j, i); const auto &ax2 = coords.GetFaceAreaX2(vg, b, k, j, i); @@ -391,6 +456,7 @@ TaskStatus FluxSourceImpl(MeshData *md, PKG &pkg, PRIM vp, CONS vcons, FAC [[maybe_unused]] const auto x1dep_ = x1dep; [[maybe_unused]] const auto x2dep_ = x2dep; [[maybe_unused]] const auto x3dep_ = x3dep; + [[maybe_unused]] const auto mu0 = mu0_; if constexpr (G != Coordinates::cartesian) { // Extract primitive weighted timestep Real wdt = vp_(b, n, k, j, i) * dt; @@ -403,17 +469,32 @@ TaskStatus FluxSourceImpl(MeshData *md, PKG &pkg, PRIM vp, CONS vcons, FAC const Real ff = std::sqrt(SQR(fx) + SQR(fy) + SQR(fz)); const Real chi = Moments::ThriceEddingtonFactor(ff); wdt *= ((chi - 1.) / (ff + Fuzz())) * hcchat_; + } else if constexpr (F == Fluid::gas) { + // Update momenta with mhd + const Real t1 = SQR(vp_(b, IVX, k, j, i) + rfv[0]) - + (mhd ? SQR(vp_(b, field::cell::B(0), k, j, i)) / mu0 : 0.0); + const Real t2 = SQR(vp_(b, IVY, k, j, i) + rfv[1]) - + (mhd ? SQR(vp_(b, field::cell::B(1), k, j, i)) / mu0 : 0.0); + const Real t3 = SQR(vp_(b, IVZ, k, j, i) + rfv[2]) - + (mhd ? SQR(vp_(b, field::cell::B(2), k, j, i)) / mu0 : 0.0); + vc_(b, IMX, k, j, i) += + x1dep_ * wdt * (dh1[0] * t1 + dh1[1] * t2 + dh1[2] * t3); + vc_(b, IMY, k, j, i) += + x2dep_ * wdt * (dh2[0] * t1 + dh2[1] * t2 + dh2[2] * t3); + vc_(b, IMZ, k, j, i) += + x3dep_ * wdt * (dh3[0] * t1 + dh3[1] * t2 + dh3[2] * t3); + } else { + // Update momenta + const Real t1 = SQR(vp_(b, IVX, k, j, i) + rfv[0]); + const Real t2 = SQR(vp_(b, IVY, k, j, i) + rfv[1]); + const Real t3 = SQR(vp_(b, IVZ, k, j, i) + rfv[2]); + vc_(b, IMX, k, j, i) += + x1dep_ * wdt * (dh1[0] * t1 + dh1[1] * t2 + dh1[2] * t3); + vc_(b, IMY, k, j, i) += + x2dep_ * wdt * (dh2[0] * t1 + dh2[1] * t2 + dh2[2] * t3); + vc_(b, IMZ, k, j, i) += + x3dep_ * wdt * (dh3[0] * t1 + dh3[1] * t2 + dh3[2] * t3); } - - // Update momenta - // clang-format off - const Real t1 = SQR(vp_(b, IVX, k, j, i) + rfv[0]); - const Real t2 = SQR(vp_(b, IVY, k, j, i) + rfv[1]); - const Real t3 = SQR(vp_(b, IVZ, k, j, i) + rfv[2]); - vc_(b, IMX, k, j, i) += x1dep_ * wdt * (dh1[0]*t1 + dh1[1]*t2 + dh1[2]*t3); - vc_(b, IMY, k, j, i) += x2dep_ * wdt * (dh2[0]*t1 + dh2[1]*t2 + dh2[2]*t3); - vc_(b, IMZ, k, j, i) += x3dep_ * wdt * (dh3[0]*t1 + dh3[1]*t2 + dh3[2]*t3); - // clang-format on } } }); @@ -477,6 +558,13 @@ TaskStatus CalculateFluxesRiemannSelect(MeshData *md, PKG &pkg, PRIM vp, F pcm); } else if (riemann_method == R::llf) { return CalculateFluxesReconSelect(md, pkg, vp, vflx, vface, vg, pcm); + } else if (riemann_method == R::hlld) { + if constexpr (F == Fluid::gas) { + return CalculateFluxesReconSelect(md, pkg, vp, vflx, vface, vg, + pcm); + } else { + PARTHENON_FAIL("HLLD solver only supports ideal MHD"); + } } else { PARTHENON_FAIL("Riemann solver not recognized!"); } diff --git a/src/utils/fluxes/reconstruction/reconstruction.hpp b/src/utils/fluxes/reconstruction/reconstruction.hpp index c1a7011b..6fcc35af 100644 --- a/src/utils/fluxes/reconstruction/reconstruction.hpp +++ b/src/utils/fluxes/reconstruction/reconstruction.hpp @@ -47,12 +47,13 @@ struct ReconGradient { //! \class TaskStatus ArtemisUtils::post_recon //! \brief Utility to apply floors, make thermodynamically consistent, or zero radiation //! fluxes when near ~round-off -template +template KOKKOS_INLINE_FUNCTION void -post_recon(const EOS &eos, const Real dfloor, const Real siefloor, +post_recon(const EOS &eos, const Real dfloor, const Real siefloor, const bool do_mhd, parthenon::team_mbr_t const &member, const int dir, const int b, const int k, - const int j, const int il, const int iu, const V &q, + const int j, const int il, const int iu, const V &q, const VC &qc, parthenon::ScratchPad2D &ql, parthenon::ScratchPad2D &qr) { + using TE = parthenon::TopologicalElement; if constexpr (F == Fluid::radiation) { const int nspecies = q.GetSize(b, rad::prim::energy()); for (int n = 0; n < nspecies; ++n) { @@ -84,7 +85,7 @@ post_recon(const EOS &eos, const Real dfloor, const Real siefloor, Real &dL = ql(IDN, ipl); Real &pL = ql(IPR, ipl); Real &eL = ql(ISE, ipl); - Real &bL = qr(IBL, ipl); + Real &bL = ql(IBL, ipl); Real &dR = qr(IDN, i); Real &pR = qr(IPR, i); Real &eR = qr(ISE, i); @@ -107,6 +108,18 @@ post_recon(const EOS &eos, const Real dfloor, const Real siefloor, } }); } + // Replace the reconstructed B on the face in direction dir with the cons face value + if (do_mhd) { + const int n = nspecies * 7 + dir - 1; + // or maybe TE fd = TE::F1 + (dir-1)? + TE fd = (dir == 1) ? TE::F1 : ((dir == 2) ? TE::F2 : TE::F3); + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, member, il, iu, [&](const int i) { + const int ipl = i + (dir == 1); + ql(n, ipl) = qc(b, fd, field::face::B(), k + (dir == 3), j + (dir == 2), ipl); + qr(n, i) = qc(b, fd, field::face::B(), k, j, i); + }); + } } } diff --git a/src/utils/fluxes/riemann/hllc.hpp b/src/utils/fluxes/riemann/hllc.hpp index 9cd922f9..c421637c 100644 --- a/src/utils/fluxes/riemann/hllc.hpp +++ b/src/utils/fluxes/riemann/hllc.hpp @@ -47,13 +47,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const EOS &eos, const Real c, const Real chat, const Real mu0, + const bool do_mhd, parthenon::team_mbr_t const &member, const int b, + const int k, const int j, const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); @@ -193,13 +193,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const EOS &eos, const Real c, const Real chat, const Real mu0, + const bool do_mhd, parthenon::team_mbr_t const &member, const int b, + const int k, const int j, const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); diff --git a/src/utils/fluxes/riemann/hlld.hpp b/src/utils/fluxes/riemann/hlld.hpp new file mode 100644 index 00000000..305e7f74 --- /dev/null +++ b/src/utils/fluxes/riemann/hlld.hpp @@ -0,0 +1,550 @@ +//======================================================================================== +// (C) (or copyright) 2026. 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") +//======================================================================================== +//! \file hlld.hpp +//! \brief Contains HLLD Riemann solver for hydrodynamics +//! +//! Computes fluxes using the Harten-Lax-vanLeer-Discontinuities (HLLD) Riemann solver. +//! +#ifndef UTILS_FLUXES_RIEMANN_HLLD_HPP_ +#define UTILS_FLUXES_RIEMANN_HLLD_HPP_ + +// C++ headers +#include +#include + +// Artemis headers +#include "artemis.hpp" +#include "radiation/moments/moments.hpp" +#include "utils/eos/eos.hpp" + +// NOTE(AMD): The following is mostly taken directly from the open-source Athena++/AthenaK +// software + +namespace ArtemisUtils { +//---------------------------------------------------------------------------------------- +//! \class ArtemisUtils::RiemannSolver +//! \brief The HLLD Riemann solver for ideal gas +template +struct RiemannSolver> { + template + KOKKOS_INLINE_FUNCTION void + operator()(const EOS &eos, const Real c, const Real chat, const Real mu0, + const bool do_mhd, parthenon::team_mbr_t const &member, const int b, + const int k, const int j, const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { + using TE = parthenon::TopologicalElement; + // Check sensibility of flux direction + PARTHENON_REQUIRE(do_mhd, "HLLD solver only implemented for MHD!"); + PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); + auto fdir = (dir == 1) ? TE::F1 : ((dir == 2) ? TE::F2 : TE::F3); + + // Obtain number of species + int nspecies = q.GetSize(b, gas::cons::density()); + + const int IBX = nspecies * 7 + (dir - 1); + const int IBY = nspecies * 7 + ((dir - 1) + 1) % 3; + const int IBZ = nspecies * 7 + ((dir - 1) + 2) % 3; + const int IBM = nspecies * 7 + 3; + const int IBXG = dir - 1; + const int IBYG = dir % 3; + const int IBZG = (dir + 1) % 3; + + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, member, il, iu, [&](const int i) { + // Create local references for L/R states (helps compiler vectorize) + const int n = 0; + const int IDN = 0; + const int ivx = nspecies + ((dir - 1)); + const int ivy = nspecies + ((dir - 1) + 1) % 3; + const int ivz = nspecies + ((dir - 1) + 2) % 3; + // Unused indices for dust hydrodynamics + const int IPR = nspecies * 4; + const int ISE = nspecies * 5; + const int IBL = nspecies * 6; + const int IEN = IPR; + const int IEG = ISE; + Real &wl_idn = wl(IDN, i); + Real &wl_ivx = wl(ivx, i); + Real &wl_ivy = wl(ivy, i); + Real &wl_ivz = wl(ivz, i); + Real &wl_ipr = wl(IPR, i); + Real &wl_ise = wl(ISE, i); + Real &wl_ibl = wl(IBL, i); + + Real &wr_idn = wr(IDN, i); + Real &wr_ivx = wr(ivx, i); + Real &wr_ivy = wr(ivy, i); + Real &wr_ivz = wr(ivz, i); + Real &wr_ipr = wr(IPR, i); + Real &wr_ise = wr(ISE, i); + Real &wr_ibl = wr(IBL, i); + + Real &wl_ibx = wl(IBX, i); + Real &wl_iby = wl(IBY, i); + Real &wl_ibz = wl(IBZ, i); + Real &wr_ibx = wr(IBX, i); + Real &wr_iby = wr(IBY, i); + Real &wr_ibz = wr(IBZ, i); + + const Real small = 1.0e-20; + const Real eps = 1.0e-12; + + Real frho = Null(); + Real fmx = Null(); + Real fmy = Null(); + Real fmz = Null(); + Real fe = Null(); + Real pface = 0.0; + Real pmag_face = 0.0; + Real fby = 0.0; + Real fbz = 0.0; + + const Real bxi = 0.5 * (wl_ibx + wr_ibx); + const Real sqrt_mu0 = std::sqrt(mu0); + const Real inv_sqrt_mu0 = 1.0 / sqrt_mu0; + + const Real bxi_n = bxi * inv_sqrt_mu0; + const Real wl_iby_n = wl_iby * inv_sqrt_mu0; + const Real wl_ibz_n = wl_ibz * inv_sqrt_mu0; + const Real wr_iby_n = wr_iby * inv_sqrt_mu0; + const Real wr_ibz_n = wr_ibz * inv_sqrt_mu0; + const Real bxsq_n = SQR(bxi_n); + const Real pbl = 0.5 * (bxsq_n + SQR(wl_iby_n) + SQR(wl_ibz_n)); + const Real pbr = 0.5 * (bxsq_n + SQR(wr_iby_n) + SQR(wr_ibz_n)); + const Real ptl = wl_ipr + pbl; + const Real ptr = wr_ipr + pbr; + const Real vdotBl = wl_ivx * bxi_n + wl_ivy * wl_iby_n + wl_ivz * wl_ibz_n; + const Real vdotBr = wr_ivx * bxi_n + wr_ivy * wr_iby_n + wr_ivz * wr_ibz_n; + const Real el = wl_idn * wl_ise + + 0.5 * wl_idn * (SQR(wl_ivx) + SQR(wl_ivy) + SQR(wl_ivz)) + pbl; + const Real er = wr_idn * wr_ise + + 0.5 * wr_idn * (SQR(wr_ivx) + SQR(wr_ivy) + SQR(wr_ivz)) + pbr; + + const Real fl_d = wl_idn * wl_ivx; + const Real fl_mx = wl_idn * SQR(wl_ivx) + ptl - bxsq_n; + const Real fl_my = wl_idn * wl_ivx * wl_ivy - bxi_n * wl_iby_n; + const Real fl_mz = wl_idn * wl_ivx * wl_ivz - bxi_n * wl_ibz_n; + const Real fl_e = (el + ptl) * wl_ivx - bxi_n * vdotBl; + const Real fl_by = wl_ivx * wl_iby - bxi * wl_ivy; + const Real fl_bz = wl_ivx * wl_ibz - bxi * wl_ivz; + + const Real fr_d = wr_idn * wr_ivx; + const Real fr_mx = wr_idn * SQR(wr_ivx) + ptr - bxsq_n; + const Real fr_my = wr_idn * wr_ivx * wr_ivy - bxi_n * wr_iby_n; + const Real fr_mz = wr_idn * wr_ivx * wr_ivz - bxi_n * wr_ibz_n; + const Real fr_e = (er + ptr) * wr_ivx - bxi_n * vdotBr; + const Real fr_by = wr_ivx * wr_iby - bxi * wr_ivy; + const Real fr_bz = wr_ivx * wr_ibz - bxi * wr_ivz; + + const Real cfl = + MHD::FastMagnetosonicSpeed(wl_ibl, wl_idn, bxi, wl_iby, wl_ibz, mu0); + const Real cfr = + MHD::FastMagnetosonicSpeed(wr_ibl, wr_idn, bxi, wr_iby, wr_ibz, mu0); + + const Real sl = std::min(wl_ivx - cfl, wr_ivx - cfr); + const Real sr = std::max(wl_ivx + cfl, wr_ivx + cfr); + + const Real hlle_qa = sr > 0.0 ? sr : small; + const Real hlle_qb = sl < 0.0 ? sl : -small; + const Real hlle_qc = 0.5 * (hlle_qa + hlle_qb) / (hlle_qa - hlle_qb); + const Real hlle_d = 0.5 * (wl_ipr + wr_ipr) + hlle_qc * (wl_ipr - wr_ipr); + const Real hlle_db = 0.5 * (pbl + pbr) + hlle_qc * (pbl - pbr); + const Real hlle_frho = + 0.5 * ((wl_ivx - hlle_qb) * wl_idn + (wr_ivx - hlle_qa) * wr_idn) + + hlle_qc * ((wl_ivx - hlle_qb) * wl_idn - (wr_ivx - hlle_qa) * wr_idn); + const Real hlle_fmx = + 0.5 * ((wl_ivx - hlle_qb) * wl_idn * wl_ivx - bxsq_n + + (wr_ivx - hlle_qa) * wr_idn * wr_ivx - bxsq_n) + + hlle_qc * ((wl_ivx - hlle_qb) * wl_idn * wl_ivx - bxsq_n - + ((wr_ivx - hlle_qa) * wr_idn * wr_ivx - bxsq_n)); + const Real hlle_fmy = + 0.5 * ((wl_ivx - hlle_qb) * wl_idn * wl_ivy - bxi_n * wl_iby_n + + (wr_ivx - hlle_qa) * wr_idn * wr_ivy - bxi_n * wr_iby_n) + + hlle_qc * ((wl_ivx - hlle_qb) * wl_idn * wl_ivy - bxi_n * wl_iby_n - + ((wr_ivx - hlle_qa) * wr_idn * wr_ivy - bxi_n * wr_iby_n)); + const Real hlle_fmz = + 0.5 * ((wl_ivx - hlle_qb) * wl_idn * wl_ivz - bxi_n * wl_ibz_n + + (wr_ivx - hlle_qa) * wr_idn * wr_ivz - bxi_n * wr_ibz_n) + + hlle_qc * ((wl_ivx - hlle_qb) * wl_idn * wl_ivz - bxi_n * wl_ibz_n - + ((wr_ivx - hlle_qa) * wr_idn * wr_ivz - bxi_n * wr_ibz_n)); + const Real hlle_fe = + 0.5 * (el * (wl_ivx - hlle_qb) + ptl * wl_ivx - bxi_n * vdotBl + + er * (wr_ivx - hlle_qa) + ptr * wr_ivx - bxi_n * vdotBr) + + hlle_qc * (el * (wl_ivx - hlle_qb) + ptl * wl_ivx - bxi_n * vdotBl - + (er * (wr_ivx - hlle_qa) + ptr * wr_ivx - bxi_n * vdotBr)); + const Real hlle_fby = + 0.5 * (fl_by - hlle_qb * wl_iby + fr_by - hlle_qa * wr_iby) + + hlle_qc * ((fl_by - hlle_qb * wl_iby) - (fr_by - hlle_qa * wr_iby)); + const Real hlle_fbz = + 0.5 * (fl_bz - hlle_qb * wl_ibz + fr_bz - hlle_qa * wr_ibz) + + hlle_qc * ((fl_bz - hlle_qb * wl_ibz) - (fr_bz - hlle_qa * wr_ibz)); + + const Real sdl = sl - wl_ivx; + const Real sdr = sr - wr_ivx; + const Real denom = wr_idn * sdr - wl_idn * sdl; + + bool use_hlle = (std::abs(bxi) <= eps) || (std::abs(denom) <= small) || + !std::isfinite(sl) || !std::isfinite(sr); + + Real sm = 0.0; + Real ptst = 0.0; + Real dlst = 0.0; + Real drst = 0.0; + Real sal = 0.0; + Real sar = 0.0; + + Real vlst_y = 0.0, vlst_z = 0.0, blst_y = 0.0, blst_z = 0.0, elst = 0.0; + Real vrst_y = 0.0, vrst_z = 0.0, brst_y = 0.0, brst_z = 0.0, erst = 0.0; + Real vdst_y = 0.0, vdst_z = 0.0, bdst_y = 0.0, bdst_z = 0.0; + Real eldst = 0.0, erdst = 0.0; + + if (!use_hlle) { + sm = (wr_idn * sdr * wr_ivx - wl_idn * sdl * wl_ivx + ptl - ptr) / denom; + const Real ptstl = ptl + wl_idn * sdl * (sm - wl_ivx); + const Real ptstr = ptr + wr_idn * sdr * (sm - wr_ivx); + ptst = 0.5 * (ptstl + ptstr); + + const Real sdml = sl - sm; + const Real sdmr = sr - sm; + if (std::abs(sdml) <= small || std::abs(sdmr) <= small) { + use_hlle = true; + } else { + dlst = wl_idn * sdl / sdml; + drst = wr_idn * sdr / sdmr; + + if (dlst <= 0.0 || drst <= 0.0 || !std::isfinite(dlst) || + !std::isfinite(drst) || ptst <= 0.0 || !std::isfinite(ptst)) { + use_hlle = true; + } else { + const Real dsl = wl_idn * sdl * sdml - bxsq_n; + const Real dsr = wr_idn * sdr * sdmr - bxsq_n; + const Real deg_tol = 1.0e-4 * ptst; + + if (std::abs(dsl) < deg_tol) { + vlst_y = wl_ivy; + vlst_z = wl_ivz; + blst_y = wl_iby; + blst_z = wl_ibz; + } else { + const Real inv_dsl = 1.0 / dsl; + const Real mfact = bxi_n * (sm - wl_ivx) * inv_dsl; + const Real bfact = (wl_idn * SQR(sdl) - bxsq_n) * inv_dsl; + vlst_y = wl_ivy - wl_iby_n * mfact; + vlst_z = wl_ivz - wl_ibz_n * mfact; + blst_y = wl_iby_n * bfact; + blst_z = wl_ibz_n * bfact; + } + + if (std::abs(dsr) < deg_tol) { + vrst_y = wr_ivy; + vrst_z = wr_ivz; + brst_y = wr_iby; + brst_z = wr_ibz; + } else { + const Real inv_dsr = 1.0 / dsr; + const Real mfact = bxi_n * (sm - wr_ivx) * inv_dsr; + const Real bfact = (wr_idn * SQR(sdr) - bxsq_n) * inv_dsr; + vrst_y = wr_ivy - wr_iby_n * mfact; + vrst_z = wr_ivz - wr_ibz_n * mfact; + brst_y = wr_iby_n * bfact; + brst_z = wr_ibz_n * bfact; + } + + const Real vbstl = sm * bxi_n + vlst_y * blst_y + vlst_z * blst_z; + const Real vbstr = sm * bxi_n + vrst_y * brst_y + vrst_z * brst_z; + elst = (sdl * el - ptl * wl_ivx + ptst * sm + bxi_n * (vdotBl - vbstl)) / + sdml; + erst = (sdr * er - ptr * wr_ivx + ptst * sm + bxi_n * (vdotBr - vbstr)) / + sdmr; + + sal = sm - std::abs(bxi_n) / std::sqrt(dlst); + sar = sm + std::abs(bxi_n) / std::sqrt(drst); + + const Real sqrtdlst = std::sqrt(dlst); + const Real sqrtdrst = std::sqrt(drst); + const Real denom_dst = sqrtdlst + sqrtdrst; + if (std::abs(denom_dst) <= small || !std::isfinite(denom_dst)) { + use_hlle = true; + } else { + const Real sgnbx = (bxi_n >= 0.0) ? 1.0 : -1.0; + const Real inv_dst = 1.0 / denom_dst; + + if (0.5 * bxsq_n < deg_tol) { + vdst_y = vlst_y; + vdst_z = vlst_z; + bdst_y = blst_y; + bdst_z = blst_z; + eldst = elst; + erdst = erst; + } else { + vdst_y = (sqrtdlst * vlst_y + sqrtdrst * vrst_y + + (brst_y - blst_y) * sgnbx) * + inv_dst; + vdst_z = (sqrtdlst * vlst_z + sqrtdrst * vrst_z + + (brst_z - blst_z) * sgnbx) * + inv_dst; + bdst_y = (sqrtdlst * brst_y + sqrtdrst * blst_y + + sqrtdlst * sqrtdrst * (vrst_y - vlst_y) * sgnbx) * + inv_dst; + bdst_z = (sqrtdlst * brst_z + sqrtdrst * blst_z + + sqrtdlst * sqrtdrst * (vrst_z - vlst_z) * sgnbx) * + inv_dst; + + const Real vbdst = sm * bxi_n + vdst_y * bdst_y + vdst_z * bdst_z; + eldst = elst - sqrtdlst * + ((sm * bxi_n + vlst_y * blst_y + vlst_z * blst_z) - + vbdst) * + sgnbx; + erdst = erst + sqrtdrst * + ((sm * bxi_n + vrst_y * brst_y + vrst_z * brst_z) - + vbdst) * + sgnbx; + } + + use_hlle = !std::isfinite(sal) || !std::isfinite(sar) || + !std::isfinite(elst) || !std::isfinite(erst) || + !std::isfinite(eldst) || !std::isfinite(erdst); + } + } + } + } + + if (use_hlle) { + frho = hlle_frho; + fmx = hlle_fmx; + fmy = hlle_fmy; + fmz = hlle_fmz; + fe = hlle_fe; + fby = hlle_fby; + fbz = hlle_fbz; + pface = hlle_d; + pmag_face = hlle_db; + } else if (sl >= 0.0) { + frho = fl_d; + fmx = fl_mx; + fmy = fl_my; + fmz = fl_mz; + fe = fl_e; + fby = fl_by; + fbz = fl_bz; + pface = wl_ipr; + pmag_face = pbl; + } else if (sal >= 0.0) { + frho = fl_d + sl * (dlst - wl_idn); + fmx = fl_mx + sl * (dlst * sm - wl_idn * wl_ivx); + fmy = fl_my + sl * (dlst * vlst_y - wl_idn * wl_ivy); + fmz = fl_mz + sl * (dlst * vlst_z - wl_idn * wl_ivz); + fe = fl_e + sl * (elst - el); + fby = fl_by + sl * (sqrt_mu0 * blst_y - wl_iby); + fbz = fl_bz + sl * (sqrt_mu0 * blst_z - wl_ibz); + pface = ptst - 0.5 * (bxsq_n + SQR(blst_y) + SQR(blst_z)); + pmag_face = 0.5 * (bxsq_n + SQR(blst_y) + SQR(blst_z)); + } else if (sm >= 0.0) { + frho = fl_d + sl * (dlst - wl_idn); + fmx = fl_mx + sl * (dlst * sm - wl_idn * wl_ivx); + fmy = fl_my + sl * (dlst * vlst_y - wl_idn * wl_ivy) + + sal * dlst * (vdst_y - vlst_y); + fmz = fl_mz + sl * (dlst * vlst_z - wl_idn * wl_ivz) + + sal * dlst * (vdst_z - vlst_z); + fe = fl_e + sl * (elst - el) + sal * (eldst - elst); + fby = fl_by + sl * (sqrt_mu0 * blst_y - wl_iby) + + sal * sqrt_mu0 * (bdst_y - blst_y); + fbz = fl_bz + sl * (sqrt_mu0 * blst_z - wl_ibz) + + sal * sqrt_mu0 * (bdst_z - blst_z); + pface = ptst - 0.5 * (bxsq_n + SQR(bdst_y) + SQR(bdst_z)); + pmag_face = 0.5 * (bxsq_n + SQR(bdst_y) + SQR(bdst_z)); + } else if (sar > 0.0) { + frho = fr_d + sr * (drst - wr_idn); + fmx = fr_mx + sr * (drst * sm - wr_idn * wr_ivx); + fmy = fr_my + sr * (drst * vrst_y - wr_idn * wr_ivy) + + sar * drst * (vdst_y - vrst_y); + fmz = fr_mz + sr * (drst * vrst_z - wr_idn * wr_ivz) + + sar * drst * (vdst_z - vrst_z); + fe = fr_e + sr * (erst - er) + sar * (erdst - erst); + fby = fr_by + sr * (sqrt_mu0 * brst_y - wr_iby) + + sar * sqrt_mu0 * (bdst_y - brst_y); + fbz = fr_bz + sr * (sqrt_mu0 * brst_z - wr_ibz) + + sar * sqrt_mu0 * (bdst_z - brst_z); + pface = ptst - 0.5 * (bxsq_n + SQR(bdst_y) + SQR(bdst_z)); + pmag_face = 0.5 * (bxsq_n + SQR(bdst_y) + SQR(bdst_z)); + } else if (sr > 0.0) { + frho = fr_d + sr * (drst - wr_idn); + fmx = fr_mx + sr * (drst * sm - wr_idn * wr_ivx); + fmy = fr_my + sr * (drst * vrst_y - wr_idn * wr_ivy); + fmz = fr_mz + sr * (drst * vrst_z - wr_idn * wr_ivz); + fe = fr_e + sr * (erst - er); + fby = fr_by + sr * (sqrt_mu0 * brst_y - wr_iby); + fbz = fr_bz + sr * (sqrt_mu0 * brst_z - wr_ibz); + pface = ptst - 0.5 * (bxsq_n + SQR(brst_y) + SQR(brst_z)); + pmag_face = 0.5 * (bxsq_n + SQR(brst_y) + SQR(brst_z)); + } else { + frho = fr_d; + fmx = fr_mx; + fmy = fr_my; + fmz = fr_mz; + fe = fr_e; + fby = fr_by; + fbz = fr_bz; + pface = wr_ipr; + pmag_face = pbr; + } + + q.flux(b, dir, IDN, k, j, i) = frho; + q.flux(b, dir, ivx, k, j, i) = fmx; + q.flux(b, dir, ivy, k, j, i) = fmy; + q.flux(b, dir, ivz, k, j, i) = fmz; + p.flux(b, dir, IPR, k, j, i) = pface; + p.flux(b, dir, IBM, k, j, i) = pmag_face; + p.flux(b, dir, field::cell::B(IBXG), k, j, i) = 0.0; + p.flux(b, dir, field::cell::B(IBYG), k, j, i) = fby; + p.flux(b, dir, field::cell::B(IBZG), k, j, i) = fbz; + + q.flux(b, dir, IEN, k, j, i) = fe; + + // Li, 2008, https://ui.adsabs.harvard.edu/abs/2008ASPC..385..273L/abstract + q.flux(b, dir, IEG, k, j, i) = frho * ((frho >= 0.0) ? wl_ise : wr_ise); + vf(b, fdir, n, k, j, i) = frho / ((frho >= 0.0) ? wl_idn : wr_idn); + }); + + // other species are just hllc + for (int n = 1; n < nspecies; ++n) { + const int IDN = n; + const int ivx = nspecies + (n * 3) + ((dir - 1)); + const int ivy = nspecies + (n * 3) + ((dir - 1) + 1) % 3; + const int ivz = nspecies + (n * 3) + ((dir - 1) + 2) % 3; + // Unused indices for dust hydrodynamics + const int IPR = nspecies * 4 + n; + const int ISE = nspecies * 5 + n; + const int IBL = nspecies * 6 + n; + const int IEN = IPR; + const int IEG = ISE; + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, member, il, iu, [&](const int i) { + // Create local references for L/R states (helps compiler vectorize) + Real &wl_idn = wl(IDN, i); + Real &wl_ivx = wl(ivx, i); + Real &wl_ivy = wl(ivy, i); + Real &wl_ivz = wl(ivz, i); + Real &wl_ipr = wl(IPR, i); + Real &wl_ise = wl(ISE, i); + Real &wl_ibl = wl(IBL, i); + + Real &wr_idn = wr(IDN, i); + Real &wr_ivx = wr(ivx, i); + Real &wr_ivy = wr(ivy, i); + Real &wr_ivz = wr(ivz, i); + Real &wr_ipr = wr(IPR, i); + Real &wr_ise = wr(ISE, i); + Real &wr_ibl = wr(IBL, i); + + // Compute middle state estimates with PVRS (Toro 10.5.2) + // define 6 registers used below + Real qa, qb, qc, qd, qe, qf; + qa = wl_ibl / wl_idn; + qb = wr_ibl / wr_idn; + Real el = wl_idn * (wl_ise + 0.5 * (SQR(wl_ivx) + SQR(wl_ivy) + SQR(wl_ivz))); + Real er = wr_idn * (wr_ise + 0.5 * (SQR(wr_ivx) + SQR(wr_ivy) + SQR(wr_ivz))); + + // NOTE(@adempsey) + // The below choices are taken from Batten et al 1997 and Fleischmann et al + // 2020 Roe averages + Real sqrtl = std::sqrt(wl_idn); + Real sqrtr = std::sqrt(wr_idn); + const Real isqrt = 1.0 / (sqrtl + sqrtr); + sqrtl *= isqrt; + sqrtr *= isqrt; + const Real vxh = sqrtl * wl_ivx + sqrtr * wr_ivx; + const Real csh = std::sqrt(sqrtl * qa + sqrtr * qb + + 0.5 * sqrtl * sqrtr * SQR(wl_ivx - wr_ivx)); + qa = std::sqrt(qa); + qb = std::sqrt(qb); + + // Compute the max/min wave speeds based on L/R + Real sl = std::min(wl_ivx - qa, vxh - csh); + Real sr = std::max(wr_ivx + qb, vxh + csh); + + // following min/max set to TINY_NUMBER to fix bug found in converging + // supersonic flow + qa = sr > 0.0 ? sr : 1.0e-20; // bp + qb = sl < 0.0 ? sl : -1.0e-20; // bm + + // Compute the contact wave speed and pressure + qe = wl_ivx - sl; // vxl + qf = wr_ivx - sr; // vxr + + qc = wl_ipr + qe * wl_idn * wl_ivx; // tl + qd = wr_ipr + qf * wr_idn * wr_ivx; // tr + + Real ml = wl_idn * qe; + Real mr = -(wr_idn * qf); + + // Determine the contact wave speed... + Real am = (qc - qd) / (ml + mr); + // ...and the pressure at the contact surface + Real cp = (ml * qd + mr * qc) / (ml + mr); + cp = cp > 0.0 ? cp : 0.0; + + // Compute L/R fluxes along the line bm (qb), bp (qa) + qe = wl_idn * (wl_ivx - qb); + qf = wr_idn * (wr_ivx - qa); + + Real fld = qe; + Real frd = qf; + Real flmx = qe * wl_ivx; // + wl_ipr; + Real frmx = qf * wr_ivx; // + wr_ipr; + Real flmy = qe * wl_ivy; + Real frmy = qf * wr_ivy; + Real flmz = qe * wl_ivz; + Real frmz = qf * wr_ivz; + Real fle = el * (wl_ivx - qb) + wl_ipr * wl_ivx; + Real fre = er * (wr_ivx - qa) + wr_ipr * wr_ivx; + + // Compute flux weights or scales. Set an approximate interface pressure for + // coordinate source terms and pressure contribution to flux. + if (am >= 0.0) { + qc = am / (am - qb); + qd = 0.0; + qe = -qb / (am - qb); + } else { + qc = 0.0; + qd = -am / (qa - am); + qe = qa / (qa - am); + } + p.flux(b, dir, IPR, k, j, i) = qc * wl_ipr + qd * wr_ipr + qe * cp; + + // Compute the HLLC flux at interface, including weighted contribution of the + // flux along the contact + const Real frho = qc * fld + qd * frd; + q.flux(b, dir, IDN, k, j, i) = frho; + q.flux(b, dir, ivx, k, j, i) = qc * flmx + qd * frmx; // + qe * cp; + q.flux(b, dir, ivy, k, j, i) = qc * flmy + qd * frmy; + q.flux(b, dir, ivz, k, j, i) = qc * flmz + qd * frmz; + q.flux(b, dir, IEN, k, j, i) = qc * fle + qd * fre + qe * cp * am; + + // Li, 2008, https://ui.adsabs.harvard.edu/abs/2008ASPC..385..273L/abstract + q.flux(b, dir, IEG, k, j, i) = frho * ((frho >= 0.0) ? wl_ise : wr_ise); + vf(b, fdir, n, k, j, i) = frho / ((frho >= 0.0) ? wl_idn : wr_idn); + }); + } + } +}; +} // namespace ArtemisUtils + +#endif // UTILS_FLUXES_RIEMANN_HLLD_HPP_ diff --git a/src/utils/fluxes/riemann/hlle.hpp b/src/utils/fluxes/riemann/hlle.hpp index 0e554747..f1999faa 100644 --- a/src/utils/fluxes/riemann/hlle.hpp +++ b/src/utils/fluxes/riemann/hlle.hpp @@ -53,13 +53,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const EOS &eos, const Real c, const Real chat, const Real mu0, + const bool do_mhd, parthenon::team_mbr_t const &member, const int b, + const int k, const int j, const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); @@ -73,7 +73,16 @@ struct RiemannSolver(); [[maybe_unused]] Real wl_ibl = Null(); [[maybe_unused]] Real wr_ibl = Null(); + [[maybe_unused]] Real wl_ibx = 0.0; + [[maybe_unused]] Real wl_iby = 0.0; + [[maybe_unused]] Real wl_ibz = 0.0; + [[maybe_unused]] Real wr_ibx = 0.0; + [[maybe_unused]] Real wr_iby = 0.0; + [[maybe_unused]] Real wr_ibz = 0.0; if constexpr (FLUID_TYPE == Fluid::gas) { wl_ipr = wl(IPR, i); wl_ise = wl(ISE, i); @@ -111,6 +126,14 @@ struct RiemannSolver(); [[maybe_unused]] Real er = Null(); [[maybe_unused]] Real hroe = Null(); + [[maybe_unused]] Real pbl = 0.0; + [[maybe_unused]] Real pbr = 0.0; + [[maybe_unused]] Real vdBl = 0.0; + [[maybe_unused]] Real vdBr = 0.0; if constexpr (FLUID_TYPE == Fluid::gas) { // Following Roe(1981), the enthalpy H=(E+P)/d is averaged for ideal gas // EOS, rather than E or P directly. sqrtdl*hl = sqrtdl*(el+pl)/dl = @@ -133,14 +160,27 @@ struct RiemannSolver(), qb = Null(); if constexpr (FLUID_TYPE == Fluid::gas) { - qa = std::sqrt(wl_ibl / wl_idn); - qb = std::sqrt(wr_ibl / wr_idn); + if (mhd) { + qa = MHD::FastMagnetosonicSpeed(wl_ibl, wl_idn, wl_ibx, wl_iby, wl_ibz, + mu0); + qb = MHD::FastMagnetosonicSpeed(wr_ibl, wr_idn, wr_ibx, wr_iby, wr_ibz, + mu0); + } else { + qa = std::sqrt(wl_ibl / wl_idn); + qb = std::sqrt(wr_ibl / wr_idn); + } } [[maybe_unused]] Real sl = Null(); @@ -190,9 +230,30 @@ struct RiemannSolver(); [[maybe_unused]] Real fr_e = Null(); + [[maybe_unused]] Real fl_by = 0.0; + [[maybe_unused]] Real fr_by = 0.0; + [[maybe_unused]] Real fl_bz = 0.0; + [[maybe_unused]] Real fr_bz = 0.0; if constexpr (FLUID_TYPE == Fluid::gas) { - fl_e = el * qa + wl_ipr * wl_ivx; - fr_e = er * qb + wr_ipr * wr_ivx; + if (mhd) { + fl_mx -= SQR(wl_ibx) / mu0; + fr_mx -= SQR(wr_ibx) / mu0; + fl_my -= wl_ibx * wl_iby / mu0; + fr_my -= wr_ibx * wr_iby / mu0; + fl_mz -= wl_ibx * wl_ibz / mu0; + fr_mz -= wr_ibx * wr_ibz / mu0; + vdBl = wl_ivx * wl_ibx + wl_ivy * wl_iby + wl_ivz * wl_ibz; + vdBr = wr_ivx * wr_ibx + wr_ivy * wr_iby + wr_ivz * wr_ibz; + fl_e = el * qa + (wl_ipr + pbl) * wl_ivx - wl_ibx * vdBl / mu0; + fr_e = er * qb + (wr_ipr + pbr) * wr_ivx - wr_ibx * vdBr / mu0; + fl_by = (wl_ivx * wl_iby - wl_ivy * wl_ibx) - bm * wl_iby; + fr_by = (wr_ivx * wr_iby - wr_ivy * wr_ibx) - bp * wr_iby; + fl_bz = (wl_ivx * wl_ibz - wl_ivz * wl_ibx) - bm * wl_ibz; + fr_bz = (wr_ivx * wr_ibz - wr_ivz * wr_ibx) - bp * wr_ibz; + } else { + fl_e = el * qa + wl_ipr * wl_ivx; + fr_e = er * qb + wr_ipr * wr_ivx; + } } // Set an approximate interface pressure for coordinate source terms and @@ -202,6 +263,9 @@ struct RiemannSolver struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const EOS &eos, const Real c, const Real chat, const Real mu0, + const bool do_mhd, parthenon::team_mbr_t const &member, const int b, + const int k, const int j, const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); diff --git a/src/utils/fluxes/riemann/llf.hpp b/src/utils/fluxes/riemann/llf.hpp index 0ae2a0fd..f9e8c011 100644 --- a/src/utils/fluxes/riemann/llf.hpp +++ b/src/utils/fluxes/riemann/llf.hpp @@ -44,13 +44,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const EOS &eos, const Real c, const Real chat, const Real mu0, bool do_mhd, + parthenon::team_mbr_t const &member, const int b, const int k, const int j, + const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction @@ -65,7 +65,15 @@ struct RiemannSolver(); [[maybe_unused]] Real wl_ibl = Null(); [[maybe_unused]] Real wr_ibl = Null(); + // note these are intentionally 0 and not Null + [[maybe_unused]] Real wl_ibx = 0.0; + [[maybe_unused]] Real wl_iby = 0.0; + [[maybe_unused]] Real wl_ibz = 0.0; + [[maybe_unused]] Real wr_ibx = 0.0; + [[maybe_unused]] Real wr_iby = 0.0; + [[maybe_unused]] Real wr_ibz = 0.0; if constexpr (FLUID_TYPE == Fluid::gas) { wl_ipr = wl(IPR, i); wl_ise = wl(ISE, i); - wl_ibl = wr(IBL, i); + wl_ibl = wl(IBL, i); wr_ipr = wr(IPR, i); wr_ise = wr(ISE, i); wr_ibl = wr(IBL, i); + if (do_mhd) { + wl_ibx = wl(IBX, i); + wl_iby = wl(IBY, i); + wl_ibz = wl(IBZ, i); + wr_ibx = wr(IBX, i); + wr_iby = wr(IBY, i); + wr_ibz = wr(IBZ, i); + } } // Compute sum of L/R fluxes @@ -116,19 +139,50 @@ struct RiemannSolver(); [[maybe_unused]] Real er = Null(); [[maybe_unused]] Real fsum_e = Null(); + [[maybe_unused]] Real pbl = 0.0; + [[maybe_unused]] Real pbr = 0.0; + [[maybe_unused]] Real vdBl = 0.0; + [[maybe_unused]] Real vdBr = 0.0; + [[maybe_unused]] Real fsum_by = 0.0; + [[maybe_unused]] Real fsum_bz = 0.0; + if constexpr (FLUID_TYPE == Fluid::gas) { el = wl_ise * wl_idn + 0.5 * wl_idn * (SQR(wl_ivx) + SQR(wl_ivy) + SQR(wl_ivz)); er = wr_ise * wr_idn + 0.5 * wr_idn * (SQR(wr_ivx) + SQR(wr_ivy) + SQR(wr_ivz)); fsum_e = (el + wl_ipr) * wl_ivx + (er + wr_ipr) * wr_ivx; + if (do_mhd) { + pbl = MHD::MagneticEnergyDensity(wl_ibx, wl_iby, wl_ibz, mu0); + pbr = MHD::MagneticEnergyDensity(wr_ibx, wr_iby, wr_ibz, mu0); + p.flux(b, dir, IBM, k, j, i) = 0.5 * (pbl + pbr); + fsum_mx -= (SQR(wl_ibx) + SQR(wr_ibx)) / mu0; + fsum_my -= (wl_ibx * wl_iby + wr_ibx * wr_iby) / mu0; + fsum_mz -= (wl_ibx * wl_ibz + wr_ibx * wr_ibz) / mu0; + vdBl = wl_ivx * wl_ibx + wl_ivy * wl_iby + wl_ivz * wl_ibz; + vdBr = wr_ivx * wr_ibx + wr_ivy * wr_iby + wr_ivz * wr_ibz; + el += pbl; + er += pbr; + fsum_e -= (wl_ibx * vdBl + wr_ibx * vdBr) / mu0; + fsum_by = (wl_ivx * wl_iby - wl_ivy * wl_ibx) + + (wr_ivx * wr_iby - wr_ivy * wr_ibx); + fsum_bz = (wl_ivx * wl_ibz - wl_ivz * wl_ibx) + + (wr_ivx * wr_ibz - wr_ivz * wr_ibx); + } } // Compute max wave speed in L/R states (see Toro eq. 10.43) Real a = Null(); if constexpr (FLUID_TYPE == Fluid::gas) { - qa = std::sqrt(wl_ibl / wl_idn); - qb = std::sqrt(wr_ibl / wr_idn); + if (do_mhd) { + qa = MHD::FastMagnetosonicSpeed(wl_ibl, wl_idn, wl_ibx, wl_iby, wl_ibz, + mu0); + qb = MHD::FastMagnetosonicSpeed(wr_ibl, wr_idn, wr_ibx, wr_iby, wr_ibz, + mu0); + } else { + qa = std::sqrt(wl_ibl / wl_idn); + qb = std::sqrt(wr_ibl / wr_idn); + } a = std::max((std::abs(wl_ivx) + qa), (std::abs(wr_ivx) + qb)); } else if constexpr (FLUID_TYPE == Fluid::dust) { a = std::max(std::abs(wl_ivx), std::abs(wr_ivx)); @@ -141,8 +195,14 @@ struct RiemannSolver(); + [[maybe_unused]] Real du_by = 0.0; + [[maybe_unused]] Real du_bz = 0.0; if constexpr (FLUID_TYPE == Fluid::gas) { du_e = a * (er - el); + if (do_mhd) { + du_by = a * (wr_iby - wl_iby); + du_bz = a * (wr_ibz - wl_ibz); + } } // Set an approximate interface pressure for coordinate source terms @@ -157,6 +217,13 @@ struct RiemannSolver struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const EOS &eos, const Real c, const Real chat, const Real mu0, + const bool do_mhd, parthenon::team_mbr_t const &member, const int b, + const int k, const int j, const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); diff --git a/src/utils/fluxes/riemann/riemann.hpp b/src/utils/fluxes/riemann/riemann.hpp index 80c42bde..037ef5e9 100644 --- a/src/utils/fluxes/riemann/riemann.hpp +++ b/src/utils/fluxes/riemann/riemann.hpp @@ -39,6 +39,7 @@ struct RiemannSolver { // Partial specializations #include "hllc.hpp" +#include "hlld.hpp" #include "hlle.hpp" #include "llf.hpp" diff --git a/src/utils/integrators/artemis_integrator.hpp b/src/utils/integrators/artemis_integrator.hpp index c7c9f02c..b4691b88 100644 --- a/src/utils/integrators/artemis_integrator.hpp +++ b/src/utils/integrators/artemis_integrator.hpp @@ -32,13 +32,17 @@ inline TaskStatus DeepCopyConservedData(MeshData *to, MeshData *from using parthenon::MakePackDescriptor; using parthenon::variable_names::any; - std::vector flags({Metadata::Conserved}); + std::vector flags({Metadata::Cell, Metadata::Conserved}); static auto desc = MakePackDescriptor(to, flags); + std::vector flags_f({Metadata::Face, Metadata::Conserved}); + static auto desc_b = MakePackDescriptor(to, flags_f); const auto vt = desc.GetPack(to); const auto vf = desc.GetPack(from); - const auto ibe = to->GetBoundsI(IndexDomain::entire); - const auto jbe = to->GetBoundsJ(IndexDomain::entire); - const auto kbe = to->GetBoundsK(IndexDomain::entire); + const auto vtf = desc_b.GetPack(to); + const auto vff = desc_b.GetPack(from); + auto ibe = to->GetBoundsI(IndexDomain::entire); + auto jbe = to->GetBoundsJ(IndexDomain::entire); + auto kbe = to->GetBoundsK(IndexDomain::entire); parthenon::par_for( DEFAULT_LOOP_PATTERN, "DeepCopyConservedData", parthenon::DevExecSpace(), 0, @@ -48,6 +52,40 @@ inline TaskStatus DeepCopyConservedData(MeshData *to, MeshData *from vt(b, n, k, j, i) = vf(b, n, k, j, i); } }); + + ibe = to->GetBoundsI(IndexDomain::entire, TE::F1); + jbe = to->GetBoundsJ(IndexDomain::entire, TE::F1); + kbe = to->GetBoundsK(IndexDomain::entire, TE::F1); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "DeepCopyConservedData::F1", parthenon::DevExecSpace(), 0, + to->NumBlocks() - 1, kbe.s, kbe.e, jbe.s, jbe.e, ibe.s, ibe.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + for (int n = vtf.GetLowerBound(b); n <= vtf.GetUpperBound(b); ++n) { + vtf(b, TE::F1, n, k, j, i) = vff(b, TE::F1, n, k, j, i); + } + }); + ibe = to->GetBoundsI(IndexDomain::entire, TE::F2); + jbe = to->GetBoundsJ(IndexDomain::entire, TE::F2); + kbe = to->GetBoundsK(IndexDomain::entire, TE::F2); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "DeepCopyConservedData::F2", parthenon::DevExecSpace(), 0, + to->NumBlocks() - 1, kbe.s, kbe.e, jbe.s, jbe.e, ibe.s, ibe.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + for (int n = vtf.GetLowerBound(b); n <= vtf.GetUpperBound(b); ++n) { + vtf(b, TE::F2, n, k, j, i) = vff(b, TE::F2, n, k, j, i); + } + }); + ibe = to->GetBoundsI(IndexDomain::entire, TE::F3); + jbe = to->GetBoundsJ(IndexDomain::entire, TE::F3); + kbe = to->GetBoundsK(IndexDomain::entire, TE::F3); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "DeepCopyConservedData::F3", parthenon::DevExecSpace(), 0, + to->NumBlocks() - 1, kbe.s, kbe.e, jbe.s, jbe.e, ibe.s, ibe.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + for (int n = vtf.GetLowerBound(b); n <= vtf.GetUpperBound(b); ++n) { + vtf(b, TE::F3, n, k, j, i) = vff(b, TE::F3, n, k, j, i); + } + }); return TaskStatus::complete; } @@ -63,12 +101,14 @@ TaskStatus ApplyUpdate(MeshData *u0, MeshData *u1, const Real g0, auto pm = u0->GetParentPointer(); // Packing and indexing - std::vector flags({Metadata::Conserved}); + std::vector flags({Metadata::Conserved, Metadata::Cell}); static auto desc = MakePackDescriptor(u0, flags, {parthenon::PDOpt::WithFluxes}); - static auto desc_g = MakePackDescriptor(u0); + static auto desc_g = MakePackDescriptor(u0); const auto v0 = desc.GetPack(u0); const auto v1 = desc.GetPack(u1); const auto vg = desc_g.GetPack(u1); + const auto ib = u0->GetBoundsI(IndexDomain::interior); const auto jb = u0->GetBoundsJ(IndexDomain::interior); const auto kb = u0->GetBoundsK(IndexDomain::interior); @@ -111,6 +151,113 @@ TaskStatus ApplyUpdate(MeshData *u0, MeshData *u1, const Real g0, return TaskStatus::complete; } +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus ArtemisUtils::ApplyFaceUpdate +//! \brief +template +TaskStatus ApplyFaceUpdate(MeshData *u0, MeshData *u1, const Real g0, + const Real g1, const Real beta_dt) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + using parthenon::variable_names::any; + auto pm = u0->GetParentPointer(); + + // Packing and indexing + static auto desc = + MakePackDescriptor(u0, {}, {parthenon::PDOpt::WithFluxes}); + static auto desc_g = MakePackDescriptor(u0); + const auto v0 = desc.GetPack(u0); + const auto v1 = desc.GetPack(u1); + const auto vg = desc_g.GetPack(u1); + + const auto ib = u0->GetBoundsI(IndexDomain::interior); + const auto jb = u0->GetBoundsJ(IndexDomain::interior); + const auto kb = u0->GetBoundsK(IndexDomain::interior); + + const bool multi_d = (pm->ndim > 1); + const bool three_d = (pm->ndim > 2); + const auto &cpars = + pm->packages.Get("artemis")->template Param("coord_params"); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ApplyFaceUpdate::X1", parthenon::DevExecSpace(), 0, + u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e + 1, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, v0.GetCoordinates(b), k, j, i); + geometry::Coords coords_kp(cpars, v0.GetCoordinates(b), k + three_d, j, i); + geometry::Coords coords_jp(cpars, v0.GetCoordinates(b), k, j + multi_d, i); + Real &v0n = v0(b, TE::F1, field::face::B(), k, j, i); + Real &v1n = v1(b, TE::F1, field::face::B(), k, j, i); + v0n = g0 * v0n + g1 * v1n; + + const Real bdt = beta_dt / coords.GetFaceAreaX1()[0]; + const auto &bnds = coords.GetBounds(); + const auto &bnds_kp = coords_kp.GetBounds(); + const auto &bnds_jp = coords_jp.GetBounds(); + const Real dq2m = bnds.x2[1] - bnds.x2[0]; + const Real dq2p = bnds_kp.x2[1] - bnds_kp.x2[0]; + const Real dq3m = bnds.x3[1] - bnds.x3[0]; + const Real dq3p = bnds_jp.x3[1] - bnds_jp.x3[0]; + + v0n -= bdt * ((dq3p * v0.flux(b, TE::E3, field::face::B(), k, j + multi_d, i) - + dq3m * v0.flux(b, TE::E3, field::face::B(), k, j, i)) + + (dq2m * v0.flux(b, TE::E2, field::face::B(), k, j, i) - + dq2p * v0.flux(b, TE::E2, field::face::B(), k + three_d, j, i))); + }); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ApplyFaceUpdate::X2", parthenon::DevExecSpace(), 0, + u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e + multi_d, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, v0.GetCoordinates(b), k, j, i); + geometry::Coords coords_kp(cpars, v0.GetCoordinates(b), k + three_d, j, i); + geometry::Coords coords_ip(cpars, v0.GetCoordinates(b), k, j, i + 1); + Real &v0n = v0(b, TE::F2, field::face::B(), k, j, i); + Real &v1n = v1(b, TE::F2, field::face::B(), k, j, i); + v0n = g0 * v0n + g1 * v1n; + + const Real bdt = beta_dt / coords.GetFaceAreaX2()[0]; + const auto &bnds = coords.GetBounds(); + const auto &bnds_kp = coords_kp.GetBounds(); + const auto &bnds_ip = coords_ip.GetBounds(); + const Real dq1m = bnds.x1[1] - bnds.x1[0]; + const Real dq1p = bnds_kp.x1[1] - bnds_kp.x1[0]; + const Real dq3m = bnds.x3[1] - bnds.x3[0]; + const Real dq3p = bnds_ip.x3[1] - bnds_ip.x3[0]; + + v0n -= bdt * ((dq3m * v0.flux(b, TE::E3, field::face::B(), k, j, i) - + dq3p * v0.flux(b, TE::E3, field::face::B(), k, j, i + 1)) + + (dq1p * v0.flux(b, TE::E1, field::face::B(), k + three_d, j, i) - + dq1m * v0.flux(b, TE::E1, field::face::B(), k, j, i))); + }); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ApplyFaceUpdate::X3", parthenon::DevExecSpace(), 0, + u0->NumBlocks() - 1, kb.s, kb.e + three_d, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, v0.GetCoordinates(b), k, j, i); + geometry::Coords coords_jp(cpars, v0.GetCoordinates(b), k, j + multi_d, i); + geometry::Coords coords_ip(cpars, v0.GetCoordinates(b), k, j, i + 1); + Real &v0n = v0(b, TE::F3, field::face::B(), k, j, i); + Real &v1n = v1(b, TE::F3, field::face::B(), k, j, i); + v0n = g0 * v0n + g1 * v1n; + + const Real bdt = beta_dt / coords.GetFaceAreaX3()[0]; + const auto &bnds = coords.GetBounds(); + const auto &bnds_jp = coords_jp.GetBounds(); + const auto &bnds_ip = coords_ip.GetBounds(); + const Real dq1m = bnds.x1[1] - bnds.x1[0]; + const Real dq1p = bnds_jp.x1[1] - bnds_jp.x1[0]; + const Real dq2m = bnds.x2[1] - bnds.x2[0]; + const Real dq2p = bnds_ip.x2[1] - bnds_ip.x2[0]; + + v0n -= bdt * ((dq2p * v0.flux(b, TE::E2, field::face::B(), k, j, i + 1) - + dq2m * v0.flux(b, TE::E2, field::face::B(), k, j, i)) + + (dq1m * v0.flux(b, TE::E1, field::face::B(), k, j, i) - + dq1p * v0.flux(b, TE::E1, field::face::B(), k, j + multi_d, i))); + }); + + return TaskStatus::complete; +} + } // namespace ArtemisUtils #endif // UTILS_INTEGRATORS_ARTEMIS_INTEGRATOR_HPP_ diff --git a/src/utils/refinement/prolongation.hpp b/src/utils/refinement/prolongation.hpp index a21c488e..9a7a1612 100644 --- a/src/utils/refinement/prolongation.hpp +++ b/src/utils/refinement/prolongation.hpp @@ -32,11 +32,69 @@ #include "geometry/geometry.hpp" namespace ArtemisUtils { +template +KOKKOS_FORCEINLINE_FUNCTION Real +GetElementAverageWeight(const geometry::Coords &coords) { + if constexpr (EL == TE::CC) { + return coords.Volume(); + } else if constexpr (EL == TE::F1) { + return coords.template GetFaceArea(); + } else if constexpr (EL == TE::F2) { + return coords.template GetFaceArea(); + } else if constexpr (EL == TE::F3) { + return coords.template GetFaceArea(); + } else if constexpr (EL == TE::E1) { + return coords.GetEdgeLengthX1(); + } else if constexpr (EL == TE::E2) { + return coords.GetEdgeLengthX2(); + } else if constexpr (EL == TE::E3) { + return coords.GetEdgeLengthX3(); + } + PARTHENON_FAIL("Unsupported topological element for geometric averaging weight!"); + return 0.0; +} + +template +KOKKOS_FORCEINLINE_FUNCTION Real GetFaceAverageWeight(const geometry::Coords &coords, + const int face_idx) { + if (face_idx == 0) { + return coords.template GetFaceArea(); + } else if (face_idx == 1) { + return coords.template GetFaceArea(); + } + return coords.template GetFaceArea(); +} + +template +KOKKOS_FORCEINLINE_FUNCTION Real GetCoordinate(const geometry::Coords &coords) { + static_assert(DIM >= 1 && DIM <= 3, "Invalid dimension!"); + if constexpr (EL == TE::CC) { + if constexpr (DIM == 1) { + return coords.x1v(); + } else if constexpr (DIM == 2) { + return coords.x2v(); + } + return coords.x3v(); + } else if constexpr (EL == TE::F1) { + const auto xf = coords.FaceCenX1(geometry::CellFace::lower); + return xf[DIM - 1]; + } else if constexpr (EL == TE::F2) { + const auto xf = coords.FaceCenX2(geometry::CellFace::lower); + return xf[DIM - 1]; + } else if constexpr (EL == TE::F3) { + const auto xf = coords.FaceCenX3(geometry::CellFace::lower); + return xf[DIM - 1]; + } + PARTHENON_FAIL( + "Artemis prolongation only supports cell-centered and face-centered fields!"); + return 0.0; +} + //---------------------------------------------------------------------------------------- //! \fn void ArtemisUtils::GetGridSpacings //! \brief compute distances from cell center to the nearest center in the + or - //! coordinate direction. Do so for both coarse and fine grids. -template +template KOKKOS_FORCEINLINE_FUNCTION void GetGridSpacings(const Coordinates_t &coords, const Coordinates_t &coarse_coords, const bool log, int k, int j, int i, int fk, int fj, int fi, Real *dxm, @@ -49,16 +107,11 @@ GetGridSpacings(const Coordinates_t &coords, const Coordinates_t &coarse_coords, gg::Coords cp(log, coarse_coords, k + (DIM == 3), j + (DIM == 2), i + (DIM == 1)); gg::Coords fm(log, coords, fk, fj, fi); gg::Coords fp(log, coords, fk + (DIM == 3), fj + (DIM == 2), fi + (DIM == 1)); - if constexpr (DIM == 1) { - xm = cm.x1v(), xc = cc.x1v(), xp = cp.x1v(); - fxm = fm.x1v(), fxp = fp.x1v(); - } else if constexpr (DIM == 2) { - xm = cm.x2v(), xc = cc.x2v(), xp = cp.x2v(); - fxm = fm.x2v(), fxp = fp.x2v(); - } else if constexpr (DIM == 3) { - xm = cm.x3v(), xc = cc.x3v(), xp = cp.x3v(); - fxm = fm.x3v(), fxp = fp.x3v(); - } + xm = ArtemisUtils::GetCoordinate(cm); + xc = ArtemisUtils::GetCoordinate(cc); + xp = ArtemisUtils::GetCoordinate(cp); + fxm = ArtemisUtils::GetCoordinate(fm); + fxp = ArtemisUtils::GetCoordinate(fp); *dxm = xc - xm; *dxp = xp - xc; *dxfm = xc - fxm; @@ -95,7 +148,9 @@ struct ProlongateShared { const Coordinates_t &coords, const Coordinates_t &coarse_coords, const ParArrayND *pcoarse, const ParArrayND *pfine) { - PARTHENON_REQUIRE((el == TE::CC), "Artemis AMR only supports cell-centered fields!"); + PARTHENON_REQUIRE( + el == TE::CC || el == TE::F1 || el == TE::F2 || el == TE::F3, + "Artemis AMR only supports cell-centered and face-centered fields!"); auto &coarse = *pcoarse; auto &fine = *pfine; @@ -106,23 +161,33 @@ struct ProlongateShared { const int fj = (DIM > 1) ? (j - cjb.s) * 2 + jb.s : jb.s; const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s; - constexpr bool INCLUDE_X1 = (DIM > 0); - constexpr bool INCLUDE_X2 = (DIM > 1); - constexpr bool INCLUDE_X3 = (DIM > 2); + constexpr bool INCLUDE_X1 = + (DIM > 0) && (el == TE::CC || el == TE::F2 || el == TE::F3); + constexpr bool INCLUDE_X2 = + (DIM > 1) && (el == TE::CC || el == TE::F3 || el == TE::F1); + constexpr bool INCLUDE_X3 = + (DIM > 2) && (el == TE::CC || el == TE::F1 || el == TE::F2); const Real fc = coarse(element_idx, l, m, n, k, j, i); + geometry::Coords cc(log, coarse_coords, k, j, i); + const Real wc = 1.0; + const Real qc = wc * fc; Real dx1fm = 0; [[maybe_unused]] Real dx1fp = 0; [[maybe_unused]] Real gx1m = 0, gx1p = 0; if constexpr (INCLUDE_X1) { Real dx1m, dx1p; - ArtemisUtils::GetGridSpacings(coords, coarse_coords, log, k, j, i, fk, fj, - fi, &dx1m, &dx1p, &dx1fm, &dx1fp); + ArtemisUtils::GetGridSpacings(coords, coarse_coords, log, k, j, i, fk, + fj, fi, &dx1m, &dx1p, &dx1fm, &dx1fp); - Real gx1c = ArtemisUtils::GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1), - coarse(element_idx, l, m, n, k, j, i + 1), - dx1m, dx1p, gx1m, gx1p); + geometry::Coords cm(log, coarse_coords, k, j, i - 1); + geometry::Coords cp(log, coarse_coords, k, j, i + 1); + const Real fm = coarse(element_idx, l, m, n, k, j, i - 1); + const Real fp = coarse(element_idx, l, m, n, k, j, i + 1); + const Real qm = fm; + const Real qp = fp; + Real gx1c = ArtemisUtils::GradMinMod(qc, qm, qp, dx1m, dx1p, gx1m, gx1p); if constexpr (use_minmod_slope) { gx1m = gx1c; gx1p = gx1c; @@ -134,11 +199,15 @@ struct ProlongateShared { [[maybe_unused]] Real gx2m = 0, gx2p = 0; if constexpr (INCLUDE_X2) { Real dx2m, dx2p; - ArtemisUtils::GetGridSpacings(coords, coarse_coords, log, k, j, i, fk, fj, - fi, &dx2m, &dx2p, &dx2fm, &dx2fp); - Real gx2c = ArtemisUtils::GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i), - coarse(element_idx, l, m, n, k, j + 1, i), - dx2m, dx2p, gx2m, gx2p); + ArtemisUtils::GetGridSpacings(coords, coarse_coords, log, k, j, i, fk, + fj, fi, &dx2m, &dx2p, &dx2fm, &dx2fp); + geometry::Coords cm(log, coarse_coords, k, j - 1, i); + geometry::Coords cp(log, coarse_coords, k, j + 1, i); + const Real fm = coarse(element_idx, l, m, n, k, j - 1, i); + const Real fp = coarse(element_idx, l, m, n, k, j + 1, i); + const Real qm = fm; + const Real qp = fp; + Real gx2c = ArtemisUtils::GradMinMod(qc, qm, qp, dx2m, dx2p, gx2m, gx2p); if constexpr (use_minmod_slope) { gx2m = gx2c; gx2p = gx2c; @@ -150,42 +219,360 @@ struct ProlongateShared { [[maybe_unused]] Real gx3m = 0, gx3p = 0; if constexpr (INCLUDE_X3) { Real dx3m, dx3p; - ArtemisUtils::GetGridSpacings(coords, coarse_coords, log, k, j, i, fk, fj, - fi, &dx3m, &dx3p, &dx3fm, &dx3fp); - Real gx3c = ArtemisUtils::GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i), - coarse(element_idx, l, m, n, k + 1, j, i), - dx3m, dx3p, gx3m, gx3p); + ArtemisUtils::GetGridSpacings(coords, coarse_coords, log, k, j, i, fk, + fj, fi, &dx3m, &dx3p, &dx3fm, &dx3fp); + geometry::Coords cm(log, coarse_coords, k - 1, j, i); + geometry::Coords cp(log, coarse_coords, k + 1, j, i); + const Real fm = coarse(element_idx, l, m, n, k - 1, j, i); + const Real fp = coarse(element_idx, l, m, n, k + 1, j, i); + const Real qm = fm; + const Real qp = fp; + Real gx3c = ArtemisUtils::GradMinMod(qc, qm, qp, dx3m, dx3p, gx3m, gx3p); if constexpr (use_minmod_slope) { gx3m = gx3c; gx3p = gx3c; } } - // KGF: add the off-centered quantities first to preserve FP symmetry - // JMM: Extraneous quantities are zero - fine(element_idx, l, m, n, fk, fj, fi) = - fc - (gx1m * dx1fm + gx2m * dx2fm + gx3m * dx3fm); + Real qfine[2][2][2] = {{{0.0}}}; + bool active[2][2][2] = {{{false}}}; + auto stage_fine_value = [&](const int ok, const int oj, const int oi, const Real qf) { + qfine[ok][oj][oi] = qf; + active[ok][oj][oi] = true; + }; + + stage_fine_value(0, 0, 0, qc - (gx1m * dx1fm + gx2m * dx2fm + gx3m * dx3fm)); if constexpr (INCLUDE_X1) - fine(element_idx, l, m, n, fk, fj, fi + 1) = - fc + (gx1p * dx1fp - gx2m * dx2fm - gx3m * dx3fm); + stage_fine_value(0, 0, 1, qc + (gx1p * dx1fp - gx2m * dx2fm - gx3m * dx3fm)); if constexpr (INCLUDE_X2) - fine(element_idx, l, m, n, fk, fj + 1, fi) = - fc - (gx1m * dx1fm - gx2p * dx2fp + gx3m * dx3fm); + stage_fine_value(0, 1, 0, qc - (gx1m * dx1fm - gx2p * dx2fp + gx3m * dx3fm)); if constexpr (INCLUDE_X2 && INCLUDE_X1) - fine(element_idx, l, m, n, fk, fj + 1, fi + 1) = - fc + (gx1p * dx1fp + gx2p * dx2fp - gx3m * dx3fm); + stage_fine_value(0, 1, 1, qc + (gx1p * dx1fp + gx2p * dx2fp - gx3m * dx3fm)); if constexpr (INCLUDE_X3) - fine(element_idx, l, m, n, fk + 1, fj, fi) = - fc - (gx1m * dx1fm + gx2m * dx2fm - gx3p * dx3fp); + stage_fine_value(1, 0, 0, qc - (gx1m * dx1fm + gx2m * dx2fm - gx3p * dx3fp)); if constexpr (INCLUDE_X3 && INCLUDE_X1) - fine(element_idx, l, m, n, fk + 1, fj, fi + 1) = - fc + (gx1p * dx1fp - gx2m * dx2fm + gx3p * dx3fp); + stage_fine_value(1, 0, 1, qc + (gx1p * dx1fp - gx2m * dx2fm + gx3p * dx3fp)); if constexpr (INCLUDE_X3 && INCLUDE_X2) - fine(element_idx, l, m, n, fk + 1, fj + 1, fi) = - fc - (gx1m * dx1fm - gx2p * dx2fp - gx3p * dx3fp); + stage_fine_value(1, 1, 0, qc - (gx1m * dx1fm - gx2p * dx2fp - gx3p * dx3fp)); if constexpr (INCLUDE_X3 && INCLUDE_X2 && INCLUDE_X1) - fine(element_idx, l, m, n, fk + 1, fj + 1, fi + 1) = - fc + (gx1p * dx1fp + gx2p * dx2fp + gx3p * dx3fp); + stage_fine_value(1, 1, 1, qc + (gx1p * dx1fp + gx2p * dx2fp + gx3p * dx3fp)); + + if constexpr (el == TE::F1 || el == TE::F2 || el == TE::F3) { + const Real coarse_flux = GetElementAverageWeight(cc) * fc; + Real fine_flux = 0.0; + Real fine_area = 0.0; + for (int ok = 0; ok < 2; ++ok) { + for (int oj = 0; oj < 2; ++oj) { + for (int oi = 0; oi < 2; ++oi) { + if (!active[ok][oj][oi]) continue; + geometry::Coords cf(log, coords, fk + ok, fj + oj, fi + oi); + const Real area = GetElementAverageWeight(cf); + fine_flux += area * qfine[ok][oj][oi]; + fine_area += area; + } + } + } + const Real delta = (coarse_flux - fine_flux) / (fine_area + Fuzz()); + for (int ok = 0; ok < 2; ++ok) { + for (int oj = 0; oj < 2; ++oj) { + for (int oi = 0; oi < 2; ++oi) { + if (active[ok][oj][oi]) qfine[ok][oj][oi] += delta; + } + } + } + } + + for (int ok = 0; ok < 2; ++ok) { + for (int oj = 0; oj < 2; ++oj) { + for (int oi = 0; oi < 2; ++oi) { + if (active[ok][oj][oi]) { + fine(element_idx, l, m, n, fk + ok, fj + oj, fi + oi) = qfine[ok][oj][oi]; + } + } + } + } + } +}; + +template +struct ProlongateTothAndRoe { + static constexpr bool OperationRequired(TopologicalElement fel, + TopologicalElement cel) { + return (cel == TE::CC) && (GetTopologicalType(fel) == TopologicalType::Face); + } + // Here, fel is the topological element on which the field is defined and + // cel is the topological element on which we are filling the internal values + // of the field. So, for instance, we could fill the fine cell values of an + // x-face field within the volume of a coarse cell. This is assumes that the + // values of the fine cells on the elements corresponding with the coarse cell + // have been filled. + template + KOKKOS_FORCEINLINE_FUNCTION static void + Do(const int l, const int m, const int n, const int k, const int j, const int i, + const IndexRange &ckb, const IndexRange &cjb, const IndexRange &cib, + const IndexRange &kb, const IndexRange &jb, const IndexRange &ib, + const Coordinates_t &coords, const Coordinates_t &coarse_coords, + const ParArrayND *, + const ParArrayND *pfine) { + if constexpr (!IsSubmanifold(fel, cel)) { + return; + } else { + if constexpr (!(fel == TE::F1)) { + return; + } else { + const int fi = (DIM > 0) ? (i - cib.s) * 2 + ib.s : ib.s; + const int fj = (DIM > 1) ? (j - cjb.s) * 2 + jb.s : jb.s; + const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s; + auto &fine = *pfine; + constexpr int g3 = (DIM > 2) ? 1 : 0; + constexpr int g2 = (DIM > 1) ? 1 : 0; + constexpr int ny = (DIM > 1) ? 2 : 1; + constexpr int nz = (DIM > 2) ? 2 : 1; + constexpr int ncell = (DIM == 1 ? 2 : (DIM == 2 ? 4 : 8)); + constexpr int neq = ncell - 1; + constexpr int nf1 = ny * nz; + constexpr int nf2 = (DIM > 1) ? (2 * nz) : 0; + constexpr int nf3 = (DIM > 2) ? 4 : 0; + constexpr int nunk = nf1 + nf2 + nf3; + constexpr int max_eq = 7; + constexpr int max_unk = 12; + + auto get_indices = [&](const int comp, const int sx, const int sy, const int sz) { + return std::array{fk + sz * g3, fj + sy * g2, fi + sx}; + }; + + auto get_face_value = [&](const int comp, const int sx, const int sy, + const int sz) -> Real & { + const auto idx = get_indices(comp, sx, sy, sz); + return fine(comp, l, m, n, idx[0], idx[1], idx[2]); + }; + + auto get_face_area = [&](const int comp, const int sx, const int sy, + const int sz) { + const auto idx = get_indices(comp, sx, sy, sz); + geometry::Coords cf(log, coords, idx[0], idx[1], idx[2]); + return ArtemisUtils::GetFaceAverageWeight(cf, comp); + }; + + auto get_face_flux = [&](const int comp, const int sx, const int sy, + const int sz) { + return get_face_area(comp, sx, sy, sz) * get_face_value(comp, sx, sy, sz); + }; + + auto set_face_flux = [&](const int comp, const int sx, const int sy, const int sz, + const Real qf) { + const Real wf = get_face_area(comp, sx, sy, sz); + get_face_value(comp, sx, sy, sz) = qf / (wf + Fuzz()); + }; + + auto get_face_normal_pos = [&](const int comp, const int sx, const int sy, + const int sz) { + const auto idx = get_indices(comp, sx, sy, sz); + geometry::Coords cf(log, coords, idx[0], idx[1], idx[2]); + if (comp == 0) { + return ArtemisUtils::GetCoordinate(cf); + } else if (comp == 1) { + return ArtemisUtils::GetCoordinate(cf); + } + return ArtemisUtils::GetCoordinate(cf); + }; + + auto idx_f1 = [&](const int sy, const int sz) { return sz * ny + sy; }; + auto idx_f2 = [&](const int sx, const int sz) { return nf1 + sz * 2 + sx; }; + auto idx_f3 = [&](const int sx, const int sy) { return nf1 + nf2 + sy * 2 + sx; }; + + if constexpr (DIM == 1) { + const Real q0 = get_face_flux(0, 0, 0, 0); + const Real q2 = get_face_flux(0, 2, 0, 0); + const Real x0 = get_face_normal_pos(0, 0, 0, 0); + const Real x1 = get_face_normal_pos(0, 1, 0, 0); + const Real x2 = get_face_normal_pos(0, 2, 0, 0); + const Real alpha = (x1 - x0) / (x2 - x0 + Fuzz()); + set_face_flux(0, 1, 0, 0, q0 + alpha * (q2 - q0)); + return; + } + + Real u0[max_unk] = {0.0}; + Real w[max_unk] = {0.0}; + Real D[max_eq][max_unk] = {{0.0}}; + Real rhs[max_eq] = {0.0}; + + for (int sz = 0; sz < nz; ++sz) { + for (int sy = 0; sy < ny; ++sy) { + const int idx = idx_f1(sy, sz); + const Real q0 = get_face_flux(0, 0, sy, sz); + const Real q2 = get_face_flux(0, 2, sy, sz); + const Real x0 = get_face_normal_pos(0, 0, sy, sz); + const Real x1 = get_face_normal_pos(0, 1, sy, sz); + const Real x2 = get_face_normal_pos(0, 2, sy, sz); + const Real alpha = (x1 - x0) / (x2 - x0 + Fuzz()); + u0[idx] = q0 + alpha * (q2 - q0); + w[idx] = get_face_area(0, 1, sy, sz); + } + } + if constexpr (DIM > 1) { + for (int sz = 0; sz < nz; ++sz) { + for (int sx = 0; sx < 2; ++sx) { + const int idx = idx_f2(sx, sz); + const Real q0 = get_face_flux(1, sx, 0, sz); + const Real q2 = get_face_flux(1, sx, 2, sz); + const Real x0 = get_face_normal_pos(1, sx, 0, sz); + const Real x1 = get_face_normal_pos(1, sx, 1, sz); + const Real x2 = get_face_normal_pos(1, sx, 2, sz); + const Real alpha = (x1 - x0) / (x2 - x0 + Fuzz()); + u0[idx] = q0 + alpha * (q2 - q0); + w[idx] = get_face_area(1, sx, 1, sz); + } + } + } + if constexpr (DIM > 2) { + for (int sy = 0; sy < 2; ++sy) { + for (int sx = 0; sx < 2; ++sx) { + const int idx = idx_f3(sx, sy); + const Real q0 = get_face_flux(2, sx, sy, 0); + const Real q2 = get_face_flux(2, sx, sy, 2); + const Real x0 = get_face_normal_pos(2, sx, sy, 0); + const Real x1 = get_face_normal_pos(2, sx, sy, 1); + const Real x2 = get_face_normal_pos(2, sx, sy, 2); + const Real alpha = (x1 - x0) / (x2 - x0 + Fuzz()); + u0[idx] = q0 + alpha * (q2 - q0); + w[idx] = get_face_area(2, sx, sy, 1); + } + } + } + + int row = 0; + for (int sz = 0; sz < nz; ++sz) { + for (int sy = 0; sy < ny; ++sy) { + for (int sx = 0; sx < 2; ++sx) { + if (sx == 1 && sy == ny - 1 && sz == nz - 1) continue; + + if (sx == 0) { + D[row][idx_f1(sy, sz)] += 1.0; + rhs[row] += get_face_flux(0, 0, sy, sz); + } else { + D[row][idx_f1(sy, sz)] -= 1.0; + rhs[row] -= get_face_flux(0, 2, sy, sz); + } + + if constexpr (DIM > 1) { + if (sy == 0) { + D[row][idx_f2(sx, sz)] += 1.0; + rhs[row] += get_face_flux(1, sx, 0, sz); + } else { + D[row][idx_f2(sx, sz)] -= 1.0; + rhs[row] -= get_face_flux(1, sx, 2, sz); + } + } + + if constexpr (DIM > 2) { + if (sz == 0) { + D[row][idx_f3(sx, sy)] += 1.0; + rhs[row] += get_face_flux(2, sx, sy, 0); + } else { + D[row][idx_f3(sx, sy)] -= 1.0; + rhs[row] -= get_face_flux(2, sx, sy, 2); + } + } + + ++row; + } + } + } + + Real M[max_eq][max_eq] = {{0.0}}; + Real residual[max_eq] = {0.0}; + for (int r = 0; r < neq; ++r) { + Real du = 0.0; + for (int u = 0; u < nunk; ++u) { + du += D[r][u] * u0[u]; + } + residual[r] = du - rhs[r]; + for (int c = 0; c < neq; ++c) { + Real sum = 0.0; + for (int u = 0; u < nunk; ++u) { + sum += D[r][u] * w[u] * D[c][u]; + } + M[r][c] = sum; + } + } + + Real A[max_eq][max_eq + 1] = {{0.0}}; + for (int r = 0; r < neq; ++r) { + for (int c = 0; c < neq; ++c) { + A[r][c] = M[r][c]; + } + A[r][neq] = residual[r]; + } + + for (int p = 0; p < neq; ++p) { + int piv = p; + Real max_abs = std::abs(A[p][p]); + for (int r = p + 1; r < neq; ++r) { + const Real cand = std::abs(A[r][p]); + if (cand > max_abs) { + max_abs = cand; + piv = r; + } + } + if (piv != p) { + for (int c = p; c <= neq; ++c) { + const Real tmp = A[p][c]; + A[piv][c] = A[p][c]; + A[p][c] = tmp; + } + } + const Real pivot = A[p][p]; + if (std::abs(pivot) <= Fuzz()) continue; + for (int r = p + 1; r < neq; ++r) { + const Real fac = A[r][p] / pivot; + for (int c = p; c <= neq; ++c) { + A[r][c] -= fac * A[p][c]; + } + } + } + + Real lambda[max_eq] = {0.0}; + for (int r = neq - 1; r >= 0; --r) { + Real sum = A[r][neq]; + for (int c = r + 1; c < neq; ++c) { + sum -= A[r][c] * lambda[c]; + } + const Real pivot = A[r][r]; + lambda[r] = (std::abs(pivot) > Fuzz()) ? (sum / pivot) : 0.0; + } + + Real u[max_unk] = {0.0}; + for (int idx = 0; idx < nunk; ++idx) { + Real corr = 0.0; + for (int r = 0; r < neq; ++r) { + corr += D[r][idx] * lambda[r]; + } + u[idx] = u0[idx] - w[idx] * corr; + } + + for (int sz = 0; sz < nz; ++sz) { + for (int sy = 0; sy < ny; ++sy) { + set_face_flux(0, 1, sy, sz, u[idx_f1(sy, sz)]); + } + } + if constexpr (DIM > 1) { + for (int sz = 0; sz < nz; ++sz) { + for (int sx = 0; sx < 2; ++sx) { + set_face_flux(1, sx, 1, sz, u[idx_f2(sx, sz)]); + } + } + } + if constexpr (DIM > 2) { + for (int sy = 0; sy < 2; ++sy) { + for (int sx = 0; sx < 2; ++sx) { + set_face_flux(2, sx, sy, 1, u[idx_f3(sx, sy)]); + } + } + } + } + } } }; diff --git a/src/utils/refinement/restriction.hpp b/src/utils/refinement/restriction.hpp index 5cbf7f73..30e174b4 100644 --- a/src/utils/refinement/restriction.hpp +++ b/src/utils/refinement/restriction.hpp @@ -54,15 +54,16 @@ struct RestrictAverage { const Coordinates_t &pco, const Coordinates_t &coarse_pco, const ParArrayND *pcoarse, const ParArrayND *pfine) { - PARTHENON_REQUIRE( - el == TE::CC || el == TE::F1 || el == TE::F2 || el == TE::F3, - "Artemis restriction only supports cell-centered and face-centered fields!"); + PARTHENON_REQUIRE(el == TE::CC || el == TE::F1 || el == TE::F2 || el == TE::F3 || + el == TE::E1 || el == TE::E2 || el == TE::E3, + "Artemis restriction only supports cell-centered, face-centered, " + "and edge-centered fields!"); constexpr bool INCLUDE_X1 = - (DIM > 0) && (el == TE::CC || el == TE::F2 || el == TE::F3); + (DIM > 0) && (el == TE::CC || el == TE::F2 || el == TE::F3 || el == TE::E1); constexpr bool INCLUDE_X2 = - (DIM > 1) && (el == TE::CC || el == TE::F3 || el == TE::F1); + (DIM > 1) && (el == TE::CC || el == TE::F3 || el == TE::F1 || el == TE::E2); constexpr bool INCLUDE_X3 = - (DIM > 2) && (el == TE::CC || el == TE::F1 || el == TE::F2); + (DIM > 2) && (el == TE::CC || el == TE::F1 || el == TE::F2 || el == TE::E3); constexpr int element_idx = static_cast(el) % 3; auto &coarse = *pcoarse; @@ -95,6 +96,12 @@ struct RestrictAverage { vol[ok][oj][oi] = coords.template GetFaceArea(); } else if constexpr (el == TE::F3) { vol[ok][oj][oi] = coords.template GetFaceArea(); + } else if constexpr (el == TE::E1) { + vol[ok][oj][oi] = coords.GetEdgeLengthX1(); + } else if constexpr (el == TE::E2) { + vol[ok][oj][oi] = coords.GetEdgeLengthX2(); + } else if constexpr (el == TE::E3) { + vol[ok][oj][oi] = coords.GetEdgeLengthX3(); } terms[ok][oj][oi] = vol[ok][oj][oi] * fine(element_idx, l, m, n, k + ok, j + oj, i + oi); @@ -107,9 +114,11 @@ struct RestrictAverage { const Real tvol = ((vol[0][0][0] + vol[0][1][0]) + (vol[0][0][1] + vol[0][1][1])) + ((vol[1][0][0] + vol[1][1][0]) + (vol[1][0][1] + vol[1][1][1])); coarse(element_idx, l, m, n, ck, cj, ci) = - (((terms[0][0][0] + terms[0][1][0]) + (terms[0][0][1] + terms[0][1][1])) + - ((terms[1][0][0] + terms[1][1][0]) + (terms[1][0][1] + terms[1][1][1]))) / - tvol; + tvol > 0.0 + ? (((terms[0][0][0] + terms[0][1][0]) + (terms[0][0][1] + terms[0][1][1])) + + ((terms[1][0][0] + terms[1][1][0]) + (terms[1][0][1] + terms[1][1][1]))) / + tvol + : 0.0; } }; diff --git a/src/utils/units.cpp b/src/utils/units.cpp index 57e89018..6c96131c 100644 --- a/src/utils/units.cpp +++ b/src/utils/units.cpp @@ -22,6 +22,7 @@ constexpr Real Year = 31536000; constexpr Real parsec = 3.0857e18; constexpr Real Rjup = 6.991100e6; constexpr Real Mjup = 1.8982e30; +constexpr Real mu0_cgs_amp = 4.0e-2 * M_PI; //---------------------------------------------------------------------------------------- //! \class Units @@ -42,6 +43,7 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { time_ = 1.; mass_ = 1.; temp_ = 1.; + current_ = 1.; } else { std::string unit_conversion = pin->GetOrAddString("artemis", "unit_conversion", "base"); @@ -50,11 +52,13 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { time_ = pin->GetOrAddReal("artemis", "time", 1.); mass_ = pin->GetOrAddReal("artemis", "mass", 1.); temp_ = pin->GetOrAddReal("artemis", "temperature", 1.); + current_ = pin->GetOrAddReal("artemis", "current", 1.); } else if (unit_conversion == "ppd") { length_ = AU; mass_ = Msolar; time_ = Year / (2. * M_PI); temp_ = 1.0; + current_ = 1.0; } else { PARTHENON_FAIL("Unit conversion not recognized! Choices are [base, ppd]"); } @@ -62,6 +66,7 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { // Remaining conversion factors energy_ = std::pow(length_, 2) * mass_ * std::pow(time_, -2); + energy_density_ = energy_ * std::pow(length_, -3); number_density_ = std::pow(length_, -3); // Store everything necessary in params for usage in analysis @@ -70,6 +75,9 @@ Units::Units(ParameterInput *pin, std::shared_ptr pkg) { pkg->AddParam("time", time_); pkg->AddParam("mass", mass_); pkg->AddParam("temp", temp_); + pkg->AddParam("current", current_); + pkg->AddParam("energy_density", energy_density_); + pkg->AddParam("magnetic_field", GetMagneticFieldCodeToPhysical()); } //---------------------------------------------------------------------------------------- @@ -90,6 +98,7 @@ Constants::Constants(Units &units) { pc_ = 1.; Year_ = 1.; Rsolar_ = 1.; + mu0_ = 1.; } else if (units.GetPhysicalUnits() == PhysicalUnits::cgs) { parthenon::constants::PhysicalConstants pc; G_ = pc.gravitational_constant; @@ -106,6 +115,7 @@ Constants::Constants(Units &units) { pc_ = parsec; Year_ = Year; Rsolar_ = Rsolar; + mu0_ = mu0_cgs_amp; } else { PARTHENON_FAIL("Unknown unit system"); } @@ -114,6 +124,7 @@ Constants::Constants(Units &units) { const Real time = units.GetTimeCodeToPhysical(); const Real mass = units.GetMassCodeToPhysical(); const Real temp = units.GetTemperatureCodeToPhysical(); + const Real current = units.GetCurrentCodeToPhysical(); const Real energy = mass * std::pow(length / time, 2); // Convert constants to code units @@ -131,6 +142,7 @@ Constants::Constants(Units &units) { pc_code_ = pc_ / length; Year_code_ = Year_ / time; Rsolar_code_ = Rsolar_ / length; + mu0_code_ = mu0_ * SQR(current) * SQR(time) / (mass * length); } } // namespace ArtemisUtils diff --git a/src/utils/units.hpp b/src/utils/units.hpp index 533df16c..0ff70dbe 100644 --- a/src/utils/units.hpp +++ b/src/utils/units.hpp @@ -58,6 +58,11 @@ class Units { return 1. / GetTemperatureCodeToPhysical(); } + KOKKOS_INLINE_FUNCTION + Real GetCurrentCodeToPhysical() const { return current_; } + KOKKOS_INLINE_FUNCTION + Real GetCurrentPhysicalToCode() const { return 1. / GetCurrentCodeToPhysical(); } + KOKKOS_INLINE_FUNCTION Real GetSpeedCodeToPhysical() const { return length_ / time_; } KOKKOS_INLINE_FUNCTION @@ -76,14 +81,26 @@ class Units { } KOKKOS_INLINE_FUNCTION - Real GetEnergyDensityCodeToPhysical() const { - return 1. / GetNumberDensityCodeToPhysical(); - } + Real GetEnergyDensityCodeToPhysical() const { return energy_density_; } KOKKOS_INLINE_FUNCTION Real GetEnergyDensityPhysicalToCode() const { return 1. / GetEnergyDensityCodeToPhysical(); } + KOKKOS_INLINE_FUNCTION + Real GetMu0CodeToPhysical() const { + return mass_ * length_ / (SQR(current_) * SQR(time_)); + } + KOKKOS_INLINE_FUNCTION + Real GetMu0PhysicalToCode() const { return 1. / GetMu0CodeToPhysical(); } + + KOKKOS_INLINE_FUNCTION + Real GetMagneticFieldCodeToPhysical() const { return mass_ / (current_ * SQR(time_)); } + KOKKOS_INLINE_FUNCTION + Real GetMagneticFieldPhysicalToCode() const { + return 1. / GetMagneticFieldCodeToPhysical(); + } + KOKKOS_INLINE_FUNCTION Real GetMassDensityCodeToPhysical() const { return mass_ * number_density_; } KOKKOS_INLINE_FUNCTION @@ -120,8 +137,10 @@ class Units { Real time_; Real mass_; Real temp_; + Real current_; Real energy_; + Real energy_density_; Real number_density_; PhysicalUnits physical_units_; @@ -190,6 +209,11 @@ class Constants { KOKKOS_INLINE_FUNCTION Real GetYearCode() const { return Year_code_; } + KOKKOS_INLINE_FUNCTION + Real GetMu0Physical() const { return mu0_; } + KOKKOS_INLINE_FUNCTION + Real GetMu0Code() const { return mu0_code_; } + private: // Physical constants in physical units Real G_; // Gravitational constant @@ -206,6 +230,7 @@ class Constants { Real pc_; // Parsec Real Year_; // Year Real Rsolar_; // Solar radius + Real mu0_; // Magnetic constant used for B^2/(2 mu0) // Physical constants in code units Real G_code_; @@ -222,6 +247,7 @@ class Constants { Real pc_code_; Real Year_code_; Real Rsolar_code_; + Real mu0_code_; }; } // namespace ArtemisUtils