Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

revert to first-order space+time flux when RK2 flux fails #221

Merged
merged 22 commits into from
Sep 12, 2023
Merged

revert to first-order space+time flux when RK2 flux fails #221

Changes from 11 commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
b0ec43d
use Roe-average wavespeeds
BenWibking Dec 22, 2022
0a372c3
fix bug in roe-averaged sound speed
BenWibking Dec 23, 2022
14ee2a3
use LLF solver with first-order flux
BenWibking Dec 23, 2022
1f4e56e
revert to first-order space+time flux
BenWibking Dec 27, 2022
a1ec044
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Dec 27, 2022
f2a6669
cherry-pick fix for SphericalCollapse
BenWibking Dec 28, 2022
1621ff1
remove Riemann solver template parameter
BenWibking Dec 28, 2022
9d9ecc7
fix typo
BenWibking Dec 28, 2022
f85f901
Merge branch 'development' into fofc-llf
BenWibking Jan 2, 2023
701c2ec
remove LLF solver (unused)
BenWibking Jan 2, 2023
d8901de
Merge branch 'development' into fofc-llf
BenWibking Jan 4, 2023
824ca16
Merge branch 'development' into fofc-llf
markkrumholz Jan 31, 2023
fc114b1
Merge branch 'development' into fofc-llf
BenWibking Feb 1, 2023
aea4a0d
Merge branch 'development' into fofc-llf
BenWibking Feb 8, 2023
15af6ba
Merge branch 'development' into fofc-llf
BenWibking Feb 13, 2023
088e6e1
Merge branch 'development' into fofc-llf
BenWibking Aug 14, 2023
dd16eb7
Update src/RadhydroSimulation.hpp
BenWibking Aug 14, 2023
fc79048
Update src/RadhydroSimulation.hpp
BenWibking Aug 14, 2023
59ef782
Merge branch 'development' into fofc-llf
BenWibking Sep 11, 2023
65ff209
Merge github.com:quokka-astro/quokka into fofc-llf
BenWibking Sep 12, 2023
9ea9dad
update HydroHighMach tolerance
BenWibking Sep 12, 2023
74900c6
Merge branch 'fofc-llf' of github.com:quokka-astro/quokka into fofc-llf
BenWibking Sep 12, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 57 additions & 20 deletions src/RadhydroSimulation.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,15 @@
#include "AMReX_MultiFab.H"
#include "AMReX_MultiFabUtil.H"
#include "AMReX_ParmParse.H"
#include "AMReX_Periodicity.H"
#include "AMReX_PhysBCFunct.H"
#include "AMReX_Print.H"
#include "AMReX_REAL.H"
#include "AMReX_Utility.H"
#include "AMReX_YAFluxRegister.H"

#include "CloudyCooling.hpp"
#include "HLLC.hpp"
#include "SimulationData.hpp"
#include "hydro_system.hpp"
#include "hyperbolic_system.hpp"
Expand Down Expand Up @@ -766,7 +768,14 @@ void RadhydroSimulation<problem_t>::advanceHydroAtLevelWithRetries(int lev, amre
break;
}
}
AMREX_ALWAYS_ASSERT(success); // crash if we exceeded max_retries

if (!success) {
// crash, we have exceeded max_retries
amrex::Print() << "\nQUOKKA FATAL ERROR\n"
<< "Hydro update exceeded max_retries on level " << lev << ". Cannot continue, crashing...\n"
<< std::endl;
amrex::Abort();
}
}

template <typename problem_t> auto RadhydroSimulation<problem_t>::isCflViolated(int lev, amrex::Real time, amrex::Real dt_actual) -> bool
Expand Down Expand Up @@ -814,21 +823,36 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
amrex::MultiFab state_inter_cc_(grids[lev], dmap[lev], Physics_Indices<problem_t>::nvarTotal_cc, nghost_cc_);
state_inter_cc_.setVal(0); // prevent assert in fillBoundaryConditions when radiation is enabled

// Stage 1 of RK2-SSP
{
// update ghost zones [old timestep]
fillBoundaryConditions(state_old_cc_tmp, state_old_cc_tmp, lev, time, quokka::centering::cc, quokka::direction::na, PreInterpState,
PostInterpState);
// create temporary multifabs for combined RK2 flux
std::array<amrex::MultiFab, AMREX_SPACEDIM> flux_rk2;
auto ba = grids[lev];
auto dm = dmap[lev];
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
auto ba_face = amrex::convert(ba, amrex::IntVect::TheDimensionVector(idim));
flux_rk2[idim] = amrex::MultiFab(ba_face, dm, ncompHydro_, 0);
flux_rk2[idim].setVal(0);
}

// update ghost zones [old timestep]
fillBoundaryConditions(state_old_cc_tmp, state_old_cc_tmp, lev, time, quokka::centering::cc, quokka::direction::na, PreInterpState, PostInterpState);

// check state validity
AMREX_ASSERT(!state_old_cc_tmp.contains_nan(0, state_old_cc_tmp.nComp()));
AMREX_ASSERT(!state_old_cc_tmp.contains_nan()); // check ghost cells

// check state validity
AMREX_ASSERT(!state_old_cc_tmp.contains_nan(0, state_old_cc_tmp.nComp()));
AMREX_ASSERT(!state_old_cc_tmp.contains_nan()); // check ghost cells
auto [FOfluxArrays, FOfaceVel] = computeFOHydroFluxes(state_old_cc_tmp, ncompHydro_, lev);

// Stage 1 of RK2-SSP
{
// advance all grids on local processor (Stage 1 of integrator)
auto const &stateOld = state_old_cc_tmp;
auto &stateNew = state_inter_cc_;
auto [fluxArrays, faceVel] = computeHydroFluxes(stateOld, ncompHydro_, lev);

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
amrex::MultiFab::Saxpy(flux_rk2[idim], 0.5, fluxArrays[idim], 0, 0, ncompHydro_, 0);
}

amrex::MultiFab rhs(grids[lev], dmap[lev], ncompHydro_, 0);
amrex::iMultiFab redoFlag(grids[lev], dmap[lev], 1, 0);
redoFlag.setVal(quokka::redoFlag::none);
Expand All @@ -843,10 +867,11 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
if (ncells_bad > 0) {
if (Verbose()) {
amrex::Print() << "[FOFC-1] flux correcting " << ncells_bad << " cells on level " << lev << "\n";
const amrex::IntVect cell_idx = redoFlag.maxIndex(0);
amrex::print_state(stateNew, cell_idx);
}

// replace fluxes around troubled cells with Godunov fluxes
auto [FOfluxArrays, FOfaceVel] = computeFOHydroFluxes(stateOld, ncompHydro_, lev);
replaceFluxes(fluxArrays, FOfluxArrays, redoFlag);
replaceFluxes(faceVel, FOfaceVel, redoFlag); // needed for dual energy

Expand All @@ -856,10 +881,14 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
HydroSystem<problem_t>::PredictStep(stateOld, stateNew, rhs, dt_lev, ncompHydro_, redoFlag);

amrex::Gpu::streamSynchronizeAll(); // just in case
amrex::Long ncells_bad = redoFlag.sum(0);
amrex::Long ncells_bad = static_cast<int>(redoFlag.sum(0));
BenWibking marked this conversation as resolved.
Show resolved Hide resolved
if (ncells_bad > 0) {
// FOFC failed
if (Verbose()) {
const amrex::IntVect cell_idx = redoFlag.maxIndex(0);
// print cell state
amrex::Print() << "[FOFC-1] Flux correction failed:\n";
amrex::print_state(stateNew, cell_idx);
amrex::Print() << "[FOFC-1] failed for " << ncells_bad << " cells on level " << lev << "\n";
}
if (abortOnFofcFailure_ != 0) {
Expand Down Expand Up @@ -898,37 +927,46 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
auto &stateFinal = state_new_cc_[lev];
auto [fluxArrays, faceVel] = computeHydroFluxes(stateInter, ncompHydro_, lev);

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
amrex::MultiFab::Saxpy(flux_rk2[idim], 0.5, fluxArrays[idim], 0, 0, ncompHydro_, 0);
}

amrex::MultiFab rhs(grids[lev], dmap[lev], ncompHydro_, 0);
amrex::iMultiFab redoFlag(grids[lev], dmap[lev], 1, 0);
redoFlag.setVal(quokka::redoFlag::none);

HydroSystem<problem_t>::ComputeRhsFromFluxes(rhs, fluxArrays, dx, ncompHydro_);
HydroSystem<problem_t>::ComputeRhsFromFluxes(rhs, flux_rk2, dx, ncompHydro_);
HydroSystem<problem_t>::AddInternalEnergyPdV(rhs, stateInter, dx, faceVel, redoFlag);
HydroSystem<problem_t>::AddFluxesRK2(stateFinal, stateOld, stateInter, rhs, dt_lev, ncompHydro_, redoFlag);
HydroSystem<problem_t>::PredictStep(stateOld, stateFinal, rhs, dt_lev, ncompHydro_, redoFlag);

// do first-order flux correction (FOFC)
amrex::Gpu::streamSynchronizeAll(); // just in case
int ncells_bad = redoFlag.sum(0);
int ncells_bad = static_cast<int>(redoFlag.sum(0));
BenWibking marked this conversation as resolved.
Show resolved Hide resolved
if (ncells_bad > 0) {
if (Verbose()) {
amrex::Print() << "[FOFC-2] flux correcting " << ncells_bad << " cells on level " << lev << "\n";
const amrex::IntVect cell_idx = redoFlag.maxIndex(0);
amrex::print_state(stateFinal, cell_idx);
}

// replace fluxes around troubled cells with Godunov fluxes
auto [FOfluxArrays, FOfaceVel] = computeFOHydroFluxes(stateInter, ncompHydro_, lev);
replaceFluxes(fluxArrays, FOfluxArrays, redoFlag);
replaceFluxes(flux_rk2, FOfluxArrays, redoFlag);
replaceFluxes(faceVel, FOfaceVel, redoFlag); // needed for dual energy

// re-do RK update
HydroSystem<problem_t>::ComputeRhsFromFluxes(rhs, fluxArrays, dx, ncompHydro_);
HydroSystem<problem_t>::ComputeRhsFromFluxes(rhs, flux_rk2, dx, ncompHydro_);
HydroSystem<problem_t>::AddInternalEnergyPdV(rhs, stateInter, dx, faceVel, redoFlag);
HydroSystem<problem_t>::AddFluxesRK2(stateFinal, stateOld, stateInter, rhs, dt_lev, ncompHydro_, redoFlag);
HydroSystem<problem_t>::PredictStep(stateOld, stateFinal, rhs, dt_lev, ncompHydro_, redoFlag);

amrex::Gpu::streamSynchronizeAll(); // just in case
amrex::Long ncells_bad = redoFlag.sum(0);
if (ncells_bad > 0) {
// FOFC failed
if (Verbose()) {
const amrex::IntVect cell_idx = redoFlag.maxIndex(0);
// print cell state
amrex::Print() << "[FOFC-2] Flux correction failed:\n";
amrex::print_state(stateFinal, cell_idx);
amrex::Print() << "[FOFC-2] failed for " << ncells_bad << " cells on level " << lev << "\n";
}
if (abortOnFofcFailure_ != 0) {
Expand Down Expand Up @@ -979,7 +1017,7 @@ void RadhydroSimulation<problem_t>::replaceFluxes(std::array<amrex::MultiFab, AM
// left of zone i, and -1.0*flux(i+1) is the flux *into* zone i through
// the interface on the right of zone i.

amrex::IntVect ng{AMREX_D_DECL(0, 0, 0)};
amrex::IntVect ng{AMREX_D_DECL(0, 0, 0)}; // TODO(ben): must include 1 ghost zone

amrex::ParallelFor(redoFlag, ng, ncomp, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k, int n) noexcept {
if (redoFlag_arrs[bx](i, j, k) == quokka::redoFlag::redo) {
Expand Down Expand Up @@ -1116,7 +1154,6 @@ auto RadhydroSimulation<problem_t>::computeFOHydroFluxes(amrex::MultiFab const &

auto ba = grids[lev];
auto dm = dmap[lev];
const int flatteningGhost = 2;
const int reconstructRange = 1;

// allocate temporary MultiFabs
Expand Down