diff --git a/inputs/disk/binary_rummy.par b/inputs/disk/binary_rummy.par new file mode 100644 index 00000000..68818679 --- /dev/null +++ b/inputs/disk/binary_rummy.par @@ -0,0 +1,224 @@ +# ======================================================================================== +# (C) (or copyright) 2023-2025. 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. +# ======================================================================================== + +# Disk setup without the pgen! +# This uses input file functions to specify the ICs + + +# Global vars for a simple PPD +r0 = 1.0 +rho0 = 1.0 + +dslope = -0.5 +tslope = -1.0 + +h0 = 0.05 +GM = 1.0 +q = 1e-5 + +Om0 = sqrt(GM/r0**3) +omf = 0.0 +alpha0 = 1e-3 + +# Constant Mdot IC +mdot0 = 3*pi*alpha0*h0**2 + +fn Omk(r) { + return Om0 * r**(-1.5) +} + +fn nu(r) { + return alpha0 * h0**2 * sqrt(r/r0) +} + +# These are the fields that will be called by Artemis +# Note the calling convention for fields for ICs is +# func(x1,x2,x3,time) + +# pressure supported vphi +fn vphi(r,phi,z,time) { + return Omk(r)*r * ( 1.0 - 2.5*h0**2) +} + +# Viscous vr +fn vr(r,phi,z,time) { + return -1.5*nu(r)/r +} + +fn vz(r,phi,z,time) { + return 0 +} + +fn density(r,phi,z,time) { + return mdot0/(3*pi*nu(r)) +} + +fn temperature(r,phi,z,time) { + return (h0*Om0)**2 * (r/r0)**tslope +} + + +# Begin the normal input blocks + + +coordinates = "cylindrical" +radial_spacing = "logarithmic" + + +problem_id = "disk" # problem ID: basename of output filenames + + +variables = "gas.prim.density", & + "gas.prim.velocity", & + "gas.prim.pressure" +file_type = "hdf5" # HDF5 data dump +dt = 0.05 * 2*pi + + +file_type = "rst" # Restart +dt = 2*pi + + +nlim = -1 # cycle limit +tlim = 10 * 2*pi +integrator = "rk2" # time integration algorithm +ncycle_out = 1 # interval for stdout summary info + + +nghost = 2 +# refinement = adaptive +# numlevel = 4 + +nx1 = 128 # Number of zones in X1-direction +x1min = log(.3) +x1max = log(3.0) +ix1_bc = "reflecting" # Inner-X1 boundary condition flag +ox1_bc = "reflecting" # Outer-X1 boundary condition flag + +nx2 = 256 # Number of zones in X2-direction +x2min = -pi # minimum value of X2 +x2max = pi +ix2_bc = "periodic" # Inner-X2 boundary condition flag +ox2_bc = "periodic" # Outer-X2 boundary condition 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 condition flag +ox3_bc = "periodic" # Outer-X3 boundary condition flag + + +nx1 = 32 # Number of cells in each MeshBlock, X1-dir +nx2 = 32 # Number of cells in each MeshBlock, X2-dir +nx3 = 1 # Number of cells in each MeshBlock, X3-dir + + +gas = true +# dust = true +gravity = true +nbody = true +drag = true +cooling = true +viscosity = true +# Turn on RF if the omf var is non-zero +# and turn on orbital advection if RF is on +rotating_frame = (omf > 0.0) +orbital_advection = rotating_frame + + +cfl = 0.3 +reconstruct = "plm" +riemann = "hllc-general" +dfloor = 1.0e-10 +siefloor = 1.0e-10 +refine_field = "pressure" +refine_type = "gradient" +refine_thr = 3.0 + +# Pass the fields to Artemis in this block + +velocity_x1 = "vr" +velocity_x2 = "vphi" +velocity_x3 = "vz" +density = "density" +temperature = "temperature" + + + +gamma = 1.00001 + + + +omega = omf +gm = GM + + +type = "alpha" +alpha = alpha0 + + +inner_x1 = 0.45 +outer_x1 = 2.8 +inner_x1_rate = 30.0 +outer_x1_rate = 30.0 + + +mtot = GM + + +dt_reb = 0.01 +integrator = "ias15" + +# central star +# Note that we can name the particle here as star + +mass = GM +couple = 1 + +# Note that the ./ is a shorthand for the last parent node (i.e., nbody) +<./soft> +type = "none" + +# the planet + +mass = q * GM +couple = 1 + +<./soft> +radius = 0.03 +type = "plummer" + +# initialize the particles as a binary +# Note that we can refer to the primary/secondary by name + +mass = GM +primary = "star" +secondary = "planet" +a = r0 +e = 0.0 +f = 0.0 + + + +type = "self" + + +type = "beta" +tref = "powerlaw" +beta0 = 1e-8 +# We can call the temperature function directly at r0 to set the tcyl value +tcyl = temperature(r0,0,0,0) +cyl_plaw = tslope + +# Note that there is no block because there is no pgen + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 51386b94..1082dcd7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -66,6 +66,8 @@ set (SRC_LIST pgen/disk.hpp pgen/gaussian_bump.hpp pgen/grav_slab.hpp + pgen/function_init.cpp + pgen/function_init.hpp pgen/kh.hpp pgen/linear_wave.hpp pgen/lw.hpp diff --git a/src/artemis.cpp b/src/artemis.cpp index 7139c290..bfe0f50a 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -22,6 +22,7 @@ #include "geometry/geometry.hpp" #include "gravity/gravity.hpp" #include "nbody/nbody.hpp" +#include "pgen/function_init.hpp" #include "radiation/moments/moments.hpp" #include "radiation/radiation.hpp" #include "radiation/raytrace/raytrace.hpp" @@ -39,7 +40,8 @@ namespace artemis { //---------------------------------------------------------------------------------------- //! \fn Packages_t Artemis::ProcessPackages //! \brief Process and initialize relevant packages -Packages_t ProcessPackages(std::unique_ptr &pin) { +Packages_t ProcessPackages(std::unique_ptr &pin, + Rummy::FullDeck *input_state) { Packages_t packages; // Extract artemis package and params @@ -47,7 +49,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { Params ¶ms = artemis->AllParams(); // Store selected pgen name - artemis->AddParam("pgen_name", pin->GetString("artemis", "problem")); + artemis->AddParam("pgen_name", pin->GetOrAddString("artemis", "problem", "unset")); artemis->AddParam("job_name", pin->GetString("parthenon/job", "problem_id")); artemis->AddParam("integrator", pin->GetString("parthenon/time", "integrator")); std::array nx{pin->GetInteger("parthenon/mesh", "nx1"), @@ -245,6 +247,8 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { const bool report = pin->GetOrAddBoolean("artemis", "print_artemis_config", true); if (report) ArtemisUtils::PrintArtemisConfiguration(packages); + // Store any functions from the input file to be called later + function_init::Configure(pin.get(), input_state, packages); return packages; } diff --git a/src/artemis_driver.hpp b/src/artemis_driver.hpp index 0ac79ce4..725ee775 100644 --- a/src/artemis_driver.hpp +++ b/src/artemis_driver.hpp @@ -64,7 +64,8 @@ class ArtemisDriver : public EvolutionDriver { //---------------------------------------------------------------------------------------- using TaskCollectionFnPtr = TaskCollection (*)(Mesh *pm, const Real time, const Real dt); -Packages_t ProcessPackages(std::unique_ptr &pin); +Packages_t ProcessPackages(std::unique_ptr &pin, + Rummy::FullDeck *input_state); } // namespace artemis diff --git a/src/artemis_params.yaml b/src/artemis_params.yaml index b15ca028..ca1912c3 100644 --- a/src/artemis_params.yaml +++ b/src/artemis_params.yaml @@ -18,6 +18,7 @@ artemis: _description: "(r,z,φ) axisymmetric coordinates." problem: _type: string + _default: "unset" _description: "Name of the problem pgen." amr_user: _type: bool diff --git a/src/dust/params.yaml b/src/dust/params.yaml index f1e9c416..df2bc0c1 100644 --- a/src/dust/params.yaml +++ b/src/dust/params.yaml @@ -59,6 +59,31 @@ dust: _default: 0 _description: "Gas scratch memory level" + initialize: + _type: node + _description: "Initialize dust primitives from functions defined in the input." + system: + _type: string + _description: "Coordinate system used for function arguments and velocity components. Defaults to artemis/coordinates." + cartesian: + _type: opt + _description: "Cartesian coordinates and vector components." + cylindrical: + _type: opt + _description: "Cylindrical coordinates and vector components." + spherical: + _type: opt + _description: "Spherical coordinates and vector components." + axisymmetric: + _type: opt + _description: "Axisymmetric coordinates and vector components." + density: + _type: string + _description: "Name of a scalar function with arguments (x1, x2, x3, time, size)." + velocity_x(1,2,3): + _type: string + _description: "Name of a scalar velocity-component function with arguments (x1, x2, x3, time, size)." + reconstruct: _type: string _default: plm diff --git a/src/gas/params.yaml b/src/gas/params.yaml index b2a579e9..7d00ec2f 100644 --- a/src/gas/params.yaml +++ b/src/gas/params.yaml @@ -28,6 +28,40 @@ gas: _default: 0 _description: "Gas scratch memory level" + initialize: + _type: node + _description: "Initialize gas primitives from functions defined in the input." + system: + _type: string + _description: "Coordinate system used for function arguments and velocity components. Defaults to artemis/coordinates." + cartesian: + _type: opt + _description: "Cartesian coordinates and vector components." + cylindrical: + _type: opt + _description: "Cylindrical coordinates and vector components." + spherical: + _type: opt + _description: "Spherical coordinates and vector components." + axisymmetric: + _type: opt + _description: "Axisymmetric coordinates and vector components." + density: + _type: string + _description: "Name of a scalar function with arguments (x1, x2, x3, time)." + temperature: + _type: string + _description: "Name of the temperature function. Exclusive with pressure and internal_energy." + pressure: + _type: string + _description: "Name of the pressure function. Exclusive with temperature and internal_energy." + internal_energy: + _type: string + _description: "Name of the specific-internal-energy function. Exclusive with temperature and pressure." + velocity_x(1,2,3): + _type: string + _description: "Name of a scalar velocity-component function with arguments (x1, x2, x3, time)." + de_switch: _type: Real _default: 0.0 @@ -380,4 +414,4 @@ cooling: tstop: _type: Real _description: "Turn off cooling at this time" - _default: "1e99" \ No newline at end of file + _default: "1e99" diff --git a/src/main.cpp b/src/main.cpp index b6d654cf..88ccc31d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -175,7 +175,10 @@ int main(int argc, char *argv[]) { return 1; } // Redefine parthenon defaults - pman.app_input->ProcessPackages = ProcessPackages; + auto *input_state = pman.GetRummyFullDeck(); + pman.app_input->ProcessPackages = [input_state](std::unique_ptr &pin) { + return ProcessPackages(pin, input_state); + }; // Use Parthenon default reflecting boundary conditions pman.app_input->RegisterDefaultReflectingBoundaryConditions(); diff --git a/src/pgen/function_init.cpp b/src/pgen/function_init.cpp new file mode 100644 index 00000000..103a4c3d --- /dev/null +++ b/src/pgen/function_init.cpp @@ -0,0 +1,417 @@ +//======================================================================================== +// (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. +//======================================================================================== + +// This file was created in part by generative AI + +#include "function_init.hpp" + +#include +#include +#include +#include +#include + +#include +#include + +namespace artemis { +namespace function_init { +namespace { + +struct DeviceModuleStorage { + ParArray1D code; + ParArray1D constants; + ParArray1D functions; + ParArray1D globals; +}; + +struct StorageRegistry { + std::vector modules; +}; + +template +ParArray1D Upload(const std::string &label, const T *data, const std::uint32_t count) { + if (count == 0) return {}; + ParArray1D view(label, count); + auto host = view.GetHostMirror(); + for (std::uint32_t n = 0; n < count; ++n) + host(n) = data[n]; + view.DeepCopy(host); + return view; +} + +std::string RequireFunctionName(ParameterInput *pin, const std::string &block, + const std::string &field) { + PARTHENON_REQUIRE(pin->DoesParameterExist(block, field), + "Initialization block <" + block + "> requires field '" + field + + "'."); + const auto name = pin->GetString(block, field); + PARTHENON_REQUIRE(!name.empty(), "Initialization field " + block + "/" + field + + " must name a function."); + return name; +} + +Coordinates ReadCoordinateSystem(ParameterInput *pin, const std::string &block) { + const auto default_system = pin->GetOrAddString("artemis", "coordinates", "cartesian"); + const auto system = pin->GetOrAddString(block, "system", default_system); + const auto coordinates = geometry::CoordSelect(system, ProblemDimension(pin)); + PARTHENON_REQUIRE(coordinates != Coordinates::null, + "Invalid coordinate system '" + system + "' in <" + block + ">."); + return coordinates; +} + +class FunctionPacker { + public: + FunctionPacker(Rummy::FullDeck *deck, std::shared_ptr storage, + std::vector &names) + : deck_(deck), storage_(std::move(storage)), names_(names) {} + + DeviceCallable Pack(const std::string &name, const int expected_arity) { + const auto cached = callables_.find(name); + if (cached != callables_.end()) { + PARTHENON_REQUIRE(arities_.at(name) == expected_arity, + "Callable initializer function '" + name + "' has arity " + + std::to_string(arities_.at(name)) + ", expected " + + std::to_string(expected_arity) + "."); + return cached->second; + } + + std::string error; + PARTHENON_REQUIRE(deck_->PackDeviceFunction(name, error), + "Unable to pack callable initializer function '" + name + + "' for the device: " + error); + const auto host = deck_->GetDeviceFunction(name); + PARTHENON_REQUIRE(host.entry_id < host.module.function_count, + "Packed callable initializer function '" + name + + "' has an invalid entry id."); + const auto arity = host.module.functions[host.entry_id].arity; + PARTHENON_REQUIRE(static_cast(arity) == expected_arity, + "Callable initializer function '" + name + "' has arity " + + std::to_string(arity) + ", expected " + + std::to_string(expected_arity) + "."); + + const int id = static_cast(names_.size()); + names_.push_back(name); + DeviceModuleStorage owner; + const auto prefix = "pips initializer " + std::to_string(id) + " "; + owner.code = Upload(prefix + "code", host.module.code, host.module.code_size); + owner.constants = + Upload(prefix + "constants", host.module.constants, host.module.constant_count); + owner.functions = + Upload(prefix + "functions", host.module.functions, host.module.function_count); + owner.globals = + Upload(prefix + "globals", host.module.globals, host.module.global_count); + + DeviceCallable callable; + callable.module = host.module; + callable.module.code = owner.code.data(); + callable.module.constants = owner.constants.data(); + callable.module.functions = owner.functions.data(); + callable.module.globals = owner.globals.data(); + callable.entry_id = host.entry_id; + callable.function_id = id; + storage_->modules.push_back(std::move(owner)); + callables_.emplace(name, callable); + arities_.emplace(name, static_cast(arity)); + return callable; + } + + private: + Rummy::FullDeck *deck_; + std::shared_ptr storage_; + std::vector &names_; + std::unordered_map callables_; + std::unordered_map arities_; +}; + +void ConfigureGas(ParameterInput *pin, FunctionPacker &packer, + FunctionInitConfig &config) { + constexpr const char *block = "gas/initialize"; + config.gas.enabled = true; + config.gas.input_system = ReadCoordinateSystem(pin, block); + config.gas.density = packer.Pack(RequireFunctionName(pin, block, "density"), 4); + for (int d = 0; d < 3; ++d) { + const auto field = "velocity_x" + std::to_string(d + 1); + config.gas.velocity[d] = packer.Pack(RequireFunctionName(pin, block, field), 4); + } + + const std::array, 3> options = { + std::pair{"temperature", ThermodynamicInput::temperature}, + std::pair{"pressure", ThermodynamicInput::pressure}, + std::pair{"internal_energy", ThermodynamicInput::internal_energy}}; + int count = 0; + for (const auto &[field, type] : options) { + if (pin->DoesParameterExist(block, field)) { + ++count; + config.gas.thermodynamic_input = type; + config.gas.thermodynamic = packer.Pack(RequireFunctionName(pin, block, field), 4); + } + } + PARTHENON_REQUIRE(count == 1, + " requires exactly one of temperature, pressure, or " + "internal_energy."); +} + +void ConfigureDust(ParameterInput *pin, FunctionPacker &packer, + FunctionInitConfig &config) { + constexpr const char *block = "dust/initialize"; + config.dust.enabled = true; + config.dust.input_system = ReadCoordinateSystem(pin, block); + config.dust.density = packer.Pack(RequireFunctionName(pin, block, "density"), 5); + for (int d = 0; d < 3; ++d) { + const auto field = "velocity_x" + std::to_string(d + 1); + config.dust.velocity[d] = packer.Pack(RequireFunctionName(pin, block, field), 5); + } +} + +void ConfigureMoment(ParameterInput *pin, FunctionPacker &packer, + FunctionInitConfig &config) { + constexpr const char *block = "radiation/moment/initialize"; + config.moment.enabled = true; + config.moment.input_system = ReadCoordinateSystem(pin, block); + config.moment.energy = packer.Pack(RequireFunctionName(pin, block, "energy"), 4); + for (int d = 0; d < 3; ++d) { + const auto field = "reduced_flux_x" + std::to_string(d + 1); + config.moment.reduced_flux[d] = + packer.Pack(RequireFunctionName(pin, block, field), 4); + } +} + +} // namespace + +void Configure(ParameterInput *pin, Rummy::FullDeck *deck, Packages_t &packages) { + if (Globals::is_restart) return; + + const bool has_gas = pin->DoesBlockExist("gas/initialize"); + const bool has_dust = pin->DoesBlockExist("dust/initialize"); + const bool has_moment = pin->DoesBlockExist("radiation/moment/initialize"); + const auto problem = pin->GetOrAddString("artemis", "problem", "unset"); + + const auto artemis_pkg = packages.Get("artemis"); + const bool do_gas = artemis_pkg->Param("do_gas"); + const bool do_dust = artemis_pkg->Param("do_dust"); + const bool do_moment = artemis_pkg->Param("do_moment"); + + if (problem == "unset") { + PARTHENON_REQUIRE(!do_gas || has_gas, + "artemis/problem=unset requires when gas is " + "enabled."); + PARTHENON_REQUIRE(!do_dust || has_dust, + "artemis/problem=unset requires when dust is " + "enabled."); + PARTHENON_REQUIRE( + !do_moment || has_moment, + "artemis/problem=unset requires when moment " + "radiation is enabled."); + } + + if (!(has_gas || has_dust || has_moment)) return; + PARTHENON_REQUIRE(deck != nullptr, + "Callable initialization requires a Rummy FullDeck input parser."); + + auto storage = std::make_shared(); + std::vector names; + FunctionPacker packer(deck, storage, names); + FunctionInitConfig config; + config.time = pin->GetOrAddReal("parthenon/time", "start_time", 0.0); + if (has_gas) ConfigureGas(pin, packer, config); + if (has_dust) ConfigureDust(pin, packer, config); + if (has_moment) ConfigureMoment(pin, packer, config); + + artemis_pkg->AddParam("function_init_config", config); + artemis_pkg->AddParam("function_init_names", names); + artemis_pkg->AddParam("function_init_storage", + std::static_pointer_cast(std::move(storage))); +} + +template +void Initialize(MeshBlock *pmb, ParameterInput *pin) { + auto artemis_pkg = pmb->packages.Get("artemis"); + if (!artemis_pkg->AllParams().hasKey("function_init_config")) return; + + const auto config = artemis_pkg->Param("function_init_config"); + if (!config.Enabled()) return; + + using parthenon::MakePackDescriptor; + 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()); + static auto desc_g = + MakePackDescriptor((pmb->resolved_packages).get()); + auto vg = desc_g.GetPack(md.get()); + + ArtemisUtils::EOS eos_d; + if (config.gas.enabled) { + eos_d = pmb->packages.Get("gas")->Param("eos_d"); + } + ParArray1D dust_sizes; + if (config.dust.enabled) { + dust_sizes = pmb->packages.Get("dust")->Param>("sizes"); + } + + ParArray1D error("function initializer error", 6); + auto h_error = error.GetHostMirror(); + for (int n = 0; n < 6; ++n) + h_error(n) = (n == 0) ? 0 : -1; + error.DeepCopy(h_error); + + const auto &cpars = artemis_pkg->template Param("coord_params"); + auto &pco = pmb->coords; + const auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + const auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + const auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + + pmb->par_for( + "callable primitive initialization", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + geometry::Coords coords(cpars, pco, k, j, i); + const auto xi = coords.GetCellCenter(vg, 0, k, j, i); + pips::device::DeviceVM vm; + + if (config.gas.enabled) { + std::array xo, ex1, ex2, ex3; + impl::ConvertCoordinates(coords, xi, config.gas.input_system, xo, ex1, ex2, + ex3); + pips::device::DeviceValue args[4] = { + pips::device::dv_number(xo[0]), pips::device::dv_number(xo[1]), + pips::device::dv_number(xo[2]), pips::device::dv_number(config.time)}; + Real rho, thermo; + std::array velocity; + if (!impl::Call(vm, config.gas.density, args, 4, rho, error, k, j, i) || + !impl::RequirePositive(rho, config.gas.density, error, k, j, i) || + !impl::Call(vm, config.gas.thermodynamic, args, 4, thermo, error, k, j, + i) || + !impl::RequirePositive(thermo, config.gas.thermodynamic, error, k, j, i)) { + return; + } + for (int d = 0; d < 3; ++d) { + if (!impl::Call(vm, config.gas.velocity[d], args, 4, velocity[d], error, k, j, + i)) { + return; + } + } + const auto projected = impl::ProjectVector(velocity, ex1, ex2, ex3); + Real sie = thermo; + if (config.gas.thermodynamic_input == ThermodynamicInput::temperature) { + sie = eos_d.InternalEnergyFromDensityTemperature(rho, thermo); + } else if (config.gas.thermodynamic_input == ThermodynamicInput::pressure) { + eos_d.InternalEnergyFromDensityPressure(rho, thermo, sie); + } + if (!std::isfinite(sie) || sie <= 0.0) { + impl::RecordError(error, config.gas.thermodynamic.function_id, + !std::isfinite(sie) + ? static_cast(impl::ErrorCode::non_finite) + : static_cast(impl::ErrorCode::non_positive), + k, j, i); + return; + } + for (int n = 0; n < v.GetSize(0, gas::prim::density()); ++n) { + v(0, gas::prim::density(n), k, j, i) = rho; + v(0, gas::prim::sie(n), k, j, i) = sie; + for (int d = 0; d < 3; ++d) { + v(0, gas::prim::velocity(ArtemisUtils::VI(n, d)), k, j, i) = projected[d]; + } + } + } + + if (config.dust.enabled) { + std::array xo, ex1, ex2, ex3; + impl::ConvertCoordinates(coords, xi, config.dust.input_system, xo, ex1, ex2, + ex3); + for (int n = 0; n < v.GetSize(0, dust::prim::density()); ++n) { + pips::device::DeviceValue args[5] = { + pips::device::dv_number(xo[0]), pips::device::dv_number(xo[1]), + pips::device::dv_number(xo[2]), pips::device::dv_number(config.time), + pips::device::dv_number(dust_sizes(n))}; + Real rho; + std::array velocity; + if (!impl::Call(vm, config.dust.density, args, 5, rho, error, k, j, i) || + !impl::RequirePositive(rho, config.dust.density, error, k, j, i)) { + return; + } + for (int d = 0; d < 3; ++d) { + if (!impl::Call(vm, config.dust.velocity[d], args, 5, velocity[d], error, k, + j, i)) { + return; + } + } + const auto projected = impl::ProjectVector(velocity, ex1, ex2, ex3); + v(0, dust::prim::density(n), k, j, i) = rho; + for (int d = 0; d < 3; ++d) { + v(0, dust::prim::velocity(ArtemisUtils::VI(n, d)), k, j, i) = projected[d]; + } + } + } + + if (config.moment.enabled) { + std::array xo, ex1, ex2, ex3; + impl::ConvertCoordinates(coords, xi, config.moment.input_system, xo, ex1, ex2, + ex3); + pips::device::DeviceValue args[4] = { + pips::device::dv_number(xo[0]), pips::device::dv_number(xo[1]), + pips::device::dv_number(xo[2]), pips::device::dv_number(config.time)}; + Real energy; + std::array reduced_flux; + if (!impl::Call(vm, config.moment.energy, args, 4, energy, error, k, j, i) || + !impl::RequirePositive(energy, config.moment.energy, error, k, j, i)) { + return; + } + for (int d = 0; d < 3; ++d) { + if (!impl::Call(vm, config.moment.reduced_flux[d], args, 4, reduced_flux[d], + error, k, j, i)) { + return; + } + } + const auto projected = impl::ProjectVector(reduced_flux, ex1, ex2, ex3); + for (int n = 0; n < v.GetSize(0, rad::prim::energy()); ++n) { + v(0, rad::prim::energy(n), k, j, i) = energy; + for (int d = 0; d < 3; ++d) { + v(0, rad::prim::flux(ArtemisUtils::VI(n, d)), k, j, i) = projected[d]; + } + } + } + }); + + h_error = error.GetHostMirrorAndCopy(); + if (h_error(0) != 0) { + const auto &names = + artemis_pkg->Param>("function_init_names"); + const int function_id = h_error(1); + const std::string function_name = + (function_id >= 0 && function_id < static_cast(names.size())) + ? names[function_id] + : ""; + std::ostringstream msg; + msg << "Callable initializer function '" << function_name << "' " + << impl::ErrorDescription(h_error(2)) << " at cell (k,j,i)=(" << h_error(3) << "," + << h_error(4) << "," << h_error(5) << ")."; + PARTHENON_FAIL(msg.str()); + } +} + +template void Initialize(MeshBlock *pmb, ParameterInput *pin); +template void Initialize(MeshBlock *pmb, ParameterInput *pin); +template void Initialize(MeshBlock *pmb, ParameterInput *pin); +template void Initialize(MeshBlock *pmb, ParameterInput *pin); +template void Initialize(MeshBlock *pmb, ParameterInput *pin); +template void Initialize(MeshBlock *pmb, ParameterInput *pin); + +} // namespace function_init +} // namespace artemis diff --git a/src/pgen/function_init.hpp b/src/pgen/function_init.hpp new file mode 100644 index 00000000..00b22e2d --- /dev/null +++ b/src/pgen/function_init.hpp @@ -0,0 +1,223 @@ +//======================================================================================== +// (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. +//======================================================================================== + +// This file was created in part by generative AI + +#ifndef PGEN_FUNCTION_INIT_HPP_ +#define PGEN_FUNCTION_INIT_HPP_ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/eos/eos.hpp" + +namespace Rummy { +class FullDeck; +} + +namespace artemis { +namespace function_init { + +struct DeviceCallable { + pips::device::DeviceModule module{}; + std::uint32_t entry_id = 0; + int function_id = -1; + + KOKKOS_INLINE_FUNCTION bool IsValid() const { return function_id >= 0; } +}; + +enum class ThermodynamicInput : int { + none = 0, + temperature = 1, + pressure = 2, + internal_energy = 3 +}; + +struct GasConfig { + bool enabled = false; + Coordinates input_system = Coordinates::null; + ThermodynamicInput thermodynamic_input = ThermodynamicInput::none; + DeviceCallable density; + DeviceCallable thermodynamic; + std::array velocity; +}; + +struct DustConfig { + bool enabled = false; + Coordinates input_system = Coordinates::null; + DeviceCallable density; + std::array velocity; +}; + +struct MomentConfig { + bool enabled = false; + Coordinates input_system = Coordinates::null; + DeviceCallable energy; + std::array reduced_flux; +}; + +struct FunctionInitConfig { + Real time = 0.0; + GasConfig gas; + DustConfig dust; + MomentConfig moment; + + bool Enabled() const { return gas.enabled || dust.enabled || moment.enabled; } +}; + +void Configure(ParameterInput *pin, Rummy::FullDeck *deck, Packages_t &packages); + +namespace impl { + +enum class ErrorCode : int { + none = 0, + vm_error_base = 100, + non_numeric = 200, + non_finite = 201, + non_positive = 202 +}; + +template +KOKKOS_INLINE_FUNCTION void RecordError(const ERROR &error, const int function_id, + const int error_code, const int k, const int j, + const int i) { + if (Kokkos::atomic_compare_exchange(&error(0), 0, 1) == 0) { + error(1) = function_id; + error(2) = error_code; + error(3) = k; + error(4) = j; + error(5) = i; + } +} + +template +KOKKOS_INLINE_FUNCTION bool +Call(pips::device::DeviceVM &vm, const DeviceCallable &callable, + const pips::device::DeviceValue *args, const std::uint32_t argc, Real &value, + const ERROR &error, const int k, const int j, const int i) { + pips::device::DeviceValue result{}; + const auto status = vm.run(callable.module, callable.entry_id, args, argc, &result); + if (status != pips::device::DeviceStatus::OK) { + RecordError(error, callable.function_id, + static_cast(ErrorCode::vm_error_base) + static_cast(status), k, + j, i); + return false; + } + if (!pips::device::dv_is_number(result)) { + RecordError(error, callable.function_id, static_cast(ErrorCode::non_numeric), k, + j, i); + return false; + } + value = static_cast(pips::device::dv_as_number(result)); + if (!std::isfinite(value)) { + RecordError(error, callable.function_id, static_cast(ErrorCode::non_finite), k, + j, i); + return false; + } + return true; +} + +template +KOKKOS_INLINE_FUNCTION bool +RequirePositive(const Real value, const DeviceCallable &callable, const ERROR &error, + const int k, const int j, const int i) { + if (value > 0.0) return true; + RecordError(error, callable.function_id, static_cast(ErrorCode::non_positive), k, + j, i); + return false; +} + +template +KOKKOS_INLINE_FUNCTION void +ConvertCoordinates(const COORDS &coords, const XV &xi, const Coordinates input_system, + std::array &xo, std::array &ex1, + std::array &ex2, std::array &ex3) { + if (input_system == Coordinates::cartesian) { + const auto converted = coords.ConvertToCartWithVec(xi); + xo = std::get<0>(converted); + ex1 = std::get<1>(converted); + ex2 = std::get<2>(converted); + ex3 = std::get<3>(converted); + } else if (input_system == Coordinates::axisymmetric) { + const auto converted = coords.ConvertToAxiWithVec(xi); + xo = std::get<0>(converted); + ex1 = std::get<1>(converted); + ex2 = std::get<2>(converted); + ex3 = std::get<3>(converted); + } else if (geometry::is_spherical(input_system)) { + const auto converted = coords.ConvertToSphWithVec(xi); + xo = std::get<0>(converted); + ex1 = std::get<1>(converted); + ex2 = std::get<2>(converted); + ex3 = std::get<3>(converted); + } else { + const auto converted = coords.ConvertToCylWithVec(xi); + xo = std::get<0>(converted); + ex1 = std::get<1>(converted); + ex2 = std::get<2>(converted); + ex3 = std::get<3>(converted); + } +} + +KOKKOS_INLINE_FUNCTION std::array ProjectVector(const std::array &vin, + const std::array &ex1, + const std::array &ex2, + const std::array &ex3) { + return {vin[0] * ex1[0] + vin[1] * ex1[1] + vin[2] * ex1[2], + vin[0] * ex2[0] + vin[1] * ex2[1] + vin[2] * ex2[2], + vin[0] * ex3[0] + vin[1] * ex3[1] + vin[2] * ex3[2]}; +} + +inline std::string ErrorDescription(const int code) { + if (code >= static_cast(ErrorCode::vm_error_base) && + code < static_cast(ErrorCode::non_numeric)) { + return "PIPS device VM status " + + std::to_string(code - static_cast(ErrorCode::vm_error_base)); + } + if (code == static_cast(ErrorCode::non_numeric)) { + return "returned a non-numeric value"; + } + if (code == static_cast(ErrorCode::non_finite)) { + return "returned a non-finite value"; + } + if (code == static_cast(ErrorCode::non_positive)) { + return "returned a non-positive density, energy, pressure, or temperature"; + } + return "failed with unknown error code " + std::to_string(code); +} + +} // namespace impl + +template +void Initialize(MeshBlock *pmb, ParameterInput *pin); + +} // namespace function_init +} // namespace artemis + +#endif // PGEN_FUNCTION_INIT_HPP_ diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 9ae84c8f..9cae2a61 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -26,6 +26,7 @@ #include "constant.hpp" #include "crooked_pipe.hpp" #include "disk.hpp" +#include "function_init.hpp" #include "gaussian_bump.hpp" #include "geometry/geometry.hpp" #include "grav_slab.hpp" @@ -47,7 +48,8 @@ namespace artemis { template void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_INSTRUMENT - std::string name = pin->GetString("artemis", "problem"); + function_init::Initialize(pmb, pin); + std::string name = pin->GetOrAddString("artemis", "problem", "unset"); if (name == "advection") { advection::ProblemGenerator(pmb, pin); } else if (name == "beam") { @@ -84,6 +86,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { thermalization::ProblemGenerator(pmb, pin); } else if (name == "crooked_pipe") { crooked_pipe::ProblemGenerator(pmb, pin); + } else if (name == "unset") { + return; } else { PARTHENON_FAIL("Invalid problem name!"); } @@ -98,7 +102,7 @@ void InitMeshBlockData(MeshBlock *pmb, ParameterInput *pin) { geometry::InitBlockGeom(pmb, pin); - std::string name = pin->GetString("artemis", "problem"); + std::string name = pin->GetOrAddString("artemis", "problem", "unset"); if (name == "beam") { beam::InitBeamParams(pmb, pin); } else if (name == "conduction") { diff --git a/src/radiation/params.yaml b/src/radiation/params.yaml index e2a98ddd..cbbd7b9f 100644 --- a/src/radiation/params.yaml +++ b/src/radiation/params.yaml @@ -87,6 +87,30 @@ radiation: _type: int _description: "Scratch level for radiation" _default: 0 + initialize: + _type: node + _description: "Initialize moment primitives from functions defined in the input." + system: + _type: string + _description: "Coordinate system used for function arguments and reduced-flux components. Defaults to artemis/coordinates." + cartesian: + _type: opt + _description: "Cartesian coordinates and vector components." + cylindrical: + _type: opt + _description: "Cylindrical coordinates and vector components." + spherical: + _type: opt + _description: "Spherical coordinates and vector components." + axisymmetric: + _type: opt + _description: "Axisymmetric coordinates and vector components." + energy: + _type: string + _description: "Name of a scalar radiation-energy function with arguments (x1, x2, x3, time)." + reduced_flux_x(1,2,3): + _type: string + _description: "Name of a scalar reduced-flux function with arguments (x1, x2, x3, time)." refine_field: _type: string _default: "none"