diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp index a40bb292..c010ce82 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -16,7 +16,7 @@ incflo::compute_MAC_projected_velocities ( Real /*time*/) { BL_PROFILE("incflo::compute_MAC_projected_velocities()"); - Real l_dt = m_dt; + Real l_dt = dt_real(); auto mac_phi = get_mac_phi(); @@ -60,7 +60,7 @@ incflo::compute_MAC_projected_velocities ( } // end m_godunov_include_diff_in_forcing if (nghost_force() > 0) { - fillpatch_force(m_cur_time, vel_forces, nghost_force()); + fillpatch_force(cur_time_real(), vel_forces, nghost_force()); } } // end m_advection_type @@ -198,7 +198,7 @@ incflo::compute_MAC_projected_velocities ( { MultiFab time_dep_inflow_vel(vel[lev]->boxArray(),vel[lev]->DistributionMap(),AMREX_SPACEDIM,1); time_dep_inflow_vel.setVal(0.); - fillphysbc_velocity(lev, m_cur_time+Real(0.5)*l_dt, time_dep_inflow_vel, 1); + fillphysbc_velocity(lev, real_time(m_cur_time + 0.5*m_dt), time_dep_inflow_vel, 1); Box domain(geom[lev].Domain()); const auto dlo = lbound(domain); diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 93808a00..77529372 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -187,7 +187,7 @@ incflo::compute_convective_term (Vector const& conv_u, } // end m_godunov_include_diff_in_forcing if (nghost_force() > 0) - fillpatch_force(m_cur_time, vel_forces, nghost_force()); + fillpatch_force(cur_time_real(), vel_forces, nghost_force()); // Note that for conservative tracers, this is forcing for (rho s) // and for non-conservative, this is forcing for s @@ -198,12 +198,12 @@ incflo::compute_convective_term (Vector const& conv_u, for (int lev = 0; lev <= finest_level; ++lev) MultiFab::Add(*tra_forces[lev], m_leveldata[lev]->laps_o, 0, 0, m_ntrac, 0); if (nghost_force() > 0) - fillpatch_force(m_cur_time, tra_forces, nghost_force()); + fillpatch_force(cur_time_real(), tra_forces, nghost_force()); } if (m_use_temperature) { - compute_tem_forces(m_cur_time, tem_forces); + compute_tem_forces(cur_time_real(), tem_forces); for (int lev = 0; lev <= finest_level; ++lev) { auto& ld = *m_leveldata[lev]; #ifdef _OPENMP @@ -233,14 +233,14 @@ incflo::compute_convective_term (Vector const& conv_u, } } if (nghost_force() > 0) { - fillpatch_force(m_cur_time, tem_forces, nghost_force()); } + fillpatch_force(cur_time_real(), tem_forces, nghost_force()); } } } // end m_advection_type for (int lev = 0; lev <= finest_level; ++lev) { - Real time_nph = m_cur_time + Real(0.5)*m_dt; + Real time_nph = real_time(m_cur_time + 0.5*m_dt); if (nghost_mac() > 0) { // FillPatch umac. @@ -672,7 +672,7 @@ incflo::compute_convective_term (Vector const& conv_u, : m_advect_momentum ? rhovel_f.const_array() : vel_forces[lev]->const_array(mfi), - geom[lev], m_dt, + geom[lev], dt_real(), get_velocity_bcrec(), get_velocity_bcrec_device_ptr(), get_velocity_iconserv_device_ptr(), @@ -710,7 +710,7 @@ incflo::compute_convective_term (Vector const& conv_u, v_mac[lev]->const_array(mfi), w_mac[lev]->const_array(mfi)), divu_arr, Array4{}, - geom[lev], m_dt, + geom[lev], dt_real(), get_density_bcrec(), get_density_bcrec_device_ptr(), get_density_iconserv_device_ptr(), @@ -778,7 +778,7 @@ incflo::compute_convective_term (Vector const& conv_u, w_mac[lev]->const_array(mfi)), divu_arr, (!tra_forces.empty()) ? tra_forces[lev]->const_array(mfi) : Array4{}, - geom[lev], m_dt, + geom[lev], dt_real(), get_tracer_bcrec(), get_tracer_bcrec_device_ptr(), get_tracer_iconserv_device_ptr(), @@ -820,7 +820,7 @@ incflo::compute_convective_term (Vector const& conv_u, w_mac[lev]->const_array(mfi)), divu_arr, (!tem_forces.empty()) ? tem_forces[lev]->const_array(mfi) : Array4{}, - geom[lev], m_dt, + geom[lev], dt_real(), get_temperature_bcrec(), get_temperature_bcrec_device_ptr(), m_iconserv_temperature_d.data(), diff --git a/src/incflo.H b/src/incflo.H index 3950fa87..eb90677b 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -405,27 +405,34 @@ private: amrex::Vector m_t_old; amrex::Vector m_t_new; - // Times - amrex::Real m_cur_time = amrex::Real( 0.0); - amrex::Real m_dt = amrex::Real(-1.0); - amrex::Real m_prev_dt = amrex::Real(-1.0); - amrex::Real m_prev_prev_dt = amrex::Real(-1.0); + // Keep time-step bookkeeping in double precision even in single-precision builds. + double m_cur_time = 0.0; + double m_dt = -1.0; + double m_prev_dt = -1.0; + double m_prev_prev_dt = -1.0; + + [[nodiscard]] static amrex::Real real_time (double a_time) noexcept { + return static_cast(a_time); + } + [[nodiscard]] amrex::Real cur_time_real () const noexcept { return real_time(m_cur_time); } + [[nodiscard]] amrex::Real dt_real () const noexcept { return real_time(m_dt); } + [[nodiscard]] amrex::Real new_time_real () const noexcept { return real_time(m_cur_time + m_dt); } // Time step counter int m_nstep = -1; // Stop simulation if cur_time reaches stop_time OR nstep reaches max_step // OR steady_state = true AND steady_state_tol is reached - amrex::Real m_stop_time = amrex::Real(-1.0); + double m_stop_time = -1.0; int m_max_step = -1; bool m_steady_state = false; amrex::Real m_steady_state_tol = amrex::Real(1.0e-5); // Options to control time stepping - amrex::Real m_cfl = amrex::Real(0.5); - amrex::Real m_fixed_dt = amrex::Real(-1.); - amrex::Real m_init_shrink = amrex::Real(0.1); - amrex::Real m_dt_change_max = amrex::Real(1.1); + double m_cfl = 0.5; + double m_fixed_dt = -1.0; + double m_init_shrink = 0.1; + double m_dt_change_max = 1.1; // Initial projection / iterations bool m_do_initial_proj = true; @@ -600,16 +607,16 @@ private: amrex::Real m_papa_reg = 0.0; //! Lower bound on the strain rate used by the power-law viscosity, so that //! eta = mu*sr^(n-1) stays finite (and non-zero) at zero strain rate - amrex::Real m_sr_floor = 1.e-9; + amrex::Real m_sr_floor = amrex::Real(1.e-9); amrex::Real m_eta_0 = 0.0; int m_plot_int = -1; // Dump plotfiles at as close as possible to the designated period *without* changing dt - amrex::Real m_plot_per_approx = -1.0; + double m_plot_per_approx = -1.0; // Dump plotfiles at exactcly the designated period by changing dt - amrex::Real m_plot_per_exact = -1.0; + double m_plot_per_exact = -1.0; int m_last_plt = -1; std::string m_plot_file{"plt"}; @@ -638,7 +645,7 @@ private: // smallplotfile allows users to output certain variables at a different frequency. // Off by default. int m_smallplot_int = -1; - amrex::Real m_smallplot_per_approx = -1.0; + double m_smallplot_per_approx = -1.0; // No smallplot_per_exact because it requires changing the timestep int m_last_smallplt = -1; bool m_smallplotfile_on_restart = false; @@ -983,8 +990,8 @@ private: void Advance (); [[nodiscard]] bool writeNow () { return writeNow(m_plot_int, m_plot_per_approx, m_plot_per_exact); } - [[nodiscard]] bool writeNow (int a_plot_int, amrex::Real a_plot_per_approx, - amrex::Real a_plot_per_exact) const; + [[nodiscard]] bool writeNow (int a_plot_int, double a_plot_per_approx, + double a_plot_per_exact) const; /////////////////////////////////////////////////////////////////////////// // diff --git a/src/incflo.cpp b/src/incflo.cpp index 8122dff1..03b68410 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -53,7 +53,7 @@ void incflo::InitData () // This is an AmrCore member function which recursively makes new levels // with MakeNewLevelFromScratch. - InitFromScratch(m_cur_time); + InitFromScratch(cur_time_real()); #ifdef AMREX_USE_EB if (!EBFactory(0).isAllRegular()) { @@ -147,7 +147,7 @@ void incflo::Evolve() if (m_regrid_int > 0 && m_nstep > 0 && m_nstep%m_regrid_int == 0) { if (m_verbose > 0) amrex::Print() << "Regridding...\n"; - regrid(0, m_cur_time); + regrid(0, cur_time_real()); if (m_verbose > 0 && ParallelDescriptor::IOProcessor()) { printGridSummary(amrex::OutStream(), 0, finest_level); } @@ -289,7 +289,7 @@ void incflo::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_gr } bool -incflo::writeNow(int a_plot_int, Real a_plot_per_approx, Real a_plot_per_exact) const +incflo::writeNow(int a_plot_int, double a_plot_per_approx, double a_plot_per_exact) const { bool write_now = false; @@ -312,8 +312,8 @@ incflo::writeNow(int a_plot_int, Real a_plot_per_approx, Real a_plot_per_exact) // the counter, because we have indeed reached the next a_plot_per_approx interval // at this point. - const Real eps = std::numeric_limits::epsilon() * Real(10.0) * std::abs(m_cur_time); - const Real next_plot_time = (num_per_old + 1) * a_plot_per_approx; + const double eps = std::numeric_limits::epsilon() * 10.0 * std::abs(m_cur_time); + const double next_plot_time = (num_per_old + 1) * a_plot_per_approx; if ((num_per_new == num_per_old) && std::abs(m_cur_time - next_plot_time) <= eps) { diff --git a/src/incflo_advance.cpp b/src/incflo_advance.cpp index 0cab48ca..1f36d3be 100644 --- a/src/incflo_advance.cpp +++ b/src/incflo_advance.cpp @@ -11,13 +11,13 @@ void incflo::Advance() 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); + fillpatch_velocity(lev, cur_time_real(), m_leveldata[lev]->velocity, ng); + fillpatch_density(lev, cur_time_real(), m_leveldata[lev]->density, ng); if (m_advect_tracer) { - fillpatch_tracer(lev, m_cur_time, m_leveldata[lev]->tracer, ng); + fillpatch_tracer(lev, cur_time_real(), m_leveldata[lev]->tracer, ng); } if (m_use_temperature) { - fillpatch_temperature(lev, m_cur_time, m_leveldata[lev]->temperature, ng); + fillpatch_temperature(lev, cur_time_real(), m_leveldata[lev]->temperature, ng); } } @@ -34,8 +34,8 @@ void incflo::Advance() // Set new and old time to correctly use in fillpatching for(int lev = 0; lev <= finest_level; lev++) { - m_t_old[lev] = m_cur_time; - m_t_new[lev] = m_cur_time + m_dt; + m_t_old[lev] = cur_time_real(); + m_t_new[lev] = new_time_real(); } if (m_verbose > 0) diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index 4d00cae4..a623136d 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -70,7 +70,7 @@ void incflo::ApplyCorrector() BL_PROFILE("incflo::ApplyCorrector"); // We use the new time value for things computed on the "*" state - Real new_time = m_cur_time + m_dt; + Real new_time = new_time_real(); // ************************************************************************************* // Allocate space for the MAC velocities @@ -163,7 +163,7 @@ void incflo::ApplyCorrector() bool incremental_projection = false; ApplyProjection(get_density_nph_const(), AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), - GetVecOfPtrs(w_mac)),new_time,m_dt,incremental_projection); + GetVecOfPtrs(w_mac)), new_time, dt_real(), incremental_projection); #ifdef AMREX_USE_EB // ********************************************************************************************** diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index 114e3a7e..ae56443d 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -79,7 +79,7 @@ void incflo::ApplyPredictor (bool incremental_projection) BL_PROFILE("incflo::ApplyPredictor"); // We use the new time value for things computed on the "*" state - Real new_time = m_cur_time + m_dt; + Real new_time = new_time_real(); // ************************************************************************************* // Allocate space for the MAC velocities @@ -137,7 +137,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* compute_viscosity(GetVecOfPtrs(vel_eta), get_density_old(), get_velocity_old(), - m_cur_time, nghost_eta); + cur_time_real(), nghost_eta); // ************************************************************************************* // Compute explicit viscous term @@ -161,7 +161,7 @@ void incflo::ApplyPredictor (bool incremental_projection) } if (m_use_temperature) { - compute_temperature_diff_coeff(m_cur_time, GetVecOfPtrs(tem_eta)); + compute_temperature_diff_coeff(cur_time_real(), GetVecOfPtrs(tem_eta)); if (need_divtau()) { compute_laps_T(get_laps_tem_old(), get_temperature_old_const(), GetVecOfConstPtrs(tem_eta)); } @@ -180,7 +180,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* compute_MAC_projected_velocities(get_velocity_old_const(), get_density_old_const(), AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), - GetVecOfPtrs(w_mac)), GetVecOfPtrs(vel_forces), m_cur_time); + GetVecOfPtrs(w_mac)), GetVecOfPtrs(vel_forces), cur_time_real()); // ************************************************************************************* // if (advection_type == "Godunov") @@ -196,7 +196,7 @@ void incflo::ApplyPredictor (bool incremental_projection) AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), GetVecOfPtrs(w_mac)), GetVecOfPtrs(vel_forces), GetVecOfPtrs(tra_forces), - GetVecOfPtrs(tem_forces), m_cur_time); + GetVecOfPtrs(tem_forces), cur_time_real()); // ************************************************************************************* // Update density @@ -223,7 +223,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // ********************************************************************************************** ApplyProjection(get_density_nph_const(), AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), - GetVecOfPtrs(w_mac)),new_time,m_dt,incremental_projection); + GetVecOfPtrs(w_mac)), new_time, dt_real(), incremental_projection); #ifdef INCFLO_USE_PARTICLES // ************************************************************************************** diff --git a/src/incflo_compute_dt.cpp b/src/incflo_compute_dt.cpp index 960b4556..30dea77c 100644 --- a/src/incflo_compute_dt.cpp +++ b/src/incflo_compute_dt.cpp @@ -1,5 +1,6 @@ #include +#include #include #include @@ -194,28 +195,30 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion) ParallelContext::CommunicatorSub()); // Combined CFL conditioner - Real comb_cfl = cd_cfl + std::sqrt(cd_cfl*cd_cfl + Real(4.0) * forc_cfl); + double const cd_cfl_d = static_cast(cd_cfl); + double const forc_cfl_d = static_cast(forc_cfl); + double const comb_cfl = cd_cfl_d + std::sqrt(cd_cfl_d*cd_cfl_d + 4.0 * forc_cfl_d); // Update dt - Real dt_new; - if (comb_cfl > 0.) + double dt_new; + if (comb_cfl > 0.0) { - dt_new = Real(2.0) * m_cfl / comb_cfl; + dt_new = 2.0 * m_cfl / comb_cfl; } else { // This is totally random but just a way to set a timestep // when the initial velocity is zero and the forcing term // is not a body force - auto const dx = geom[finest_level].CellSizeArray(); - dt_new = std::min(dx[0],dx[1]); + auto const dx = geom[finest_level].CellSizeArray(); + dt_new = static_cast(std::min(dx[0], dx[1])); #if (AMREX_SPACEDIM == 3) - dt_new = std::min(dt_new,dx[2]); + dt_new = std::min(dt_new, static_cast(dx[2])); #endif } // Optionally reduce CFL for initial step - if(initialization) + if (initialization) { dt_new *= m_init_shrink; } @@ -223,59 +226,61 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion) // Protect against very small comb_cfl // This may happen, for example, when the initial velocity field // is zero for an inviscid flow with no external forcing - Real eps = std::numeric_limits::epsilon(); - if(! initialization && comb_cfl <= eps) + double const real_eps = static_cast(std::numeric_limits::epsilon()); + double const time_eps = std::numeric_limits::epsilon(); + if (!initialization && comb_cfl <= real_eps) { - dt_new = Real(0.5) * m_dt; + dt_new = 0.5 * m_dt; } // Don't let the timestep grow by more than m_dt_change_max per step // unless the previous time step was unduly shrunk to match m_plot_per_exact - Real allowed_change_factor = m_dt_change_max; - if( (m_dt > Real(0.0)) && !(m_plot_per_exact > 0 && m_last_plt == m_nstep && m_nstep > 0) ) + double const allowed_change_factor = m_dt_change_max; + if ((m_dt > 0.0) && !(m_plot_per_exact > 0.0 && m_last_plt == m_nstep && m_nstep > 0)) { - dt_new = amrex::min(dt_new, allowed_change_factor * m_prev_dt); + dt_new = std::min(dt_new, allowed_change_factor * m_prev_dt); } - else if ( (m_dt > Real(0.0)) && (m_plot_per_exact > 0 && m_last_plt == m_nstep && m_nstep > 0) ) + else if ((m_dt > 0.0) && (m_plot_per_exact > 0.0 && m_last_plt == m_nstep && m_nstep > 0)) { - dt_new = amrex::min( dt_new, allowed_change_factor * amrex::max(m_prev_dt, m_prev_prev_dt) ); + dt_new = std::min(dt_new, allowed_change_factor * std::max(m_prev_dt, m_prev_prev_dt)); } // Don't overshoot specified plot times - if(m_plot_per_exact > Real(0.0) && - (std::trunc((m_cur_time + dt_new + eps) / m_plot_per_exact) > std::trunc((m_cur_time + eps) / m_plot_per_exact))) + if (m_plot_per_exact > 0.0 && + (std::trunc((m_cur_time + dt_new + time_eps) / m_plot_per_exact) > + std::trunc((m_cur_time + time_eps) / m_plot_per_exact))) { dt_new = std::trunc((m_cur_time + dt_new) / m_plot_per_exact) * m_plot_per_exact - m_cur_time; } // Don't overshoot the final time if not running to steady state - if(!m_steady_state && m_stop_time > Real(0.0)) + if (!m_steady_state && m_stop_time > 0.0) { - if(m_cur_time + dt_new > m_stop_time) + if (m_cur_time + dt_new > m_stop_time) { dt_new = m_stop_time - m_cur_time; } } // Make sure the timestep is not set to zero after a m_plot_per_exact stop - if (dt_new < eps) + if (dt_new < time_eps) { - dt_new = Real(0.5) * m_dt; + dt_new = 0.5 * m_dt; } // If using fixed time step, check CFL condition and give warning if not satisfied - if (m_fixed_dt > Real(0.0)) + if (m_fixed_dt > 0.0) { - if(dt_new < m_fixed_dt) - { - amrex::Print() << "WARNING: fixed_dt does not satisfy CFL condition: \n" - << "max dt by CFL : " << dt_new << "\n" - << "fixed dt specified: " << m_fixed_dt << "\n"; - } - m_dt = m_fixed_dt; + if (dt_new < m_fixed_dt) + { + amrex::Print() << "WARNING: fixed_dt does not satisfy CFL condition: \n" + << "max dt by CFL : " << dt_new << "\n" + << "fixed dt specified: " << m_fixed_dt << "\n"; + } + m_dt = m_fixed_dt; } else { - m_dt = dt_new; + m_dt = dt_new; } } diff --git a/src/incflo_explicit_update.cpp b/src/incflo_explicit_update.cpp index bb636e61..3e4f57af 100644 --- a/src/incflo_explicit_update.cpp +++ b/src/incflo_explicit_update.cpp @@ -7,7 +7,7 @@ void incflo::tracer_explicit_update (Vector const& tra_forces) if (m_advect_tracer == 0) { return; } constexpr Real m_half = Real(0.5); - Real l_dt = m_dt; + Real l_dt = dt_real(); int l_ntrac = m_ntrac; for (int lev = 0; lev <= finest_level; lev++) { @@ -98,7 +98,7 @@ void incflo::tracer_explicit_update_corrector (Vector const& tra_force if (m_advect_tracer == 0) { return; } constexpr Real m_half = Real(0.5); - Real l_dt = m_dt; + Real l_dt = dt_real(); int l_ntrac = m_ntrac; for (int lev = 0; lev <= finest_level; lev++) { diff --git a/src/incflo_redistribute.cpp b/src/incflo_redistribute.cpp index 34391962..5f6b3e80 100644 --- a/src/incflo_redistribute.cpp +++ b/src/incflo_redistribute.cpp @@ -90,7 +90,7 @@ incflo::redistribute_term ( MFIter const& mfi, scratch, flag, AMREX_D_DECL(apx, apy, apz), vfrac, AMREX_D_DECL(fcx, fcy, fcz), ccc, - bc, geom[lev], m_dt, m_redistribution_type); + bc, geom[lev], dt_real(), m_redistribution_type); } else { diff --git a/src/incflo_update_density.cpp b/src/incflo_update_density.cpp index 3ee93394..d377e5e3 100644 --- a/src/incflo_update_density.cpp +++ b/src/incflo_update_density.cpp @@ -8,7 +8,7 @@ void incflo::update_density (StepType step_type) int ng = (step_type == StepType::Corrector) ? 0 : 1; - Real l_dt = m_dt; + Real l_dt = dt_real(); if (!m_constant_density) { diff --git a/src/incflo_update_temperature.cpp b/src/incflo_update_temperature.cpp index 04efba63..cbafe972 100644 --- a/src/incflo_update_temperature.cpp +++ b/src/incflo_update_temperature.cpp @@ -8,8 +8,8 @@ void incflo::update_temperature (StepType step_type, Vector& tem_eta, if (!m_use_temperature) { return; } - Real const new_time = m_cur_time + m_dt; - Real const half_time = m_cur_time + m_dt * Real(0.5); + Real const new_time = new_time_real(); + Real const half_time = real_time(m_cur_time + 0.5*m_dt); // ************************************************************************************* // Compute the temperature forcing terms @@ -32,7 +32,7 @@ void incflo::update_temperature (StepType step_type, Vector& tem_eta, // ************************************************************************************* if (step_type == StepType::Predictor) { constexpr Real m_half = Real(0.5); - Real l_dt = m_dt; + Real l_dt = dt_real(); for (int lev = 0; lev <= finest_level; lev++) { @@ -107,7 +107,7 @@ void incflo::update_temperature (StepType step_type, Vector& tem_eta, m_leveldata[lev]->temperature.FillBoundary(geom[lev].periodicity()); fillphysbc_temperature(lev, new_time, m_leveldata[lev]->temperature, ng_diffusion); } - Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : Real(0.5)*m_dt; + Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? dt_real() : real_time(0.5*m_dt); // scratch holds rhoCp diffuse_temperature(get_temperature_new(), GetVecOfPtrs(scratch), GetVecOfConstPtrs(tem_eta), dt_diff); diff --git a/src/incflo_update_tracer.cpp b/src/incflo_update_tracer.cpp index bd547402..5ecb0d0a 100644 --- a/src/incflo_update_tracer.cpp +++ b/src/incflo_update_tracer.cpp @@ -6,7 +6,7 @@ void incflo::update_tracer (StepType step_type, Vector& tra_eta, Vecto { BL_PROFILE("incflo::update_tracer"); - Real new_time = m_cur_time + m_dt; + Real new_time = new_time_real(); if (m_advect_tracer) { @@ -53,7 +53,7 @@ void incflo::update_tracer (StepType step_type, Vector& tra_eta, Vecto fillphysbc_tracer(lev, new_time, m_leveldata[lev]->tracer, ng_diffusion); } - Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : Real(0.5)*m_dt; + Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? dt_real() : real_time(0.5*m_dt); diffuse_scalar(get_tracer_new(), get_density_new(), GetVecOfConstPtrs(tra_eta), dt_diff); } else diff --git a/src/incflo_update_velocity.cpp b/src/incflo_update_velocity.cpp index f4457e42..90f1d6d2 100644 --- a/src/incflo_update_velocity.cpp +++ b/src/incflo_update_velocity.cpp @@ -6,9 +6,9 @@ void incflo::update_velocity (StepType step_type, Vector& vel_eta, Vec { BL_PROFILE("incflo::update_velocity"); - Real new_time = m_cur_time + m_dt; + Real new_time = new_time_real(); - Real l_dt = m_dt; + Real l_dt = dt_real(); Real l_half = Real(0.5); if (step_type == StepType::Predictor) { @@ -337,7 +337,7 @@ void incflo::update_velocity (StepType step_type, Vector& vel_eta, Vec fillphysbc_density (lev, new_time, m_leveldata[lev]->density , ng_diffusion); } - Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : l_half*m_dt; + Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? dt_real() : real_time(0.5*m_dt); diffuse_velocity(get_velocity_new(), get_density_new(), GetVecOfConstPtrs(vel_eta), dt_diff); } } diff --git a/src/particles/incflo_Tracers.cpp b/src/particles/incflo_Tracers.cpp index 81cd031a..a5dd519e 100644 --- a/src/particles/incflo_Tracers.cpp +++ b/src/particles/incflo_Tracers.cpp @@ -73,7 +73,7 @@ void incflo::evolveTracerParticles (AMREX_D_DECL(Vector const& if (m_use_tracer_particles) { for (int lev = 0; lev <= finest_level; ++lev) { - particleData[incfloParticleNames::tracers]->EvolveParticles(lev, m_dt, + particleData[incfloParticleNames::tracers]->EvolveParticles(lev, dt_real(), AMREX_D_DECL(u_mac[lev],v_mac[lev],w_mac[lev])); } diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 5470993c..1edfc4f8 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -57,7 +57,7 @@ void incflo::ReadParameters () // This limits dt growth per time step pp.query("dt_change_max", m_dt_change_max); - if ( m_dt_change_max < 1.0_rt || m_dt_change_max > 1.1_rt ) { + if ( m_dt_change_max < 1.0 || m_dt_change_max > 1.1 ) { amrex::Abort("We require 1. <= dt_change_max <= 1.1"); } @@ -452,13 +452,13 @@ void incflo::InitialIterations () 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); + fillpatch_velocity(lev, cur_time_real(), m_leveldata[lev]->velocity, ng); + fillpatch_density(lev, cur_time_real(), m_leveldata[lev]->density, ng); if (m_advect_tracer) { - fillpatch_tracer(lev, m_cur_time, m_leveldata[lev]->tracer, ng); + fillpatch_tracer(lev, cur_time_real(), m_leveldata[lev]->tracer, ng); } if (m_use_temperature) { - fillpatch_temperature(lev, m_cur_time, m_leveldata[lev]->temperature, ng); + fillpatch_temperature(lev, cur_time_real(), m_leveldata[lev]->temperature, ng); } } @@ -494,8 +494,8 @@ void incflo::InitialIterations () } // Reset dt to get initial step as specified, otherwise we can see increase to dt - m_prev_dt = Real(-1.0); - m_dt = Real(-1.0); + m_prev_dt = -1.0; + m_dt = -1.0; } // Project velocity field to make sure initial velocity is divergence-free @@ -532,7 +532,7 @@ void incflo::InitialProjection() ApplyProjection(get_density_new_const(), AMREX_D_DECL(GetVecOfPtrs(u_mac_tmp), GetVecOfPtrs(v_mac_tmp), - GetVecOfPtrs(w_mac_tmp)),m_cur_time,dummy_dt,incremental_projection); + GetVecOfPtrs(w_mac_tmp)), cur_time_real(), dummy_dt, incremental_projection); // We set p and gp back to zero (p0 may still be still non-zero) @@ -621,7 +621,7 @@ void incflo::InitialPressureProjection() // FIXME FIXME FIXME - THIS ONLY WORKS RIGHT FOR NODAL PROJ ApplyProjection(get_density_new_const(), GetVecOfPtrs(vel), Source, - m_cur_time, dummy_dt, false /*incremental*/, + cur_time_real(), dummy_dt, false /*incremental*/, true /*set_inflow_bc*/); } diff --git a/src/utilities/incflo_steady_state.cpp b/src/utilities/incflo_steady_state.cpp index d0353569..335d671e 100644 --- a/src/utilities/incflo_steady_state.cpp +++ b/src/utilities/incflo_steady_state.cpp @@ -27,7 +27,7 @@ bool incflo::SteadyStateReached() int condition2[finest_level + 1]; // Make sure velocity is up to date - incflo_set_velocity_bcs(m_cur_time, vel); + incflo_set_velocity_bcs(cur_time_real(), vel); // Use temporaries to store the difference between current and previous solution Vector> diff_vel; diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index 7051030c..19f51815 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -208,7 +208,7 @@ void incflo::ReadCheckpointFile() // Create distribution mapping DistributionMapping dm{ba, ParallelDescriptor::NProcs()}; - MakeNewLevelFromScratch(lev, m_cur_time, ba, dm); + MakeNewLevelFromScratch(lev, cur_time_real(), ba, dm); } /*************************************************************************** @@ -251,7 +251,7 @@ void incflo::ReadCheckpointFile() #endif if ( m_regrid_on_restart ) { - regrid(0, m_cur_time); + regrid(0, cur_time_real()); } amrex::Print() << "Restart complete" << "\n"; @@ -391,9 +391,9 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo #else const int ng = 1; #endif - fillpatch_velocity(lev, m_cur_time, m_leveldata[lev]->velocity, ng); - fillpatch_density(lev, m_cur_time, m_leveldata[lev]->density, ng); - fillpatch_tracer(lev, m_cur_time, m_leveldata[lev]->tracer, ng); + fillpatch_velocity(lev, cur_time_real(), m_leveldata[lev]->velocity, ng); + fillpatch_density(lev, cur_time_real(), m_leveldata[lev]->density, ng); + fillpatch_tracer(lev, cur_time_real(), m_leveldata[lev]->tracer, ng); // Whether temperature fillpatch is needed depends on form of forcing term // fillpatch_temperature(lev, m_cur_time, m_leveldata[lev]->temperature, ng); } @@ -538,7 +538,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy(mf[lev], m_leveldata[lev]->velocity, 0, icomp, 1, 0); - DiffFromExact(lev, Geom(lev), m_cur_time, m_dt, mf[lev], icomp, icomp_err_u); + DiffFromExact(lev, Geom(lev), cur_time_real(), dt_real(), mf[lev], icomp, icomp_err_u); amrex::Print() << "Norm0 / Norm2 of u error " << mf[lev].norm0(icomp) << " " << mf[lev].norm2(icomp) / std::sqrt(mf[lev].boxArray().numPts()) << "\n"; } @@ -550,7 +550,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy(mf[lev], m_leveldata[lev]->velocity, 1, icomp, 1, 0); - DiffFromExact(lev, Geom(lev), m_cur_time, m_dt, mf[lev], icomp, icomp_err_v); + DiffFromExact(lev, Geom(lev), cur_time_real(), dt_real(), mf[lev], icomp, icomp_err_v); amrex::Print() << "Norm0 / Norm2 of v error " << mf[lev].norm0(icomp) << " " << mf[lev].norm2(icomp) / std::sqrt(mf[lev].boxArray().numPts()) << "\n"; } @@ -563,7 +563,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy(mf[lev], m_leveldata[lev]->velocity, 2, icomp, 1, 0); - DiffFromExact(lev, Geom(lev), m_cur_time, m_dt, mf[lev], icomp, icomp_err_w); + DiffFromExact(lev, Geom(lev), cur_time_real(), dt_real(), mf[lev], icomp, icomp_err_w); amrex::Print() << "Norm0 / Norm2 of w error " << mf[lev].norm0(icomp) << " " << mf[lev].norm2(icomp) / std::sqrt(mf[lev].boxArray().numPts()) << "\n"; } @@ -589,7 +589,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo for (int lev = 0; lev <= finest_level; ++lev) { mf[lev].plus(-offset, icomp, 1); - DiffFromExact(lev, Geom(lev), m_cur_time, m_dt, mf[lev], icomp, icomp_err_p); + DiffFromExact(lev, Geom(lev), cur_time_real(), dt_real(), mf[lev], icomp, icomp_err_p); amrex::Print() << "Norm0 / Norm2 of p error " << mf[lev].norm0(icomp) << " " << mf[lev].norm2(icomp) / std::sqrt(mf[lev].boxArray().numPts()) << "\n"; } @@ -608,7 +608,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo for (int lev = 0; lev <= finest_level; ++lev) { mf[lev].plus(-offset, icomp, 1); - DiffFromExact(lev, Geom(lev), m_cur_time, m_dt, mf[lev], icomp, icomp_err_mac_p); + DiffFromExact(lev, Geom(lev), cur_time_real(), dt_real(), mf[lev], icomp, icomp_err_mac_p); amrex::Print() << "Norm0 / Norm2 of mac_p error " << mf[lev].norm0(icomp) << " " << mf[lev].norm2(icomp) / std::sqrt(mf[lev].boxArray().numPts()) << "\n"; } @@ -623,7 +623,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo &m_leveldata[lev]->density, &m_leveldata[lev]->velocity, Geom(lev), - m_cur_time, 0); + cur_time_real(), 0); } pltscaVarsName.push_back("eta"); ++icomp; @@ -631,7 +631,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo else if (vars[n] == "magvel") { for (int lev = 0; lev <= finest_level; ++lev) { MultiFab magvel(mf[lev], amrex::make_alias, icomp, 1); - ComputeMagVel(lev, m_cur_time, magvel, m_leveldata[lev]->velocity); + ComputeMagVel(lev, cur_time_real(), magvel, m_leveldata[lev]->velocity); } pltscaVarsName.push_back("magvel"); ++icomp; @@ -641,7 +641,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo for (int lev = 0; lev <= finest_level; ++lev) { (m_leveldata[lev]->velocity).FillBoundary(geom[lev].periodicity()); MultiFab vort(mf[lev], amrex::make_alias, icomp, 1); - ComputeVorticity(lev, m_cur_time, vort, m_leveldata[lev]->velocity); + ComputeVorticity(lev, cur_time_real(), vort, m_leveldata[lev]->velocity); } pltscaVarsName.push_back("vort"); ++icomp; @@ -667,7 +667,7 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo &strainrate, &m_leveldata[lev]->velocity, Geom(lev), - m_cur_time, 0); + cur_time_real(), 0); } pltscaVarsName.push_back("strainrate"); ++icomp; @@ -719,5 +719,5 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo // Write the plotfile amrex::WriteMultiLevelPlotfile(plotfilename, finest_level + 1, GetVecOfConstPtrs(mf), - pltscaVarsName, Geom(), m_cur_time, istep, refRatio()); + pltscaVarsName, Geom(), cur_time_real(), istep, refRatio()); }