diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H index 7a17c69f..54b47f92 100644 --- a/src/diffusion/DiffusionScalarOp.H +++ b/src/diffusion/DiffusionScalarOp.H @@ -21,7 +21,8 @@ public: amrex::Vector const& eb_dirichlet, amrex::Vector const& use_chi, amrex::Vector bcrec, - amrex::Real dt); + amrex::Real dt, + bool is_temperature = false); void diffuse_vel_components ( amrex::Vector const& vel, @@ -33,7 +34,8 @@ public: amrex::Vector const& a_scalar, amrex::Vector const& eta, amrex::Vector const& eb_dirichlet, - amrex::Vector bcrec); + amrex::Vector bcrec, + bool is_temperature = false); void compute_divtau (amrex::Vector const& a_divtau, amrex::Vector const& a_vel, @@ -49,6 +51,11 @@ private: #ifdef AMREX_USE_EB std::unique_ptr m_eb_scal_solve_op; std::unique_ptr m_eb_scal_apply_op; + // The temperature gets its own EB ops: MLEBABecLap cannot revert an EB + // Dirichlet BC to Neumann, so the tracer and temperature passes must not + // share an operator when only one of them has an EB Dirichlet value. + std::unique_ptr m_eb_tem_solve_op; + std::unique_ptr m_eb_tem_apply_op; std::unique_ptr m_eb_vel_solve_op; std::unique_ptr m_eb_vel_apply_op; #endif diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp index c78f31ca..bf5454fe 100644 --- a/src/diffusion/DiffusionScalarOp.cpp +++ b/src/diffusion/DiffusionScalarOp.cpp @@ -33,6 +33,15 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo) info_solve, ebfact); m_eb_scal_solve_op->setMaxOrder(m_mg_maxorder); + if (m_incflo->m_use_temperature) + { + m_eb_tem_solve_op = std::make_unique(m_incflo->Geom(0,finest_level), + m_incflo->boxArray(0,finest_level), + m_incflo->DistributionMap(0,finest_level), + info_solve, ebfact); + m_eb_tem_solve_op->setMaxOrder(m_mg_maxorder); + } + if (!m_incflo->useTensorSolve()) { m_eb_vel_solve_op = std::make_unique(m_incflo->Geom(0,finest_level), @@ -51,6 +60,15 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo) m_incflo->DistributionMap(0,finest_level), info_apply, ebfact); m_eb_scal_apply_op->setMaxOrder(m_mg_maxorder); + + if (m_incflo->m_use_temperature) + { + m_eb_tem_apply_op = std::make_unique(m_incflo->Geom(0,finest_level), + m_incflo->boxArray(0,finest_level), + m_incflo->DistributionMap(0,finest_level), + info_apply, ebfact); + m_eb_tem_apply_op->setMaxOrder(m_mg_maxorder); + } } if ( (m_incflo->need_divtau() && !m_incflo->useTensorSolve()) || @@ -135,8 +153,10 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, Vector const& eb_dirichlet, amrex::Vector const& use_chi, amrex::Vector bcrec, - Real dt) + Real dt, + bool is_temperature) { + amrex::ignore_unused(is_temperature); // // Solves // [alpha a - beta div ( b grad )] sca = RHS @@ -168,6 +188,12 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, const int finest_level = m_incflo->finestLevel(); +#ifdef AMREX_USE_EB + // Tracer and temperature must not share an EB op (see the header) + MLEBABecLap* eb_op = is_temperature ? m_eb_tem_solve_op.get() + : m_eb_scal_solve_op.get(); +#endif + Vector rhs_c(finest_level+1); // Note only conservative uses this rhs_c container for (int lev = 0; lev <= finest_level; ++lev) { @@ -175,14 +201,14 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, } #ifdef AMREX_USE_EB - if (m_eb_scal_solve_op) + if (eb_op) { - m_eb_scal_solve_op->setScalars(1.0, dt); + eb_op->setScalars(1.0, dt); for (int lev = 0; lev <= finest_level; ++lev) { if ( use_chi[0] ) { - m_eb_scal_solve_op->setACoeffs(lev, *chi[lev]); + eb_op->setACoeffs(lev, *chi[lev]); } else { - m_eb_scal_solve_op->setACoeffs(lev, 1.0); + eb_op->setACoeffs(lev, 1.0); } } } @@ -202,16 +228,16 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, for (int comp = 0; comp < a_scalar[0]->nComp(); ++comp) { #ifdef AMREX_USE_EB - if (m_eb_scal_solve_op) + if (eb_op) { - m_eb_scal_solve_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, - bcrec[0].lo()), - m_incflo->get_diffuse_scalar_bc(Orientation::high, - bcrec[0].hi())); + eb_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, + bcrec[0].lo()), + m_incflo->get_diffuse_scalar_bc(Orientation::high, + bcrec[0].hi())); if ( m_incflo->m_has_mixedBC && comp>0 ) { // Must reset scalars (and Acoef, done below) to reuse solver with Robin BC - m_eb_scal_solve_op->setScalars(1.0, dt); + eb_op->setScalars(1.0, dt); } for (int lev = 0; lev <= finest_level; ++lev) { @@ -219,9 +245,9 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, if (comp > 0 && ( m_incflo->m_has_mixedBC || (use_chi[comp] != use_chi[comp-1]) ) ) { if ( use_chi[comp] ) { - m_eb_scal_solve_op->setACoeffs(lev, *chi[lev]); + eb_op->setACoeffs(lev, *chi[lev]); } else { - m_eb_scal_solve_op->setACoeffs(lev, 1.0); + eb_op->setACoeffs(lev, 1.0); } } @@ -230,12 +256,12 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, // single component of eta that we are solving for MultiFab phi (*eb_dirichlet[lev], amrex::make_alias, comp, 1); MultiFab beta(*eta[lev] , amrex::make_alias, comp, 1); - m_eb_scal_solve_op->setEBDirichlet(lev, phi, beta); + eb_op->setEBDirichlet(lev, phi, beta); } // else use default homogeneous Neumann on EB Array b = m_incflo->average_scalar_eta_to_faces(lev, comp, *eta[lev]); - m_eb_scal_solve_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid); + eb_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid); } } else @@ -286,19 +312,19 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, } #ifdef AMREX_USE_EB - if (m_eb_scal_solve_op) { + if (eb_op) { if ( m_incflo->m_has_mixedBC ) { auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]); - m_eb_scal_solve_op->setLevelBC(lev, &phi[lev], - &robin[0], &robin[1], &robin[2]); + eb_op->setLevelBC(lev, &phi[lev], + &robin[0], &robin[1], &robin[2]); } else { - m_eb_scal_solve_op->setLevelBC(lev, &phi[lev]); + eb_op->setLevelBC(lev, &phi[lev]); } // For when we use the stencil for centroid values - // m_eb_scal_solve_op->setPhiOnCentroid(); + // eb_op->setPhiOnCentroid(); } else #endif { @@ -307,7 +333,7 @@ DiffusionScalarOp::diffuse_scalar (Vector const& a_scalar, } #ifdef AMREX_USE_EB - MLMG mlmg(m_eb_scal_solve_op ? static_cast(*m_eb_scal_solve_op) : static_cast(*m_reg_scal_solve_op)); + MLMG mlmg(eb_op ? static_cast(*eb_op) : static_cast(*m_reg_scal_solve_op)); #else MLMG mlmg(*m_reg_scal_solve_op); #endif @@ -492,15 +518,20 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, Vector const& a_scalar, Vector const& a_eta, Vector const& eb_dirichlet, - amrex::Vector bcrec) + amrex::Vector bcrec, + bool is_temperature) { BL_PROFILE("DiffusionScalarOp::compute_laps"); + amrex::ignore_unused(is_temperature); int finest_level = m_incflo->finestLevel(); int n_comp = a_laps[0]->nComp(); #ifdef AMREX_USE_EB - if (m_eb_scal_apply_op) + // Tracer and temperature must not share an EB op (see the header) + MLEBABecLap* eb_op = is_temperature ? m_eb_tem_apply_op.get() + : m_eb_scal_apply_op.get(); + if (eb_op) { Vector laps_tmp(finest_level+1); int tmp_comp = (m_incflo->m_redistribution_type == "StateRedist") ? 3 : 2; @@ -514,22 +545,22 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, // We want to return div (mu grad)) phi - m_eb_scal_apply_op->setScalars(0.0, -1.0); + eb_op->setScalars(0.0, -1.0); // For when we use the stencil for centroid values - // m_eb_scal_apply_op->setPhiOnCentroid(); + // eb_op->setPhiOnCentroid(); for (int comp = 0; comp < n_comp; ++comp) { - m_eb_scal_apply_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, - bcrec[comp].lo()), - m_incflo->get_diffuse_scalar_bc(Orientation::high, - bcrec[comp].hi())); + eb_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low, + bcrec[comp].lo()), + m_incflo->get_diffuse_scalar_bc(Orientation::high, + bcrec[comp].hi())); int eta_comp = comp; if ( m_incflo->m_has_mixedBC && comp>0 ){ // Must reset scalars to reuse solver with Robin BC - m_eb_scal_apply_op->setScalars(0.0, -1.0); + eb_op->setScalars(0.0, -1.0); } Vector laps_comp; @@ -540,35 +571,31 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, // Use the same EB BC as the implicit solve in diffuse_scalar, // otherwise the explicit and implicit diffusion terms are - // inconsistent at the EB. NOTE that, exactly as for the solve - // op, this op is shared by the tracer and temperature passes and - // MLEBABecLap has no way to revert an EB Dirichlet BC back to - // Neumann, so if only one of tracer_eb/temperature_eb is defined - // the other pass will see the EB phi left by the first. + // inconsistent at the EB. if (lev < static_cast(eb_dirichlet.size()) && !eb_dirichlet[lev]->empty()) { MultiFab phi_eb (*eb_dirichlet[lev], amrex::make_alias, comp , 1); MultiFab beta_eb(*a_eta[lev] , amrex::make_alias, eta_comp, 1); - m_eb_scal_apply_op->setEBDirichlet(lev, phi_eb, beta_eb); + eb_op->setEBDirichlet(lev, phi_eb, beta_eb); } // else use default homogeneous Neumann on EB Array b = m_incflo->average_scalar_eta_to_faces(lev, eta_comp, *a_eta[lev]); - m_eb_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid); + eb_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid); if ( m_incflo->m_has_mixedBC ) { auto const robin = m_incflo->make_robinBC_MFs(lev, &scalar_comp[lev]); - m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev], - &robin[0], &robin[1], &robin[2]); + eb_op->setLevelBC(lev, &scalar_comp[lev], + &robin[0], &robin[1], &robin[2]); } else { - m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev]); + eb_op->setLevelBC(lev, &scalar_comp[lev]); } } - MLMG mlmg(*m_eb_scal_apply_op); + MLMG mlmg(*eb_op); mlmg.apply(GetVecOfPtrs(laps_comp), GetVecOfPtrs(scalar_comp)); } diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp index ad77063a..d4d5cb5e 100644 --- a/src/diffusion/incflo_diffusion.cpp +++ b/src/diffusion/incflo_diffusion.cpp @@ -80,7 +80,7 @@ incflo::compute_laps_T(Vector const& laps, Vector const& eta) { get_diffusion_scalar_op()->compute_laps(laps, scalar, eta, get_temperature_eb(), - get_temperature_bcrec()); + get_temperature_bcrec(), true); } void @@ -103,7 +103,7 @@ incflo::diffuse_temperature(Vector const& temperature, get_diffusion_scalar_op()->diffuse_scalar(temperature, rhocp, eta, get_temperature_eb(), {1} /* use rhocp */, - get_temperature_bcrec(), dt_diff); + get_temperature_bcrec(), dt_diff, true); } void diff --git a/src/incflo.cpp b/src/incflo.cpp index ed47f4e5..db497da1 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -79,7 +79,7 @@ void incflo::InitData () // xxxxx TODO averagedown ??? - if (m_check_int > 0) { WriteCheckPointFile(); } + if (m_check_int > 0) { WriteCheckPointFile(); m_last_chk = 0; } // Plot initial distribution if (m_plot_int > 0 || m_plot_per_exact > 0 || m_plot_per_approx > 0) @@ -98,6 +98,11 @@ void incflo::InitData () // Read starting configuration from chk file. ReadCheckpointFile(); + // The checkpoint just read is the checkpoint for this step; without this + // the final-output block in Evolve() rewrites it (renaming the original to + // *.old.*) when the restart does not take any steps. + m_last_chk = m_nstep; + #ifdef INCFLO_USE_PARTICLES particleData.Redistribute(); #endif @@ -105,12 +110,12 @@ void incflo::InitData () if (m_plotfile_on_restart) { WritePlotFile(); - m_last_plt = 0; + m_last_plt = m_nstep; } if (m_smallplotfile_on_restart) { WriteSmallPlotFile(); - m_last_smallplt = 0; + m_last_smallplt = m_nstep; } } diff --git a/src/incflo_update_density.cpp b/src/incflo_update_density.cpp index 3ee93394..6bee4c60 100644 --- a/src/incflo_update_density.cpp +++ b/src/incflo_update_density.cpp @@ -6,7 +6,9 @@ void incflo::update_density (StepType step_type) { BL_PROFILE("incflo::update_density"); - int ng = (step_type == StepType::Corrector) ? 0 : 1; + // One ghost cell in both passes: ApplyCCProjection averages density_nph to + // faces and reads its ghost cells, so the corrector must refresh them too. + const int ng = 1; Real l_dt = m_dt; diff --git a/src/incflo_update_temperature.cpp b/src/incflo_update_temperature.cpp index 04efba63..dc2df40c 100644 --- a/src/incflo_update_temperature.cpp +++ b/src/incflo_update_temperature.cpp @@ -23,7 +23,7 @@ void incflo::update_temperature (StepType step_type, Vector& tem_eta, { compute_temperature_diff_coeff(new_time, GetVecOfPtrs(tem_eta)); if (m_diff_type == DiffusionType::Explicit) { - compute_laps_T(get_laps_new(), get_temperature_new_const(), GetVecOfConstPtrs(tem_eta)); + compute_laps_T(get_laps_tem_new(), get_temperature_new_const(), GetVecOfConstPtrs(tem_eta)); } } diff --git a/src/particles/incflo_PCInit.cpp b/src/particles/incflo_PCInit.cpp index 1a28141e..4ef53e3c 100644 --- a/src/particles/incflo_PCInit.cpp +++ b/src/particles/incflo_PCInit.cpp @@ -238,11 +238,12 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par Real z_ctr = cyl_center[2]; #endif - // The cylinder axis matters here exactly as it does in AdvectWithFlow; - // assuming a z-parallel axis culls against the wrong axis for direction 0/1. + // The cylinder axis matters here exactly as it does in AdvectWithFlow, + // which accepts 0, 1 and 2 in both 2D and 3D; the culling below must + // use the same distance for each direction. int cyl_direction; pp.get("direction",cyl_direction); - AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction < AMREX_SPACEDIM); + AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= 2); // Remove particles that are outside of the cylinder for (ParIterType pti(*this, lev); pti.isValid(); ++pti) @@ -264,8 +265,11 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par : (cyl_direction == 1) ? std::sqrt(x*x + z*z) : std::sqrt(x*x + y*y); #else - // In 2D the only cylinder axis is the out-of-plane one - Real r = std::sqrt(x*x + y*y); + // Same convention as AdvectWithFlow and amrex::EB2::CylinderIF: + // direction 2 is a disk, directions 0 and 1 are slabs. + Real r = (cyl_direction == 2) ? std::sqrt(x*x + y*y) + : (cyl_direction == 1) ? std::abs(x) + : std::abs(y); #endif if (r > cyl_radius) { diff --git a/src/particles/incflo_Tracers.cpp b/src/particles/incflo_Tracers.cpp index 81cd031a..19af633d 100644 --- a/src/particles/incflo_Tracers.cpp +++ b/src/particles/incflo_Tracers.cpp @@ -50,7 +50,9 @@ void incflo::initializeTracerParticles ( ParGDBBase* a_gdb if (m_use_tracer_particles) namelist_unalloc.remove( incfloParticleNames::tracers ); - // If we've already allocated, then we probably want to resize + // If we've already allocated, then we probably want to resize, and the + // particles (seeded on level 0 only) must be moved onto the level that was + // just created; the restart path gets this from ParticleContainer::Restart. const auto& particles_namelist( particleData.getNames() ); for (auto it = particles_namelist.begin(); it != particles_namelist.end(); ++it) { std::string species_name( *it ); @@ -59,6 +61,7 @@ void incflo::initializeTracerParticles ( ParGDBBase* a_gdb if (!particleData[incfloParticleNames::tracers]->OK()) { particleData[incfloParticleNames::tracers]->resizeData(); } + particleData[incfloParticleNames::tracers]->Redistribute(); } } diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index b2aef967..5d1e79ca 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -1166,8 +1166,10 @@ void incflo::init_burggraf (Box const& vbx, Box const& /*gbx*/, { ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real x = (Real(i)+Real(0.5))*dx[0] - Real(0.5); - Real y = (Real(j)+Real(0.5))*dx[1] - Real(0.5); + // Same cell coordinates as the probtype-16 forcing (incflo_compute_forces.cpp), + // the lid profile (prob_bc.H) and DiffFromExact: the cavity is [0,1]^2. + Real x = (Real(i)+Real(0.5))*dx[0]; + Real y = (Real(j)+Real(0.5))*dx[1]; vel(i,j,k,0) = Real(8) * (x*x*x*x - Real(2) * x*x*x + x*x) * (Real(4)*y*y*y - Real(2)*y); vel(i,j,k,1) = -Real(8) * (Real(4)*x*x*x - Real(6) * x*x + Real(2)*x) * (y*y*y*y - y*y); #if (AMREX_SPACEDIM == 3) diff --git a/src/projection/incflo_apply_cc_projection.cpp b/src/projection/incflo_apply_cc_projection.cpp index 8b29c463..5db553a3 100644 --- a/src/projection/incflo_apply_cc_projection.cpp +++ b/src/projection/incflo_apply_cc_projection.cpp @@ -341,14 +341,14 @@ void incflo::ApplyCCProjection (Vector density, } } - Vector cc_phi(finest_level+1); - Vector cc_gphi(finest_level+1); + Vector cc_phi(finest_level+1); + Vector cc_gphi(finest_level+1); for (int lev = 0; lev <= finest_level; ++lev ) { - cc_phi[lev] = new MultiFab(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); - cc_gphi[lev] = new MultiFab(grids[lev], dmap[lev], AMREX_SPACEDIM, 0, MFInfo(), Factory(lev)); - cc_phi[lev]->setVal(0.); - cc_gphi[lev]->setVal(0.); + cc_phi[lev].define(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); + cc_gphi[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, 0, MFInfo(), Factory(lev)); + cc_phi[lev].setVal(0.); + cc_gphi[lev].setVal(0.); } Vector > mac_vec(finest_level+1); @@ -382,7 +382,7 @@ void incflo::ApplyCCProjection (Vector density, // // Perform MAC projection: - del dot (dt/rho) grad phi = div(U) // - macproj->project(cc_phi,m_mac_mg_rtol,m_mac_mg_atol); + macproj->project(GetVecOfPtrs(cc_phi),m_mac_mg_rtol,m_mac_mg_atol); // // After the projection we grab the dt/rho (grad phi) used in the projection @@ -403,9 +403,9 @@ void incflo::ApplyCCProjection (Vector density, // Note that "fluxes" comes back as MINUS (dt/rho) Gphi // #ifdef AMREX_USE_EB - macproj->getFluxes(amrex::GetVecOfArrOfPtrs(m_fluxes), cc_phi, MLMG::Location::FaceCentroid); + macproj->getFluxes(amrex::GetVecOfArrOfPtrs(m_fluxes), GetVecOfPtrs(cc_phi), MLMG::Location::FaceCentroid); #else - macproj->getFluxes(amrex::GetVecOfArrOfPtrs(m_fluxes), cc_phi, MLMG::Location::FaceCenter); + macproj->getFluxes(amrex::GetVecOfArrOfPtrs(m_fluxes), GetVecOfPtrs(cc_phi), MLMG::Location::FaceCenter); #endif for (int lev=0; lev <= finest_level; ++lev) @@ -413,7 +413,7 @@ void incflo::ApplyCCProjection (Vector density, #ifdef AMREX_USE_EB amrex::Abort("Haven't written mac_to_ccvel for EB"); #else - average_mac_to_ccvel(GetArrOfPtrs(m_fluxes[lev]),*cc_gphi[lev]); + average_mac_to_ccvel(GetArrOfPtrs(m_fluxes[lev]),cc_gphi[lev]); #endif } @@ -427,8 +427,8 @@ void incflo::ApplyCCProjection (Vector density, Box const& tbx = mfi.tilebox(); Array4 const& gp_cc = ld.gp.array(mfi); Array4 const& p_cc = ld.p_cc.array(mfi); - Array4 const& gphi = cc_gphi[lev]->const_array(mfi); - Array4 const& phi = cc_phi[lev]->const_array(mfi); + Array4 const& gphi = cc_gphi[lev].const_array(mfi); + Array4 const& phi = cc_phi[lev].const_array(mfi); Array4 const& u = ld.velocity.array(mfi); Array4 const& rho = density[lev]->const_array(mfi); diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 6337cb7b..b379805a 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -15,6 +15,11 @@ void incflo::ReadParameters () pp.query("stop_time", m_stop_time); pp.query("max_step", m_max_step); pp.query("steady_state", m_steady_state); + if (m_steady_state) { + // SteadyStateReached() is a stub that Aborts, so refuse the option up + // front instead of taking one step (with the stop_time cap disabled) first. + amrex::Abort("steady_state = 1: SteadyStateReached() is not implemented yet"); + } } { // Prefix amr diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index 7051030c..e536cb1e 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -149,6 +149,12 @@ void incflo::ReadCheckpointFile() // Finest level is >> finest_level; GotoNextLine(is); + if (finest_level > max_level) { + // grids/dmap/m_leveldata/m_t_new are sized max_level+1 + amrex::Abort("ReadCheckpointFile: checkpoint finest_level (" + + std::to_string(finest_level) + ") exceeds amr.max_level (" + + std::to_string(max_level) + ")"); + } // Step count is >> m_nstep;