From 304eb2a1467330014b435ca9c415226a4236985d Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 18 May 2026 09:38:00 -0700 Subject: [PATCH 1/2] more precision changes --- src/embedded_boundaries/eb_annulus.cpp | 2 +- src/embedded_boundaries/eb_box.cpp | 18 +++++++++++------- src/embedded_boundaries/eb_cylinder.cpp | 2 +- src/embedded_boundaries/eb_sphere.cpp | 2 +- src/setup/init.cpp | 2 +- src/utilities/io.cpp | 4 ++-- 6 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/embedded_boundaries/eb_annulus.cpp b/src/embedded_boundaries/eb_annulus.cpp index f8c96a83a..b7eed7152 100644 --- a/src/embedded_boundaries/eb_annulus.cpp +++ b/src/embedded_boundaries/eb_annulus.cpp @@ -39,7 +39,7 @@ void incflo::make_eb_annulus() // Compute distance between cylinder centres Real offset = 0.0; for(int i = 0; i < AMREX_SPACEDIM; i++) - offset += std::pow(outer_center[i] - inner_center[i], 2); + offset += amrex::Math::powi<2>(outer_center[i] - inner_center[i]); offset = std::sqrt(offset); // Check that the inner cylinder is fully contained in the outer one diff --git a/src/embedded_boundaries/eb_box.cpp b/src/embedded_boundaries/eb_box.cpp index 3c361e1d8..9a372ae09 100644 --- a/src/embedded_boundaries/eb_box.cpp +++ b/src/embedded_boundaries/eb_box.cpp @@ -32,7 +32,11 @@ void incflo::make_eb_box() ************************************************************************/ Vector boxLo(AMREX_SPACEDIM), boxHi(AMREX_SPACEDIM); - Real offset = 1.0e-15; +#ifdef AMREX_USE_FLOAT + Real offset = Real(1e-8); +#else + Real offset = Real(1e-15); +#endif bool inside = true; for(int i = 0; i < AMREX_SPACEDIM; i++) @@ -54,8 +58,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(0)) { - xlo = 2.0 * geom[0].ProbLo(0) - geom[0].ProbHi(0); - xhi = 2.0 * geom[0].ProbHi(0) - geom[0].ProbLo(0); + xlo = Real(2) * geom[0].ProbLo(0) - geom[0].ProbHi(0); + xhi = Real(2) * geom[0].ProbHi(0) - geom[0].ProbLo(0); } Real ylo = boxLo[1] + offset; @@ -65,8 +69,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(1)) { - ylo = 2.0 * geom[0].ProbLo(1) - geom[0].ProbHi(1); - yhi = 2.0 * geom[0].ProbHi(1) - geom[0].ProbLo(1); + ylo = Real(2) * geom[0].ProbLo(1) - geom[0].ProbHi(1); + yhi = Real(2) * geom[0].ProbHi(1) - geom[0].ProbLo(1); } #if (AMREX_SPACEDIM > 2) @@ -77,8 +81,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(2)) { - zlo = 2.0 * geom[0].ProbLo(2) - geom[0].ProbHi(2); - zhi = 2.0 * geom[0].ProbHi(2) - geom[0].ProbLo(2); + zlo = Real(2) * geom[0].ProbLo(2) - geom[0].ProbHi(2); + zhi = Real(2) * geom[0].ProbHi(2) - geom[0].ProbLo(2); } #endif diff --git a/src/embedded_boundaries/eb_cylinder.cpp b/src/embedded_boundaries/eb_cylinder.cpp index 02afa9815..cc5ba9c03 100644 --- a/src/embedded_boundaries/eb_cylinder.cpp +++ b/src/embedded_boundaries/eb_cylinder.cpp @@ -33,7 +33,7 @@ void incflo::make_eb_cylinder() pp.getarr("center", centervec, 0, 3); Array center = {AMREX_D_DECL(centervec[0], centervec[1], centervec[2])}; - rotation = (rotation/180.)*M_PI; + rotation = (rotation/180.)*Real(M_PI); // Print info about cylinder amrex::Print() << " " << "\n"; diff --git a/src/embedded_boundaries/eb_sphere.cpp b/src/embedded_boundaries/eb_sphere.cpp index cdf735e06..09ef3563b 100644 --- a/src/embedded_boundaries/eb_sphere.cpp +++ b/src/embedded_boundaries/eb_sphere.cpp @@ -16,7 +16,7 @@ void incflo::make_eb_sphere() { // Initialise sphere parameters bool inside = true; - Real radius = 0.0002; + Real radius = Real(0.0002); Vector centervec(3); // Get sphere information from inputs file. * diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 2a3db0709..db7259667 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -222,7 +222,7 @@ void incflo::ReadParameters () amrex::Real tol_deg(0.); pp_eb_flow.query("normal_tol", tol_deg); - m_eb_flow.normal_tol = tol_deg*M_PI/amrex::Real(180.); + m_eb_flow.normal_tol = tol_deg*amrex::Real(M_PI)/amrex::Real(180); } if (m_advect_tracer && m_eb_flow.enabled && m_eb_flow.tracer.empty()) { diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index 801011e10..4f3096bb9 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -175,7 +175,7 @@ void incflo::ReadCheckpointFile() int i = 0; while(lis >> word) { - prob_lo[i++] = std::stod(word); // NOLINT(clang-analyzer-security.ArrayBound) + prob_lo[i++] = amrex::Real(std::stod(word)); // NOLINT(clang-analyzer-security.ArrayBound) } } @@ -186,7 +186,7 @@ void incflo::ReadCheckpointFile() int i = 0; while(lis >> word) { - prob_hi[i++] = std::stod(word); // NOLINT(clang-analyzer-security.ArrayBound) + prob_hi[i++] = amrex::Real(std::stod(word)); // NOLINT(clang-analyzer-security.ArrayBound) } } From a876fe8ad46e55502e2249962e3b54b1fb404509 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 21 Sep 2026 13:28:27 -0700 Subject: [PATCH 2/2] Fix issues #240-#249 - #240 ApplyCCProjection leaked cc_phi/cc_gphi on every call (raw new with no delete). Hold them in Vector and pass GetVecOfPtrs to Hydro. - #241 The single EB scalar op was shared by the tracer and temperature passes, and MLEBABecLap cannot revert an EB Dirichlet BC to Neumann, so an EB value set by one field was silently applied to the other. Give the temperature its own solve/apply ops, selected by an is_temperature flag. - #242 init_burggraf evaluated the exact solution at x-0.5, y-0.5 while the probtype-16 forcing, lid profile and DiffFromExact all use [0,1]^2 (regression from PR #160). Drop the shift. - #243 The 2D cylinder cull in incflo_PCInit rejected direction = 2 (the only value its disk formula fit) and culled a disk for directions 0/1, where AdvectWithFlow and EB2::CylinderIF use slabs. Use the same convention as AdvectWithFlow. - #244 A fresh start seeded tracer particles on level 0 and never redistributed, so initial particle tagging counted zero on level >= 1 and the first advection used level-0 MAC velocities. Redistribute after each level is created from scratch. - #245 steady_state = 1 was accepted but SteadyStateReached() is an Abort stub, so the run took one step and then died. Refuse it at read time, as PR #212 did for amr.KE_int. - #246 InitData never recorded m_last_chk (and restart set m_last_plt = 0), so Evolve's final-output block rewrote the step's checkpoint/plotfile, renaming the restart checkpoint to *.old.*. - #247 ReadCheckpointFile took finest_level from the header without checking amr.max_level, writing grids[]/m_leveldata[] out of bounds when restarting with a smaller max_level. Abort with a diagnostic instead. - #248 The update_temperature corrector passed the tracer laps (get_laps_new) to compute_laps_T instead of get_laps_tem_new. Latent today because MOL + temperature is refused at read time. - #249 The update_density corrector refreshed no ghost cells of density/density_nph, but ApplyCCProjection averages density_nph to faces through them, mixing predictor and corrector values on inter-grid faces. Use one ghost cell in both passes. Built clean (no warnings) for 2D/3D with and without EB. Co-Authored-By: Claude Opus 5 (1M context) --- src/diffusion/DiffusionScalarOp.H | 11 +- src/diffusion/DiffusionScalarOp.cpp | 109 +++++++++++------- src/diffusion/incflo_diffusion.cpp | 4 +- src/incflo.cpp | 11 +- src/incflo_update_density.cpp | 4 +- src/incflo_update_temperature.cpp | 2 +- src/particles/incflo_PCInit.cpp | 14 ++- src/particles/incflo_Tracers.cpp | 5 +- src/prob/prob_init_fluid.cpp | 6 +- src/projection/incflo_apply_cc_projection.cpp | 24 ++-- src/setup/init.cpp | 5 + src/utilities/io.cpp | 6 + 12 files changed, 131 insertions(+), 70 deletions(-) diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H index 7a17c69fd..54b47f92e 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 c78f31cae..bf5454feb 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 ad77063aa..d4d5cb5e4 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 ed47f4e5b..db497da18 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 3ee933944..6bee4c608 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 04efba638..dc2df40c6 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 1a28141e6..4ef53e3c7 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 81cd031ad..19af633d7 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 b2aef967b..5d1e79ca4 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 8b29c4635..5db553a39 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 6337cb7b6..b379805ab 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 7051030ca..e536cb1ee 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;