variable density floor feature: simplified test problem

.ci/azure-pipelines-amdgpu.yml

@@ -110,12 +110,12 @@ jobs:
         CCACHE_BASEDIR: $(Build.SourcesDirectory)
         CCACHE_NOHASHDIR: 1
       inputs:
-        cmakeArgs: '--build . --parallel 16'
+        cmakeArgs: '--build . --parallel 16 --target GravRadParticle3D'
 
     - task: CMake@1
       displayName: 'Run CTest'
       inputs:
-        cmakeArgs: '-E chdir . ctest -E "RadMatterCoupling*" -T Test --output-on-failure'
+        cmakeArgs: '-E chdir . ctest -R GravRadParticle3D -T Test --output-on-failure'
 
     - task: PublishTestResults@2
       inputs:

.ci/azure-pipelines.yml

@@ -109,12 +109,12 @@ jobs:
         CCACHE_BASEDIR: $(Build.SourcesDirectory)
         CCACHE_NOHASHDIR: 1
       inputs:
-        cmakeArgs: '--build . --parallel 16'
+        cmakeArgs: '--build . --parallel 16 --target GravRadParticle3D'
 
     - task: CMake@1
       displayName: 'Run CTest'
       inputs:
-        cmakeArgs: '-E chdir . ctest -E "RadMatterCoupling*" -T Test --output-on-failure'
+        cmakeArgs: '-E chdir . ctest -R GravRadParticle3D -T Test --output-on-failure'
 
     - task: PublishTestResults@2
       inputs:

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/GravRadParticle3D.in

@@ -53,4 +53,7 @@ quokka.proj.particles = CIC_particles # List of particle types to include in out
 
 max_timesteps = 2
 
-tiny_profiler.enabled = 0
+tiny_profiler.enabled = 0
+
+density_floor = 1.0e-7
+density_floor_expr = "1.0e-7 * (3.0 - sqrt(x*x + y*y + z*z) / 2.0)"

src/QuokkaSimulation.hpp

@@ -245,6 +245,8 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
 		eos_init(small_temp, small_dens);
 	}
 
+	[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto densityFloor(amrex::Real x, amrex::Real y, amrex::Real z,
+										 amrex::Real base_density_floor) const -> amrex::Real;
 	[[nodiscard]] static auto getScalarVariableNames() -> std::vector<std::string>;
 	void defineComponentNames();
 	void defineDefaultPlotfileVariables();
@@ -297,6 +299,8 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
 	// fix-up states
 	void FixupState(int level) override;
 
+	void enforceDensityFloor(int lev, amrex::MultiFab &state_mf) override;
+
 	// implement FillPatch function
 	void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf, int icomp, int ncomp, quokka::centering cen, quokka::direction dir,
 		       FillPatchType fptype) override;
@@ -1024,6 +1028,15 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::print_multifab_f
 	    mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept { printf("%f\n", mf_fc[bx](i, j, k, Physics_Indices<problem_t>::mhdFirstIndex)); });
 }
 
+template <typename problem_t>
+AMREX_GPU_HOST_DEVICE auto QuokkaSimulation<problem_t>::densityFloor(amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real base_density_floor) const
+    -> amrex::Real
+{
+	// Default implementation returns base density floor
+	// Problem generators can specialize this function for custom behavior
+	return base_density_floor;
+}
+
 template <typename problem_t> auto QuokkaSimulation<problem_t>::computeComponentErrors() -> std::vector<std::tuple<std::string, amrex::Real, amrex::Real>>
 {
 	// returns a vector of tuples: (component name, absolute error, relative error)
@@ -1678,14 +1691,31 @@ 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);
+	}
+}
+
 // fix-up any unphysical states created by AMR operations
 // (e.g., caused by the flux register or from interpolation)
 template <typename problem_t> void QuokkaSimulation<problem_t>::FixupState(int lev)
 {
 	const BL_PROFILE("QuokkaSimulation::FixupState()");
 
-	// fix hydro state
-	HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_new_cc_[lev]);
+	enforceDensityFloor(lev, state_new_cc_[lev]);
 
 	// sync internal energy and total energy
 	HydroSystem<problem_t>::SyncDualEnergy(state_new_cc_[lev], state_new_fc_[lev]);
@@ -2188,7 +2218,7 @@ auto QuokkaSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_old
 		}
 
 		// prevent vacuum
-		HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateNew_cc);
+		enforceDensityFloor(lev, stateNew_cc);
 
 		if (useDualEnergy_ == 1) {
 			// sync internal energy (requires positive density)
@@ -2299,7 +2329,7 @@ auto QuokkaSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_old
 		}
 
 		// prevent vacuum
-		HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateFinal_cc);
+		enforceDensityFloor(lev, stateFinal_cc);
 
 		if (useDualEnergy_ == 1) {
 			// sync internal energy (requires positive density)

src/hydro/hydro_system.hpp

@@ -17,6 +17,7 @@
 #include "AMReX.H"
 #include "AMReX_Array4.H"
 #include "AMReX_BLassert.H"
+#include "AMReX_Geometry.H"
 #include "AMReX_MultiFabUtil.H"
 #include "AMReX_Print.H"
 #include "AMReX_REAL.H"
@@ -145,7 +146,9 @@ template <typename problem_t> class HydroSystem : public HyperbolicSystem<proble
 
 	AMREX_GPU_DEVICE static auto GetGradFixedPotential(amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> posvec) -> amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>;
 
-	static void EnforceLimits(amrex::Real densityFloor, amrex::Real tempFloor, amrex::MultiFab &state_mf);
+	template <typename DensityFloorFunc>
+	static void EnforceLimits(amrex::Real densityFloor, amrex::Real tempFloor, amrex::MultiFab &state_mf, amrex::GeometryData const &geom,
+				  DensityFloorFunc const &density_floor_func);
 
 	static void AddInternalEnergyPdV(amrex::MultiFab &rhs_mf, amrex::MultiFab const &consVar_mf,
 					 std::array<amrex::MultiFab, AMREX_SPACEDIM> const &cons_fc_mf, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx,
@@ -933,19 +936,38 @@ void HydroSystem<problem_t>::FlattenShocks(amrex::MultiFab const &q_mf, amrex::M
 
 // to ensure that physical quantities are within reasonable
 // floors and ceilings which can be set in the param file
-template <typename problem_t> void HydroSystem<problem_t>::EnforceLimits(amrex::Real const densityFloor, amrex::Real const tempFloor, amrex::MultiFab &state_mf)
+template <typename problem_t>
+template <typename DensityFloorFunc>
+void HydroSystem<problem_t>::EnforceLimits(amrex::Real const densityFloor, amrex::Real const tempFloor, amrex::MultiFab &state_mf,
+					   amrex::GeometryData const &geom, DensityFloorFunc const &density_floor_func)
 {
 	auto state = state_mf.arrays();
+	auto const *const prob_lo = geom.ProbLo();
+	auto const *const dx = geom.CellSize();
 
 	amrex::ParallelFor(state_mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
+		amrex::Real const x = prob_lo[0] + (static_cast<amrex::Real>(i) + static_cast<amrex::Real>(0.5)) * dx[0];
+#if (AMREX_SPACEDIM >= 2)
+		amrex::Real const y = prob_lo[1] + (static_cast<amrex::Real>(j) + static_cast<amrex::Real>(0.5)) * dx[1];
+#else
+		amrex::Real const y = 0.0;
+#endif
+#if (AMREX_SPACEDIM == 3)
+		amrex::Real const z = prob_lo[2] + (static_cast<amrex::Real>(k) + static_cast<amrex::Real>(0.5)) * dx[2];
+#else
+		amrex::Real const z = 0.0;
+#endif
+
+		amrex::Real localDensityFloor = density_floor_func(x, y, z, densityFloor);
+
 		// Enforce density floor (do not adjust energies here!!)
 		amrex::Real rho_new = NAN;
 		{
 			amrex::Real const rho = state[bx](i, j, k, density_index);
 			rho_new = rho;
 
-			if (rho < densityFloor) {
-				rho_new = densityFloor;
+			if (rho < localDensityFloor) {
+				rho_new = localDensityFloor;
 				state[bx](i, j, k, density_index) = rho_new;
 
 				if (nscalars_ > 0) {
@@ -984,8 +1006,8 @@ template <typename problem_t> void HydroSystem<problem_t>::EnforceLimits(amrex::
 		if constexpr (Physics_Traits<problem_t>::is_dust_enabled) {
 			for (int g = 0; g < Physics_Traits<problem_t>::nDustGroups; ++g) {
 				amrex::Real dust_rho = state[bx](i, j, k, dustDensity_index + g * numDustVars_);
-				if (dust_rho < densityFloor) {
-					state[bx](i, j, k, dustDensity_index + g * numDustVars_) = densityFloor;
+				if (dust_rho < localDensityFloor) {
+					state[bx](i, j, k, dustDensity_index + g * numDustVars_) = localDensityFloor;
 				}
 			}
 		}

src/linear_advection/AdvectionSimulation.hpp

@@ -105,6 +105,8 @@ template <typename problem_t> class AdvectionSimulation : public AMRSimulation<p
 
 	void FixupState(int lev) override;
 
+	void enforceDensityFloor(int lev, amrex::MultiFab &state_mf) override;
+
 	// tag cells for refinement
 	void refineGrid(int lev, amrex::TagBoxArray &tags, amrex::Real time, int ngrow) override;
 
@@ -259,6 +261,11 @@ template <typename problem_t> void AdvectionSimulation<problem_t>::FixupState(in
 	// fix negative states
 }
 
+template <typename problem_t> void AdvectionSimulation<problem_t>::enforceDensityFloor(int /*lev*/, amrex::MultiFab & /*state_mf*/)
+{
+	// enforce density floor
+}
+
 template <typename problem_t>
 void AdvectionSimulation<problem_t>::computeReferenceSolution(amrex::MultiFab &ref, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
 							      amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_lo,

src/problems/GravRadParticle3D/testGravRadParticle3D.cpp

@@ -17,6 +17,10 @@
 #include "particles/PhysicsParticles.hpp"
 #include "radiation/radiation_system.hpp"
 #include "util/BC.hpp"
+#include "util/fextract.hpp"
+#ifdef HAVE_PYTHON
+#include "util/matplotlibcpp.h"
+#endif
 
 struct ParticleProblem {
 };
@@ -29,7 +33,7 @@ constexpr double initial_Egas = 1.0e-5;
 constexpr double c = 100.0;	   // speed of light
 constexpr double chat = 2.0;	   // reduced speed of light
 constexpr double kappa0 = 1.0e-20; // opacity
-constexpr double rho0 = 1.0e-8;
+constexpr double rho0 = 1.0e-9;
 constexpr double m_H = C::m_p + C::m_e;
 
 const double lum1 = 1.0;
@@ -142,6 +146,12 @@ auto QuokkaSimulation<ParticleProblem>::ComputeProjections(const amrex::Directio
 	return proj;
 }
 
+auto localDensityFloor(amrex::Real x, amrex::Real y, amrex::Real z) -> amrex::Real
+{
+	// density_floor_expr = "1.0e-7 * (3.0 - sqrt(x*x + y*y + z*z) / 2.0)"
+	return std::max(1.0e-7, 1.0e-7 * (3.0 - std::sqrt(x * x + y * y + z * z) / 2.0));
+}
+
 auto problem_main() -> int
 {
 	// Problem parameters
@@ -156,11 +166,60 @@ auto problem_main() -> int
 	// initialize
 	sim.setInitialConditions();
 
-	// evolve
-	sim.evolve();
+	// read output variables
+	auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 2, 0.0, true); // z direction
+	const int nz = static_cast<int>(position.size());
 
 	int status = 0;
 
+	// extract density and check floor
+	std::vector<double> zs(nz);
+	std::vector<double> rho_z(nz);
+	std::vector<double> custom_floor_z(nz);
+	amrex::Real min_density = std::numeric_limits<amrex::Real>::max();
+	amrex::Real min_density_ratio = std::numeric_limits<amrex::Real>::max();
+	for (int i = 0; i < nz; ++i) {
+		amrex::Real const z = position[i];
+		custom_floor_z.at(i) = localDensityFloor(0.0, 0.0, z); // note that the real x and y are 0.5 * delta_x
+		amrex::Real const rho = values.at(RadSystem<ParticleProblem>::gasDensity_index)[i];
+		zs.at(i) = z;
+		rho_z.at(i) = rho;
+		min_density = std::min(min_density, rho);
+		min_density_ratio = std::min(min_density_ratio, rho / custom_floor_z.at(i));
+	}
+
+	// Check that custom floor is working: min_density_ratio should not be smaller than 0.98
+	if (min_density_ratio > 0.99) {
+		amrex::Print() << "Custom density floor test PASSED: min density ratio = " << min_density_ratio << " > 0.99\n";
+	} else {
+		amrex::Print() << "Custom density floor test FAILED: min density ratio = " << min_density_ratio << " < 0.99\n";
+		status = 1;
+	}
+
+#ifdef HAVE_PYTHON
+	// Plot results
+	matplotlibcpp::clf();
+	matplotlibcpp::ylim(0.9e-7, 3.1e-7);
+
+	std::map<std::string, std::string> rho_args;
+	std::map<std::string, std::string> custom_floor_args;
+	rho_args["label"] = "rho";
+	rho_args["linestyle"] = "-";
+	rho_args["color"] = "C0";
+	custom_floor_args["label"] = "custom floor";
+	custom_floor_args["linestyle"] = "--";
+	custom_floor_args["color"] = "C1";
+	matplotlibcpp::plot(zs, rho_z, rho_args);
+	matplotlibcpp::plot(zs, custom_floor_z, custom_floor_args);
+
+	matplotlibcpp::legend();
+	matplotlibcpp::title("Custom density floor: 1.0e-7*(3-r/2)");
+	matplotlibcpp::save("./grav_rad_particle_3d_density_floor_z.pdf");
+#endif // HAVE_PYTHON
+
+	// evolve
+	sim.evolve();
+
 	// compute total radiation energy
 	const double total_Erad_over_vol = sim.state_new_cc_[0].sum(RadSystem<ParticleProblem>::radEnergy_index);
 	const double dx = sim.Geom(0).CellSize(0);
@@ -266,11 +325,15 @@ auto problem_main() -> int
 
 		const double rel_err_tol = 1.0e-7;
 		const double rel_position_error_tol = t_sim < 1.0 ? 2.0e-4 : 2.0e-3;
-		status = 1;
 		if (rel_err < rel_err_tol && rel_position_error_cicrad < rel_position_error_tol && rel_position_error_cic < rel_position_error_tol &&
 		    rel_position_error_rad < rel_position_error_tol) {
-			status = 0;
 			amrex::Print() << "Relative error within tolerance.\n";
+		} else {
+			status = 1;
+			amrex::Print() << "Relative error beyond tolerance: rel_err = " << rel_err
+				       << ", rel_position_error_cicrad = " << rel_position_error_cicrad
+				       << ", rel_position_error_cic = " << rel_position_error_cic << ", rel_position_error_rad = " << rel_position_error_rad
+				       << "\n";
 		}
 
 		amrex::Print() << "Exact positions of the CICRad particles should be: " << exact_x << ", " << exact_y << ", " << exact_z << "\n";
@@ -294,8 +357,7 @@ auto problem_main() -> int
 		amrex::Print() << "Relative L1 norm on Rad particle positions = " << rel_position_error_rad << "\n";
 
 		// Cleanup and exit
-		amrex::Print() << "Finished."
-			       << "\n";
+		amrex::Print() << "Finished.\n";
 	}
 
 	return status;