From 311f0f861e7fa7fb5c81c57b04f3a2568eef9985 Mon Sep 17 00:00:00 2001 From: cgilet Date: Wed, 19 Aug 2026 14:30:49 -0400 Subject: [PATCH] Sync old/new-time state before ComputeDt to fix Boussinesq forcing and thus the computed dt on restart and regrid. ComputeDt calls compute_vel_forces_on_level, which uses both the old and new tracer fields for m_use_boussinesq. Previously, immediately after a restart (old-time fields aren't saved to checkpoints) or a regrid (new LevelData arrays start uninitialized) ComputeDt could read bad tracer_o data. Moving the fillpatch + copy-new-to-old block to the top of Advance and InitialIterations ensures old and new are both good before ComputeDt. Note this changes the adaptive dt for boussinesq. However, the dt is now computed with the same velocity forcing that is used to compute the advective velocity. This also ensures that copy_from_*_temperature checks m_use_temperature before trying to copy an mf that might not be defined. --- src/incflo_advance.cpp | 34 +++++++++++++++++----------------- src/incflo_utils.cpp | 2 ++ src/setup/init.cpp | 34 +++++++++++++++++----------------- 3 files changed, 36 insertions(+), 34 deletions(-) diff --git a/src/incflo_advance.cpp b/src/incflo_advance.cpp index 51f9099a..5344b76a 100644 --- a/src/incflo_advance.cpp +++ b/src/incflo_advance.cpp @@ -9,6 +9,23 @@ void incflo::Advance() // Start timing current time step Real strt_step = static_cast(ParallelDescriptor::second()); + int ng = nghost_state(); + for (int lev = 0; lev <= finest_level; ++lev) { + fillpatch_velocity(lev, m_cur_time, m_leveldata[lev]->velocity, ng); + fillpatch_density(lev, m_cur_time, m_leveldata[lev]->density, ng); + if (m_advect_tracer) { + fillpatch_tracer(lev, m_cur_time, m_leveldata[lev]->tracer, ng); + } + if (m_use_temperature) { + fillpatch_temperature(lev, m_cur_time, m_leveldata[lev]->temperature, ng); + } + } + + copy_from_new_to_old_velocity(IntVect(ng)); + copy_from_new_to_old_density(IntVect(ng)); + copy_from_new_to_old_tracer(IntVect(ng)); + copy_from_new_to_old_temperature(IntVect(ng)); + // Compute time step size int initialisation = ( m_dt < 0 ); bool explicit_diffusion = (m_diff_type == DiffusionType::Explicit); @@ -29,23 +46,6 @@ void incflo::Advance() << " with dt = " << m_dt << ".\n" << "\n"; } - copy_from_new_to_old_velocity(); - copy_from_new_to_old_density(); - copy_from_new_to_old_tracer(); - copy_from_new_to_old_temperature(); - - int ng = nghost_state(); - for (int lev = 0; lev <= finest_level; ++lev) { - fillpatch_velocity(lev, m_t_old[lev], m_leveldata[lev]->velocity_o, ng); - fillpatch_density(lev, m_t_old[lev], m_leveldata[lev]->density_o, ng); - if (m_advect_tracer) { - fillpatch_tracer(lev, m_t_old[lev], m_leveldata[lev]->tracer_o, ng); - } - if (m_use_temperature) { - fillpatch_temperature(lev, m_t_old[lev], m_leveldata[lev]->temperature_o, ng); - } - } - #ifdef AMREX_USE_EB if (m_eb_flow.enabled) { for (int lev = 0; lev <= finest_level; ++lev) { diff --git a/src/incflo_utils.cpp b/src/incflo_utils.cpp index c3cf793d..36a9c3ca 100644 --- a/src/incflo_utils.cpp +++ b/src/incflo_utils.cpp @@ -474,6 +474,7 @@ void incflo::copy_from_old_to_new_tracer (int lev, IntVect const& ng) void incflo::copy_from_new_to_old_temperature (IntVect const& ng) { + if (!m_use_temperature) { return; } for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy(m_leveldata[lev]->temperature_o, m_leveldata[lev]->temperature, 0, 0, 1, ng); @@ -482,6 +483,7 @@ void incflo::copy_from_new_to_old_temperature (IntVect const& ng) void incflo::copy_from_old_to_new_temperature (IntVect const& ng) { + if (!m_use_temperature) { return; } for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy(m_leveldata[lev]->temperature, m_leveldata[lev]->temperature_o, 0, 0, 1, ng); diff --git a/src/setup/init.cpp b/src/setup/init.cpp index db725966..dcc18a7d 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -437,10 +437,22 @@ void incflo::InitialIterations () { BL_PROFILE("incflo::InitialIterations()"); - copy_from_new_to_old_velocity(); - copy_from_new_to_old_density(); - copy_from_new_to_old_tracer(); - copy_from_new_to_old_temperature(); + int ng = nghost_state(); + for (int lev = 0; lev <= finest_level; ++lev) { + fillpatch_velocity(lev, m_cur_time, m_leveldata[lev]->velocity, ng); + fillpatch_density(lev, m_cur_time, m_leveldata[lev]->density, ng); + if (m_advect_tracer) { + fillpatch_tracer(lev, m_cur_time, m_leveldata[lev]->tracer, ng); + } + if (m_use_temperature) { + fillpatch_temperature(lev, m_cur_time, m_leveldata[lev]->temperature, ng); + } + } + + copy_from_new_to_old_velocity(IntVect(ng)); + copy_from_new_to_old_density(IntVect(ng)); + copy_from_new_to_old_tracer(IntVect(ng)); + copy_from_new_to_old_temperature(IntVect(ng)); int initialisation = 1; bool explicit_diffusion = (m_diff_type == DiffusionType::Explicit); @@ -456,23 +468,11 @@ void incflo::InitialIterations () for (int lev = 0; lev <= finest_level; ++lev) m_t_old[lev] = m_t_new[lev]; for (int lev = 0; lev <= finest_level; ++lev) mac_phi[lev]->setVal(0.); - int ng = nghost_state(); - for (int lev = 0; lev <= finest_level; ++lev) { - fillpatch_velocity(lev, m_t_old[lev], m_leveldata[lev]->velocity_o, ng); - fillpatch_density(lev, m_t_old[lev], m_leveldata[lev]->density_o, ng); - if (m_advect_tracer) { - fillpatch_tracer(lev, m_t_old[lev], m_leveldata[lev]->tracer_o, ng); - } - if (m_use_temperature) { - fillpatch_temperature(lev, m_t_old[lev], m_leveldata[lev]->temperature_o, ng); - } - } - for (int iter = 0; iter < m_initial_iterations; ++iter) { if (m_verbose) amrex::Print() << "\n In initial_iterations: iter = " << iter << "\n"; - ApplyPredictor(true); + ApplyPredictor(true); copy_from_old_to_new_velocity(); copy_from_old_to_new_density();