use FloorParser

Author: chongchonghe

docs/markdown/parameters.md

@@ -26,6 +26,7 @@ These parameters are read in the `AMRSimulation<problem_t>::readParameters()` fu
 | derived_vars                | String        | A list of the names of derived variables that should be included in the plotfile and Ascent outputs.                                                                                                            |
 | regrid_interval             | Integer       | The number of timesteps between AMR regridding.                                                                                                                                                                 |
 | density_floor               | Float         | The minimum density value allowed in the simulation. Enforced through EnforceLimits.                                                                                                                            |
+| density_floor_expr          | String        | Optional AMReX parser expression for a spatially varying density floor. Variables: x, y, z, base_density_floor. When set, this overrides the constant floor.                                             |
 | temperature_floor           | Float         | The minimum temperature value allowed in the simulation. Enforced through EnforceLimits.                                                                                                                        |
 | max_walltime                | String        | The maximum walltime for the simulation in the format DD:HH:SS (days/hours/seconds). After 90% of this walltime elapses, the simulation will automatically stop and exit.                                       |
 | dt_cutoff                   | Float         | Timestep drop detector threshold. If the timestep drops below dt_cutoff \* current_time, the simulation aborts with an error message. This helps detect numerical instabilities early. Default: 0.0 (disabled). |

inputs/StarCluster.in

@@ -46,3 +46,5 @@ derived_vars = log_density
 
 amrex.throw_exception = 0
 amrex.signal_handling = 1
+
+density_floor_expr = "1.0 - (sqrt(x*x + y*y + z*z)/10.0)"

src/QuokkaSimulation.hpp

@@ -1693,27 +1693,20 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::postInitializati
 
 template <typename problem_t> void QuokkaSimulation<problem_t>::enforceDensityFloor(int lev, amrex::MultiFab &state_mf)
 {
-	// if (this->useDensityFloorParser_) {
-	// 	auto const density_floor_parser = this->densityFloorParserExe_.value();
-	// 	auto const density_floor_func = [=] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
-	// 									amrex::Real base_density_floor) -> amrex::Real {
-	// 		return density_floor_parser(x, y, z, base_density_floor);
-	// 	};
-	// 	HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_mf, geom[lev].data(), density_floor_func);
-	// } else {
-	// 	auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
-	// 											amrex::Real base_density_floor) -> amrex::Real {
-	// 		return densityFloor(x, y, z, base_density_floor);
-	// 	};
-	// 	HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_mf, geom[lev].data(), density_floor_func);
-	// }
-
-	// enforce density floor
-	auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
-								     amrex::Real base_density_floor) -> amrex::Real {
-		return densityFloor(x, y, z, base_density_floor);
-	};
-	HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_mf, geom[lev].data(), density_floor_func);
+	if (this->useDensityFloorParser_) {
+		auto const density_floor_parser = this->densityFloorParserExe_.value();
+		auto const density_floor_func = [=] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
+										amrex::Real base_density_floor) -> amrex::Real {
+			return density_floor_parser(x, y, z, base_density_floor);
+		};
+		HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_mf, geom[lev].data(), density_floor_func);
+	} else {
+		auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
+												amrex::Real base_density_floor) -> amrex::Real {
+			return densityFloor(x, y, z, base_density_floor);
+		};
+		HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_mf, geom[lev].data(), density_floor_func);
+	}
 }
 
 // fix-up any unphysical states created by AMR operations

src/problems/StarCluster/testStarCluster.cpp

@@ -213,7 +213,7 @@ template <> void QuokkaSimulation<StarCluster>::ComputeDerivedVar(int lev, std::
 }
 
 // A GPU helper function to set up the density floor
-AMREX_GPU_HOST_DEVICE auto localDensityFloor(amrex::Real x, amrex::Real y, amrex::Real z) -> amrex::Real
+auto localDensityFloor(amrex::Real x, amrex::Real y, amrex::Real z) -> amrex::Real
 {
 	amrex::Real r = std::sqrt(x * x + y * y + z * z);
 	const Real rho_base_floor = 0.001;
@@ -223,14 +223,6 @@ AMREX_GPU_HOST_DEVICE auto localDensityFloor(amrex::Real x, amrex::Real y, amrex
 	return custom_floor;
 }
 
-// Custom density floor for StarCluster: max(0.001, 0.1 - 0.01 * r)
-template <>
-AMREX_GPU_HOST_DEVICE auto QuokkaSimulation<StarCluster>::densityFloor(amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real /*base_density_floor*/) const
-    -> amrex::Real
-{
-	return localDensityFloor(x, y, z);
-}
-
 auto problem_main() -> int
 {
 	// read problem parameters

src/simulation.hpp

@@ -62,6 +62,7 @@ namespace filesystem = experimental::filesystem;
 #include "AMReX_Orientation.H"
 #include "AMReX_ParallelDescriptor.H"
 #include "AMReX_ParmParse.H"
+#include "AMReX_Parser.H"
 #include "AMReX_PlotFileUtil.H"
 #include "AMReX_Print.H"
 #include "AMReX_REAL.H"
@@ -213,6 +214,10 @@ template <typename problem_t> class AMRSimulation : public amrex::AmrCore
 
 	amrex::Real densityFloor_ = 0.0; // default
 	amrex::Real tempFloor_ = 0.0;	 // default
+	bool useDensityFloorParser_ = false;
+	std::string densityFloorExpr_;
+	std::optional<amrex::Parser> densityFloorParser_;
+	std::optional<amrex::ParserExecutor<4>> densityFloorParserExe_;
 
 	mutable YAML::Node simulationMetadata_;
 
@@ -857,6 +862,19 @@ template <typename problem_t> void AMRSimulation<problem_t>::readParameters()
 	// read temperature floor in K
 	pp.query("temperature_floor", tempFloor_);
 
+	// optional density floor expression (variables: x, y, z, base_density_floor)
+	densityFloorExpr_.clear();
+	pp.query("density_floor_expr", densityFloorExpr_);
+	useDensityFloorParser_ = !densityFloorExpr_.empty();
+	if (useDensityFloorParser_) {
+		densityFloorParser_.emplace(densityFloorExpr_);
+		densityFloorParser_->registerVariables({"x", "y", "z", "base_density_floor"});
+		densityFloorParserExe_ = densityFloorParser_->compile<4>();
+	} else {
+		densityFloorParser_.reset();
+		densityFloorParserExe_.reset();
+	}
+
 	// specify maximum walltime in HH:MM:SS format
 	std::string maxWalltimeInput;
 	pp.query("max_walltime", maxWalltimeInput);