Skip to content

Commit d548a80

Browse files
BenWibkingpre-commit-ci[bot]
BenWibking
and
pre-commit-ci[bot]
authored
Add spatially-variable density floor (#1520)
### Description Allow specifying a spatially-variable density floor as a runtime ParmParse expression. Usage example: ``` density_floor = 1.0 density_floor_expr = "1.0e-2*base_density_floor*exp(-x/0.25)" ``` When evaluating `density_floor_expr` in each cell, `base_density_floor` is replaced with the value of `density_floor` in the input file. The spatial coordinates are available as `x`, `y`, and `z`. Multiple subexpressions can be separated using semicolons, e.g., the radius can be used with `r=sqrt(x*x+y*y+z*z); base_density_floor / pow(r,2)` ### Related issues N/A ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [ ] I have added a description (see above). - [ ] I have added a link to any related issues (if applicable; see above). - [ ] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [ ] *(For quokka-astro org members)* I have manually triggered the GPU tests with the magic comment `/azp run`. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
1 parent 10549da commit d548a80

11 files changed

Lines changed: 510 additions & 24 deletions

File tree

‎docs/markdown/parameters.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@ These parameters are read in the `AMRSimulation<problem_t>::readParameters()` fu
2626
| derived_vars | String | A list of the names of derived variables that should be included in the plotfile and Ascent outputs. |
2727
| regrid_interval | Integer | The number of timesteps between AMR regridding. |
2828
| density_floor | Float | The minimum density value allowed in the simulation. Enforced through EnforceLimits. |
29+
| 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. |
30+
| debug_density_floor_plot | Boolean (0/1) | If set to 1, adds a derived field `density_floor_dbg` to plotfiles to visualize the spatially varying density floor. Default: 0 (disabled). |
2931
| temperature_floor | Float | The minimum temperature value allowed in the simulation. Enforced through EnforceLimits. |
3032
| 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. |
3133
| 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). |

‎extern/amrex‎

‎inputs/HydrostaticAtmosphere.in‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
# Hydrostatic exponential atmosphere density floor unit test inputs
2+
3+
# *****************************************************************
4+
# Problem size and geometry
5+
# *****************************************************************
6+
geometry.prob_lo = 0.0 0.0 0.0
7+
geometry.prob_hi = 1.0 1.0 1.0
8+
quokka.bc = ext_dir periodic periodic
9+
10+
# *****************************************************************
11+
# VERBOSITY
12+
# *****************************************************************
13+
amr.v = 0 # verbosity in Amr
14+
15+
# *****************************************************************
16+
# Resolution and refinement
17+
# *****************************************************************
18+
amr.n_cell = 8 8 8
19+
amr.max_level = 0 # number of levels = max_level + 1
20+
amr.blocking_factor = 1 # grid size must be divisible by this
21+
22+
do_reflux = 0
23+
do_subcycle = 0
24+
25+
density_floor = 1.0
26+
atmosphere_scale_height = 0.25
27+
density_floor_expr = "1.0e-2*base_density_floor*exp(-x/0.25)"

‎src/.clang-tidy‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
---
2-
Checks: 'clang-diagnostic-*,clang-analyzer-*,*,-fuchsia-*,-llvmlibc-*,-altera-*,-misc-non-private-member-variables-in-classes,-cppcoreguidelines-non-private-member-variables-in-classes,-google-runtime-references,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers,-bugprone-easily-swappable-parameters,-readability-identifier-length,-cppcoreguidelines-pro-bounds-constant-array-index,-readability-function-cognitive-complexity,-llvm-header-guard,-performance-enum-size,-misc-include-cleaner,-cppcoreguidelines-misleading-capture-default-by-value,-readability-math-missing-parentheses,-misc-use-internal-linkage,-boost-*,-cppcoreguidelines-pro-bounds-pointer-arithmetic,-abseil-*,-llvm-*,-modernize-use-ranges'
2+
Checks: 'clang-diagnostic-*,clang-analyzer-*,*,-fuchsia-*,-llvmlibc-*,-altera-*,-misc-non-private-member-variables-in-classes,-cppcoreguidelines-non-private-member-variables-in-classes,-google-runtime-references,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers,-bugprone-easily-swappable-parameters,-readability-identifier-length,-cppcoreguidelines-pro-bounds-constant-array-index,-readability-function-cognitive-complexity,-llvm-header-guard,-performance-enum-size,-misc-include-cleaner,-cppcoreguidelines-misleading-capture-default-by-value,-readability-math-missing-parentheses,-misc-use-internal-linkage,-boost-*,-cppcoreguidelines-pro-bounds-pointer-arithmetic,-abseil-*,-llvm-*,-modernize-use-ranges,-portability-template-virtual-member-function'
33
WarningsAsErrors: ''
44
HeaderFilterRegex: ''
55
FormatStyle: none

‎src/QuokkaSimulation.hpp‎

Lines changed: 130 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ namespace filesystem = experimental::filesystem;
3030
#include <limits>
3131
#include <map>
3232
#include <memory>
33+
#include <source_location>
3334
#include <string>
3435
#include <tuple>
3536
#include <unordered_map>
@@ -245,13 +246,17 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
245246
eos_init(small_temp, small_dens);
246247
}
247248

249+
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto densityFloor(amrex::Real x, amrex::Real y, amrex::Real z,
250+
amrex::Real base_density_floor) const -> amrex::Real;
248251
[[nodiscard]] static auto getScalarVariableNames() -> std::vector<std::string>;
249252
void defineComponentNames();
250253
void defineDefaultPlotfileVariables();
251254
void readParmParse();
252255
void rereadRuntimeParameters(); // Re-read parameters to ensure runtime values override compile-time settings
253256

254257
void checkHydroStates(amrex::MultiFab &mf, std::array<amrex::MultiFab, AMREX_SPACEDIM> &mf_fc, char const *file, int line);
258+
void CheckHydroStates(amrex::MultiFab &mf, std::array<amrex::MultiFab, AMREX_SPACEDIM> &mf_fc,
259+
std::source_location const &location = std::source_location::current());
255260
void computeMaxSignalLocal(int level) override;
256261
void printCellProperties(int lev, amrex::IntVect const &index) override;
257262
void preCalculateInitialConditions() override;
@@ -287,6 +292,7 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
287292

288293
// compute derived variables
289294
void ComputeDerivedVar(int lev, std::string const &dname, amrex::MultiFab &mf, int ncomp) const override;
295+
void ComputeDensityFloorDebug(int lev, amrex::MultiFab &mf, int ncomp) const override;
290296

291297
// compute projected vars
292298
[[nodiscard]] auto ComputeProjections(amrex::Direction dir) const -> std::unordered_map<std::string, amrex::BaseFab<amrex::Real>> override;
@@ -783,11 +789,18 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::printCellPropert
783789
}
784790
}
785791

786-
#if !defined(NDEBUG)
787-
#define CHECK_HYDRO_STATES(mf, mf_fc) checkHydroStates(mf, mf_fc, __FILE__, __LINE__)
792+
template <typename problem_t>
793+
void QuokkaSimulation<problem_t>::CheckHydroStates(amrex::MultiFab &mf, std::array<amrex::MultiFab, AMREX_SPACEDIM> &mf_fc,
794+
std::source_location const &location)
795+
{
796+
#ifndef NDEBUG
797+
checkHydroStates(mf, mf_fc, location.file_name(), static_cast<int>(location.line()));
788798
#else
789-
#define CHECK_HYDRO_STATES(mf, mf_fc)
799+
static_cast<void>(mf);
800+
static_cast<void>(mf_fc);
801+
static_cast<void>(location);
790802
#endif
803+
}
791804

792805
template <typename problem_t>
793806
void QuokkaSimulation<problem_t>::checkHydroStates(amrex::MultiFab &mf, std::array<amrex::MultiFab, AMREX_SPACEDIM> &mf_fc, char const *file, int line)
@@ -967,6 +980,66 @@ auto QuokkaSimulation<problem_t>::addStrangSplitSourcesWithBuiltin(amrex::MultiF
967980
template <typename problem_t> void QuokkaSimulation<problem_t>::ComputeDerivedVar(int lev, std::string const &dname, amrex::MultiFab &mf, const int ncomp) const
968981
{
969982
// compute derived variables and save in 'mf' -- user should implement
983+
(void)lev;
984+
(void)dname;
985+
(void)mf;
986+
(void)ncomp;
987+
}
988+
989+
template <typename problem_t> void QuokkaSimulation<problem_t>::ComputeDensityFloorDebug(int lev, amrex::MultiFab &mf, int ncomp) const
990+
{
991+
auto const ncomp_out = ncomp;
992+
auto const geom_data = geom[lev].data();
993+
auto const prob_lo = geom_data.ProbLo();
994+
auto const dx = geom_data.CellSize();
995+
auto const density_floor = densityFloor_;
996+
auto const ngrow = mf.nGrow();
997+
998+
if (this->useDensityFloorParser_) {
999+
auto const density_floor_parser = this->densityFloorParserExe_.value();
1000+
for (amrex::MFIter iter(mf); iter.isValid(); ++iter) {
1001+
amrex::Box const &box = iter.growntilebox(ngrow);
1002+
auto const &arr = mf.array(iter);
1003+
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
1004+
amrex::Real const x = prob_lo[0] + (static_cast<amrex::Real>(i) + static_cast<amrex::Real>(0.5)) * dx[0];
1005+
#if (AMREX_SPACEDIM >= 2)
1006+
amrex::Real const y = prob_lo[1] + (static_cast<amrex::Real>(j) + static_cast<amrex::Real>(0.5)) * dx[1];
1007+
#else
1008+
amrex::Real const y = 0.0;
1009+
#endif
1010+
#if (AMREX_SPACEDIM == 3)
1011+
amrex::Real const z = prob_lo[2] + (static_cast<amrex::Real>(k) + static_cast<amrex::Real>(0.5)) * dx[2];
1012+
#else
1013+
amrex::Real const z = 0.0;
1014+
#endif
1015+
arr(i, j, k, ncomp_out) = density_floor_parser(x, y, z, density_floor);
1016+
});
1017+
}
1018+
} else {
1019+
auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
1020+
amrex::Real base_density_floor) -> amrex::Real {
1021+
return densityFloor(x, y, z, base_density_floor);
1022+
};
1023+
for (amrex::MFIter iter(mf); iter.isValid(); ++iter) {
1024+
amrex::Box const &box = iter.growntilebox(ngrow);
1025+
auto const &arr = mf.array(iter);
1026+
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
1027+
amrex::Real const x = prob_lo[0] + (static_cast<amrex::Real>(i) + static_cast<amrex::Real>(0.5)) * dx[0];
1028+
#if (AMREX_SPACEDIM >= 2)
1029+
amrex::Real const y = prob_lo[1] + (static_cast<amrex::Real>(j) + static_cast<amrex::Real>(0.5)) * dx[1];
1030+
#else
1031+
amrex::Real const y = 0.0;
1032+
#endif
1033+
#if (AMREX_SPACEDIM == 3)
1034+
amrex::Real const z = prob_lo[2] + (static_cast<amrex::Real>(k) + static_cast<amrex::Real>(0.5)) * dx[2];
1035+
#else
1036+
amrex::Real const z = 0.0;
1037+
#endif
1038+
arr(i, j, k, ncomp_out) = density_floor_func(x, y, z, density_floor);
1039+
});
1040+
}
1041+
}
1042+
amrex::Gpu::streamSynchronize();
9701043
}
9711044

9721045
template <typename problem_t>
@@ -1024,6 +1097,14 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::print_multifab_f
10241097
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)); });
10251098
}
10261099

1100+
template <typename problem_t>
1101+
AMREX_GPU_HOST_DEVICE auto QuokkaSimulation<problem_t>::densityFloor(amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real base_density_floor) const
1102+
-> amrex::Real
1103+
{
1104+
amrex::ignore_unused(x, y, z);
1105+
return base_density_floor;
1106+
}
1107+
10271108
template <typename problem_t> auto QuokkaSimulation<problem_t>::computeComponentErrors() -> std::vector<std::tuple<std::string, amrex::Real, amrex::Real>>
10281109
{
10291110
// returns a vector of tuples: (component name, absolute error, relative error)
@@ -1288,7 +1369,7 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::advanceSingleTim
12881369
}
12891370

12901371
// check hydro states before update (this can be caused by the flux register!)
1291-
CHECK_HYDRO_STATES(state_old_cc_[lev], state_old_fc_[lev]);
1372+
CheckHydroStates(state_old_cc_[lev], state_old_fc_[lev]);
12921373

12931374
// advance hydro
12941375
if constexpr (Physics_Traits<problem_t>::is_hydro_enabled) {
@@ -1300,21 +1381,21 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::advanceSingleTim
13001381
}
13011382

13021383
// check hydro states after hydro update
1303-
CHECK_HYDRO_STATES(state_new_cc_[lev], state_new_fc_[lev]);
1384+
CheckHydroStates(state_new_cc_[lev], state_new_fc_[lev]);
13041385

13051386
// subcycle radiation
13061387
if constexpr (Physics_Traits<problem_t>::is_radiation_enabled) {
13071388
subcycleRadiationAtLevel(lev, time, dt_lev, fr_as_crse, fr_as_fine);
13081389
}
13091390

13101391
// check hydro states after radiation update
1311-
CHECK_HYDRO_STATES(state_new_cc_[lev], state_new_fc_[lev]);
1392+
CheckHydroStates(state_new_cc_[lev], state_new_fc_[lev]);
13121393

13131394
// compute any operator-split terms here (user-defined)
13141395
computeAfterLevelAdvance(lev, time, dt_lev, ncycle);
13151396

13161397
// check hydro states after user work
1317-
CHECK_HYDRO_STATES(state_new_cc_[lev], state_new_fc_[lev]);
1398+
CheckHydroStates(state_new_cc_[lev], state_new_fc_[lev]);
13181399

13191400
// check state validity
13201401
AMREX_ASSERT(!state_new_cc_[lev].contains_nan(0, state_new_cc_[lev].nComp()));
@@ -1685,7 +1766,20 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::FixupState(int l
16851766
const BL_PROFILE("QuokkaSimulation::FixupState()");
16861767

16871768
// fix hydro state
1688-
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_new_cc_[lev]);
1769+
if (this->useDensityFloorParser_) {
1770+
auto const density_floor_parser = this->densityFloorParserExe_.value();
1771+
auto const density_floor_func = [=] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
1772+
amrex::Real base_density_floor) -> amrex::Real {
1773+
return density_floor_parser(x, y, z, base_density_floor);
1774+
};
1775+
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_new_cc_[lev], geom[lev], density_floor_func);
1776+
} else {
1777+
auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
1778+
amrex::Real base_density_floor) -> amrex::Real {
1779+
return densityFloor(x, y, z, base_density_floor);
1780+
};
1781+
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, state_new_cc_[lev], geom[lev], density_floor_func);
1782+
}
16891783

16901784
// sync internal energy and total energy
16911785
HydroSystem<problem_t>::SyncDualEnergy(state_new_cc_[lev], state_new_fc_[lev]);
@@ -2188,7 +2282,20 @@ auto QuokkaSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_old
21882282
}
21892283

21902284
// prevent vacuum
2191-
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateNew_cc);
2285+
if (this->useDensityFloorParser_) {
2286+
auto const density_floor_parser = this->densityFloorParserExe_.value();
2287+
auto const density_floor_func = [=] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
2288+
amrex::Real base_density_floor) -> amrex::Real {
2289+
return density_floor_parser(x, y, z, base_density_floor);
2290+
};
2291+
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateNew_cc, geom[lev], density_floor_func);
2292+
} else {
2293+
auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
2294+
amrex::Real base_density_floor) -> amrex::Real {
2295+
return densityFloor(x, y, z, base_density_floor);
2296+
};
2297+
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateNew_cc, geom[lev], density_floor_func);
2298+
}
21922299

21932300
if (useDualEnergy_ == 1) {
21942301
// sync internal energy (requires positive density)
@@ -2299,7 +2406,20 @@ auto QuokkaSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_old
22992406
}
23002407

23012408
// prevent vacuum
2302-
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateFinal_cc);
2409+
if (this->useDensityFloorParser_) {
2410+
auto const density_floor_parser = this->densityFloorParserExe_.value();
2411+
auto const density_floor_func = [=] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
2412+
amrex::Real base_density_floor) -> amrex::Real {
2413+
return density_floor_parser(x, y, z, base_density_floor);
2414+
};
2415+
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateFinal_cc, geom[lev], density_floor_func);
2416+
} else {
2417+
auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
2418+
amrex::Real base_density_floor) -> amrex::Real {
2419+
return densityFloor(x, y, z, base_density_floor);
2420+
};
2421+
HydroSystem<problem_t>::EnforceLimits(densityFloor_, tempFloor_, stateFinal_cc, geom[lev], density_floor_func);
2422+
}
23032423

23042424
if (useDualEnergy_ == 1) {
23052425
// sync internal energy (requires positive density)

‎src/hydro/hydro_system.hpp‎

Lines changed: 28 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#include "AMReX.H"
1818
#include "AMReX_Array4.H"
1919
#include "AMReX_BLassert.H"
20+
#include "AMReX_Geometry.H"
2021
#include "AMReX_MultiFabUtil.H"
2122
#include "AMReX_Print.H"
2223
#include "AMReX_REAL.H"
@@ -145,7 +146,9 @@ template <typename problem_t> class HydroSystem : public HyperbolicSystem<proble
145146

146147
AMREX_GPU_DEVICE static auto GetGradFixedPotential(amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> posvec) -> amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>;
147148

148-
static void EnforceLimits(amrex::Real densityFloor, amrex::Real tempFloor, amrex::MultiFab &state_mf);
149+
template <typename DensityFloorFunc>
150+
static void EnforceLimits(amrex::Real densityFloor, amrex::Real tempFloor, amrex::MultiFab &state_mf, amrex::Geometry const &geom,
151+
DensityFloorFunc const &density_floor_func);
149152

150153
static void AddInternalEnergyPdV(amrex::MultiFab &rhs_mf, amrex::MultiFab const &consVar_mf,
151154
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
933936

934937
// to ensure that physical quantities are within reasonable
935938
// floors and ceilings which can be set in the param file
936-
template <typename problem_t> void HydroSystem<problem_t>::EnforceLimits(amrex::Real const densityFloor, amrex::Real const tempFloor, amrex::MultiFab &state_mf)
939+
template <typename problem_t>
940+
template <typename DensityFloorFunc>
941+
void HydroSystem<problem_t>::EnforceLimits(amrex::Real const densityFloor, amrex::Real const tempFloor, amrex::MultiFab &state_mf, amrex::Geometry const &geom,
942+
DensityFloorFunc const &density_floor_func)
937943
{
938944
auto state = state_mf.arrays();
945+
auto const prob_lo = geom.ProbLoArray();
946+
auto const dx = geom.CellSizeArray();
939947

940948
amrex::ParallelFor(state_mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
949+
amrex::Real const x = prob_lo[0] + (static_cast<amrex::Real>(i) + static_cast<amrex::Real>(0.5)) * dx[0];
950+
#if (AMREX_SPACEDIM >= 2)
951+
amrex::Real const y = prob_lo[1] + (static_cast<amrex::Real>(j) + static_cast<amrex::Real>(0.5)) * dx[1];
952+
#else
953+
amrex::Real const y = 0.0;
954+
#endif
955+
#if (AMREX_SPACEDIM == 3)
956+
amrex::Real const z = prob_lo[2] + (static_cast<amrex::Real>(k) + static_cast<amrex::Real>(0.5)) * dx[2];
957+
#else
958+
amrex::Real const z = 0.0;
959+
#endif
960+
961+
amrex::Real localDensityFloor = density_floor_func(x, y, z, densityFloor);
962+
941963
// Enforce density floor (do not adjust energies here!!)
942964
amrex::Real rho_new = NAN;
943965
{
944966
amrex::Real const rho = state[bx](i, j, k, density_index);
945967
rho_new = rho;
946968

947-
if (rho < densityFloor) {
948-
rho_new = densityFloor;
969+
if (rho < localDensityFloor) {
970+
rho_new = localDensityFloor;
949971
state[bx](i, j, k, density_index) = rho_new;
950972

951973
if (nscalars_ > 0) {
@@ -984,8 +1006,8 @@ template <typename problem_t> void HydroSystem<problem_t>::EnforceLimits(amrex::
9841006
if constexpr (Physics_Traits<problem_t>::is_dust_enabled) {
9851007
for (int g = 0; g < Physics_Traits<problem_t>::nDustGroups; ++g) {
9861008
amrex::Real dust_rho = state[bx](i, j, k, dustDensity_index + g * numDustVars_);
987-
if (dust_rho < densityFloor) {
988-
state[bx](i, j, k, dustDensity_index + g * numDustVars_) = densityFloor;
1009+
if (dust_rho < localDensityFloor) {
1010+
state[bx](i, j, k, dustDensity_index + g * numDustVars_) = localDensityFloor;
9891011
}
9901012
}
9911013
}

‎src/io/DiagPlotfile.H‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -129,9 +129,11 @@ template <typename problem_t> void DiagPlotfile::processDiag(int a_nstep, const
129129
// Write physics particles
130130
if (m_particleTypes.empty()) {
131131
// Write all particle types
132+
sim->particleRegister_.redistribute(0, 0);
132133
sim->particleRegister_.writePlotFile(plotfilename);
133134
} else {
134135
// Write only specified particle types
136+
sim->particleRegister_.redistribute(0, 0);
135137
sim->particleRegister_.writePlotFileFiltered(plotfilename, m_particleTypes);
136138
}
137139
#endif

0 commit comments

Comments
 (0)