Skip to content

Commit

Permalink
Add redo if hydro check does not pass (#20)
Browse files Browse the repository at this point in the history 
- A validity check is placed after each step
- redo is triggered when negative density or pressure are found
  • Loading branch information
@chengcli
chengcli committed Mar 9, 2024
1 parent e9e9587 commit 78d31ee
Show file tree
Hide file tree
Showing 5 changed files with 123 additions and 98 deletions.
14 changes: 14 additions & 0 deletions main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,18 +116,32 @@ int main(int argc, char **argv) {
}

std::uint64_t mbcnt = 0;
bool redo = false;

while ((pmesh->time < pmesh->tlim) &&
(pmesh->nlim < 0 || pmesh->ncycle < pmesh->nlim)) {
if (Globals::my_rank == 0)
pmesh->OutputCycleDiagnostics();

if (!redo) {
pmesh->SaveAllStates();
}

for (int stage=1; stage<=ptlist->nstages; ++stage) {
ptlist->DoTaskListOneStage(pmesh, stage);
}

pmesh->UserWorkInLoop();

if (pmesh->CheckAllValid()) {
redo = false;
} else {
pmesh->LoadAllStates();
pmesh->DecreaseTimeStep();
redo = true;
continue;
}

pmesh->ncycle++;
pmesh->time += pmesh->dt;
mbcnt += pmesh->nbtotal;
Expand Down
2 changes: 2 additions & 0 deletions src/impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,8 @@ class MeshBlock::Impl {

private:
MeshBlock const *pmy_block_;
std::vector<char> mbdata_;
bool state_valid_flag_;
};

#endif // SRC_IMPL_HPP_
46 changes: 46 additions & 0 deletions src/mesh_redo.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
// athena
#include <athena/mesh/mesh.hpp>

// canoe
#include <configure.hpp>
#include <impl.hpp>

#ifdef MPI_PARALLEL
#include <mpi.h>
#endif // MPI_PARALLEL

void Mesh::SaveAllStates()
{
for (int b = 0; b < nblocal; ++b) {
my_blocks(b)->pimpl->SaveAllStates();
}
}

void Mesh::LoadAllStates()
{
for (int b = 0; b < nblocal; ++b) {
my_blocks(b)->pimpl->LoadAllStates();
}
}

void Mesh::DecreaseTimeStep()
{
dt /= 2;
}

bool Mesh::CheckAllValid() const
{
bool valid = true;
for (int b = 0; b < nblocal; ++b)
if (!my_blocks(b)->pimpl->IsStateValid())
valid = false;

bool all_valid = valid;

#ifdef MPI_PARALLEL
// MPI call to check if all blocks are valid
MPI_Allreduce(&valid, &all_valid, 1, MPI_C_BOOL, MPI_LAND, MPI_COMM_WORLD);
#endif // MPI_PARALLEL

return all_valid;
}
87 changes: 37 additions & 50 deletions src/snap/eos/ideal_moist_hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#include <cmath> // sqrt()
#include <sstream>
#include <stdexcept>
#include <iomanip>

// application
#include <application/exceptions.hpp>
Expand Down Expand Up @@ -36,11 +37,17 @@
std::string str_grid_ij(AthenaArray<Real> const& var, int n, int k, int j,
int i, int il, int iu, int jl, int ju) {
std::stringstream msg;
for (int ii = std::max(i - 3, il); ii <= std::min(i + 3, iu); ++ii) {
msg << "i = " << ii << " ";
for (int jj = std::max(j - 3, jl); jj <= std::min(j + 3, ju); ++jj)
msg << var(n, k, jj, ii) << " ";
msg << std::endl;
msg << std::endl << std::setw(8) << " ";
for (int jj = std::max(j - NGHOST, jl); jj <= std::min(j + NGHOST, ju); ++jj) {
msg << "j = " << jj << ", ";
}
msg << std::endl;

for (int ii = std::max(i - NGHOST, il); ii <= std::min(i + NGHOST, iu); ++ii) {
msg << "i = " << ii << " |";
for (int jj = std::max(j - NGHOST, jl); jj <= std::min(j + NGHOST, ju); ++jj)
msg << std::setw(8) << var(n, k, jj, ii) << ", ";
msg << " |" << std::endl;
}

return msg.str();
Expand Down Expand Up @@ -71,12 +78,13 @@ void EquationOfState::ConservedToPrimitive(
AthenaArray<Real>& cons, const AthenaArray<Real>& prim_old,
const FaceField& b, AthenaArray<Real>& prim, AthenaArray<Real>& bcc,
Coordinates* pco, int il, int iu, int jl, int ju, int kl, int ku) {
Real gm1 = GetGamma() - 1.0;
std::stringstream msg;
auto pthermo = Thermodynamics::GetInstance();
auto pmb = pmy_block_;
pmb->pimpl->SetStateValid();

apply_vapor_limiter(&cons, pmy_block_);

Real gm1 = GetGamma() - 1.0;
for (int k = kl; k <= ku; ++k)
for (int j = jl; j <= ju; ++j)
for (int i = il; i <= iu; ++i) {
Expand All @@ -86,18 +94,21 @@ void EquationOfState::ConservedToPrimitive(
Real& u_m3 = cons(IM3, k, j, i);
Real& u_e = cons(IEN, k, j, i);

if (u_e < 0.) {
throw RuntimeError("ideal_moist_hydro", "negative energy");
}

Real& w_d = prim(IDN, k, j, i);
Real& w_vx = prim(IVX, k, j, i);
Real& w_vy = prim(IVY, k, j, i);
Real& w_vz = prim(IVZ, k, j, i);
Real& w_p = prim(IPR, k, j, i);

#ifdef ENABLED_GLOG
Real density = 0.;
for (int n = 0; n <= NVAPOR; ++n) density += cons(n, k, j, i);
w_d = density;
Real di = 1. / density;

if ((w_d < 0.) || std::isnan(w_d)) pmb->pimpl->SetStateInvalid();
#ifdef ENABLE_GLOG
LOG_IF(ERROR, std::isnan(w_d) || (w_d < density_floor_))
<< "cycle = " << pmb->pmy_mesh->ncycle << ", "
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", w_d = " << w_d << std::endl;
Expand All @@ -106,14 +117,6 @@ void EquationOfState::ConservedToPrimitive(
<< str_grid_ij(cons, IDN, k, j, i, il, iu, jl, ju);
#endif

// apply density floor, without changing momentum or energy
u_d = (u_d > density_floor_) ? u_d : density_floor_;

Real density = 0.;
for (int n = 0; n <= NVAPOR; ++n) density += cons(n, k, j, i);
w_d = density;
Real di = 1. / density;

// mass mixing ratio
for (int n = 1; n <= NVAPOR; ++n)
prim(n, k, j, i) = cons(n, k, j, i) * di;
Expand Down Expand Up @@ -141,43 +144,27 @@ void EquationOfState::ConservedToPrimitive(
static_cast<GnomonicEquiangle*>(pco)->GetCosineCell(k, j);
#endif

int decay_factor = 1, iter = 0, max_iter = 10;
do {
iter++;
u_m1 /= decay_factor;
u_m2 /= decay_factor;
u_m3 /= decay_factor;

// Real di = 1.0/u_d;
w_vx = u_m1 * di;
w_vy = u_m2 * di;
w_vz = u_m3 * di;
// Real di = 1.0/u_d;
w_vx = u_m1 * di;
w_vy = u_m2 * di;
w_vz = u_m3 * di;

#ifdef CUBED_SPHERE
cs::CovariantToContravariant(prim.at(k, j, i), cos_theta);
cs::CovariantToContravariant(prim.at(k, j, i), cos_theta);
#endif

// internal energy
KE = 0.5 * (u_m1 * w_vx + u_m2 * w_vy + u_m3 * w_vz);
w_p = gm1 * (u_e - KE) * feps / fsig;
// internal energy
KE = 0.5 * (u_m1 * w_vx + u_m2 * w_vy + u_m3 * w_vz);
w_p = gm1 * (u_e - KE) * feps / fsig;

if ((w_p < 0.) || std::isnan(w_p)) pmb->pimpl->SetStateInvalid();
#ifdef ENABLE_GLOG
LOG_IF(ERROR, w_p < 0.)
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", w_p = " << w_p << std::endl;
LOG_IF(WARNING, (w_p < 0.) || std::isnan(w_p))
<< "cycle = " << pmb->pmy_mesh->ncycle << ", "
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", w_p = " << w_p << std::endl;
#endif
decay_factor = 2;
if (iter > max_iter) {
throw RuntimeError("ideal_moist_hydro", "negative pressure");
}
} while (w_p < 0.);

// apply pressure floor, correct total energy
u_e = (w_p > pressure_floor_)
? u_e
: ((pressure_floor_ / gm1) * fsig / feps + KE);
w_p = (w_p > pressure_floor_) ? w_p : pressure_floor_;
}
}

Expand Down
72 changes: 24 additions & 48 deletions src/snap/implicit/solve_implicit_3d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,57 +119,33 @@ void ImplicitSolver::SolveImplicit3D(AthenaArray<Real> &du, AthenaArray<Real> &w
}
#endif // CUBED_SPHERE

// save a copy of du for redo
for (int n = 0; n < NHYDRO; ++n)
for (int i = 0; i < du_.GetDim1(); ++i)
du(n, k, j, i) = du_(n, k, j, i);

bool redo;
int factor = 1, iter = 0, max_iter = 5;
do {
redo = false;
iter++;
// reset du
if (factor > 1) {
for (int n = 0; n < NHYDRO; ++n)
for (int i = 0; i < du_.GetDim1(); ++i)
du_(n, k, j, i) = du(n, k, j, i);
}

// do implicit
if (implicit_flag_ & (1 << 3))
FullCorrection(du_, w, dt * factor, k, j, is, ie);
else
PartialCorrection(du_, w, dt * factor, k, j, is, ie);

// check for negative density and internal energy
for (int i = is; i <= ie; i++) {
// do implicit
if (implicit_flag_ & (1 << 3))
FullCorrection(du_, w, dt, k, j, is, ie);
else
PartialCorrection(du_, w, dt, k, j, is, ie);

#ifdef ENABLE_GLOG
LOG_IF(ERROR, ph->u(IEN,k,j,i) + du_(IEN,k,j,i) < 0.)
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", u[IEN] = " << ph->u(IEN,k,j,i) + du_(IEN,k,j,i)
<< ", u[IVX] = " << ph->u(IVX,k,j,i) << std::endl;

LOG_IF(ERROR, ph->u(IDN,k,j,i) + du_(IDN,k,j,i) < 0.)
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", u[IDN] = " << ph->u(IDN,k,j,i) + du_(IDN,k,j,i)
<< ", u[IVX] = " << ph->u(IVX,k,j,i) << std::endl;
// check for negative density and internal energy
for (int i = is; i <= ie; i++) {
LOG_IF(WARNING, ph->u(IEN,k,j,i) + du_(IEN,k,j,i) < 0.)
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", u[IEN](before) = " << ph->u(IEN,k,j,i) + du(IEN,k,j,i)
<< ", u[IEN](after) = " << ph->u(IEN,k,j,i) + du_(IEN,k,j,i)
<< ", u[IVX](before) = " << ph->u(IVX,k,j,i) + du(IVX,k,j,i) << std::endl
<< ", u[IVX](after) = " << ph->u(IVX,k,j,i) + du_(IVX,k,j,i) << std::endl;

LOG_IF(WARNING, ph->u(IDN,k,j,i) + du_(IDN,k,j,i) < 0.)
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", u[IDN](before) = " << ph->u(IDN,k,j,i) + du(IDN,k,j,i)
<< ", u[IDN](after) = " << ph->u(IDN,k,j,i) + du_(IDN,k,j,i)
<< ", u[IVX](before) = " << ph->u(IVX,k,j,i) + du(IVX,k,j,i) << std::endl
<< ", u[IVX](after) = " << ph->u(IVX,k,j,i) + du_(IVX,k,j,i) << std::endl;
}
#endif // ENABLE_GLOG
if ((ph->u(IEN,k,j,i) + du_(IEN,k,j,i) < 0.) ||
(ph->u(IDN,k,j,i) + du_(IDN,k,j,i) < 0.)) {
factor *= 2;
redo = true;
break;
}
}

if (iter > max_iter) {
throw RuntimeError("solver_implicit_3d", "negative density or internal energy");
}
} while (redo);

// project back
#ifdef CUBED_SPHERE
for (int i = 0; i < w.GetDim1(); ++i) {
Expand Down

0 comments on commit 78d31ee

Please sign in to comment.