From 0645f477eb2bc7a265d8f9236565d581fc6d641e Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 15 Jun 2024 08:30:56 -0700 Subject: [PATCH] separate routines into separate files so predictor and corrector aren't so long --- src/Make.package | 4 +- src/incflo.H | 12 +++- src/incflo_apply_corrector.cpp | 100 ++------------------------- src/incflo_apply_predictor.cpp | 120 ++------------------------------- src/incflo_update_density.cpp | 80 ++++++++++++++++++++++ src/incflo_update_tracer.cpp | 67 ++++++++++++++++++ src/setup/incflo_arrays.cpp | 29 ++++---- 7 files changed, 186 insertions(+), 226 deletions(-) create mode 100644 src/incflo_update_density.cpp create mode 100644 src/incflo_update_tracer.cpp diff --git a/src/Make.package b/src/Make.package index e98e6a5a..3c2b87e0 100644 --- a/src/Make.package +++ b/src/Make.package @@ -5,8 +5,10 @@ CEXE_sources += incflo_apply_corrector.cpp CEXE_sources += incflo_compute_dt.cpp CEXE_sources += incflo_compute_forces.cpp CEXE_sources += incflo_explicit_update.cpp -CEXE_sources += incflo_tagging.cpp CEXE_sources += incflo_regrid.cpp +CEXE_sources += incflo_tagging.cpp +CEXE_sources += incflo_update_density.cpp +CEXE_sources += incflo_update_tracer.cpp CEXE_sources += main.cpp ifeq ($(USE_EB), TRUE) diff --git a/src/incflo.H b/src/incflo.H index 16061dba..1006118d 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -21,6 +21,10 @@ #include #include +enum struct StepType { + Predictor, Corrector +}; + class incflo : public amrex::AmrCore { public: @@ -131,8 +135,8 @@ public: // /////////////////////////////////////////////////////////////////////////// - void ApplyPredictor(bool incremental_projection = false); - void ApplyCorrector(); + void ApplyPredictor (bool incremental_projection = false); + void ApplyCorrector (); void compute_convective_term (amrex::Vector const& conv_u, amrex::Vector const& conv_r, amrex::Vector const& conv_t, @@ -158,6 +162,9 @@ public: void tracer_explicit_update(amrex::Vector const& tra_forces); void tracer_explicit_update_corrector(amrex::Vector const& tra_forces); + void update_density (StepType step_type); + void update_tracer (StepType step_type, amrex::Vector& tra_forces); + /////////////////////////////////////////////////////////////////////////// // // derive @@ -589,6 +596,7 @@ private: amrex::MultiFab density; amrex::MultiFab density_eb; amrex::MultiFab density_o; + amrex::MultiFab density_nph; amrex::MultiFab tracer; amrex::MultiFab tracer_eb; amrex::MultiFab tracer_o; diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index da9d0300..aa388010 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -107,7 +107,7 @@ void incflo::ApplyCorrector() // in constructing the advection term // ********************************************************************************************** Vector vel_forces, tra_forces; - Vector vel_eta, tra_eta; + Vector vel_eta; for (int lev = 0; lev <= finest_level; ++lev) { vel_forces.emplace_back(grids[lev], dmap[lev], AMREX_SPACEDIM, nghost_force(), MFInfo(), Factory(lev)); @@ -116,9 +116,6 @@ void incflo::ApplyCorrector() MFInfo(), Factory(lev)); } vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); - if (m_advect_tracer) { - tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev)); - } } // ************************************************************************************* @@ -146,7 +143,6 @@ void incflo::ApplyCorrector() compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 1); - compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); // Here we create divtau of the (n+1,*) state that was computed in the predictor; // we use this laps only if DiffusionType::Explicit @@ -156,11 +152,6 @@ void incflo::ApplyCorrector() get_density_new_const(), GetVecOfConstPtrs(vel_eta)); } - if (m_advect_tracer && m_diff_type == DiffusionType::Explicit) { - compute_laps(get_laps_new(), get_tracer_new_const(), - get_density_new_const(), GetVecOfConstPtrs(tra_eta)); - } - // ************************************************************************************* // Define local variables for lambda to capture. // ************************************************************************************* @@ -169,95 +160,12 @@ void incflo::ApplyCorrector() // ************************************************************************************* // Update density first // ************************************************************************************* - if (m_constant_density) - { - for (int lev = 0; lev <= finest_level; lev++) - MultiFab::Copy(density_nph[lev], m_leveldata[lev]->density_o, 0, 0, 1, 0); - } else { - for (int lev = 0; lev <= finest_level; lev++) - { - auto& ld = *m_leveldata[lev]; -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& bx = mfi.tilebox(); - Array4 const& rho_o = ld.density_o.const_array(mfi); - Array4 const& rho_n = ld.density.array(mfi); - Array4 const& rho_nph = density_nph[lev].array(mfi); - Array4 const& drdt_o = ld.conv_density_o.const_array(mfi); - Array4 const& drdt = ld.conv_density.const_array(mfi); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - const Real rho_old = rho_o(i,j,k); - - Real rho_new = rho_old + l_dt * m_half*(drdt(i,j,k)+drdt_o(i,j,k)); - rho_nph(i,j,k) = m_half * (rho_old + rho_new); - - rho_n (i,j,k) = rho_new; - }); - } // mfi - } // lev - - // Average down solution - for (int lev = finest_level-1; lev >= 0; --lev) { -#ifdef AMREX_USE_EB - amrex::EB_average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density, - 0, 1, refRatio(lev)); -#else - amrex::average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density, - 0, 1, refRatio(lev)); -#endif - } - } // not constant density - - // ************************************************************************************* - // Compute the tracer forcing terms (forcing for (rho s), not for s) - // ************************************************************************************* - if (m_advect_tracer) - compute_tra_forces(GetVecOfPtrs(tra_forces), GetVecOfConstPtrs(density_nph)); + update_density(StepType::Corrector); // ************************************************************************************* - // Update the tracer next (note that dtdt already has rho in it) - // (rho trac)^new = (rho trac)^old + dt * ( - // div(rho trac u) + div (mu grad trac) + rho * f_t + // Update tracer next // ************************************************************************************* - if (m_advect_tracer) - { - tracer_explicit_update_corrector(tra_forces); - } - - // ************************************************************************************* - // Solve diffusion equation for tracer - // ************************************************************************************* - if ( m_advect_tracer ) - { - if (m_diff_type == DiffusionType::Crank_Nicolson || m_diff_type == DiffusionType::Implicit) - { - const int ng_diffusion = 1; - for (int lev = 0; lev <= finest_level; ++lev) - fillphysbc_tracer(lev, new_time, m_leveldata[lev]->tracer, ng_diffusion); - - Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : m_half*m_dt; - diffuse_scalar(get_tracer_new(), get_density_new(), GetVecOfConstPtrs(tra_eta), dt_diff); - } - else - { - // Need to average down tracer since the diffusion solver didn't do it for us. - for (int lev = finest_level-1; lev >= 0; --lev) { -#ifdef AMREX_USE_EB - amrex::EB_average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer, - 0, m_ntrac, refRatio(lev)); -#else - amrex::average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer, - 0, m_ntrac, refRatio(lev)); -#endif - } - } - } // if (m_advect_tracer) - + update_tracer(StepType::Corrector, tra_forces); // ************************************************************************************* // Define the forcing terms to use in the final update (using half-time density) diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index ca161c24..205bc08a 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -106,7 +106,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // Forcing terms Vector vel_forces, tra_forces; - Vector vel_eta, tra_eta; + Vector vel_eta; // ************************************************************************************* // Allocate space for the forcing terms @@ -120,9 +120,6 @@ void incflo::ApplyPredictor (bool incremental_projection) MFInfo(), Factory(lev)); } vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); - if (m_advect_tracer) { - tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev)); - } } // ************************************************************************************* @@ -135,27 +132,16 @@ void incflo::ApplyPredictor (bool incremental_projection) compute_viscosity(GetVecOfPtrs(vel_eta), get_density_old(), get_velocity_old(), m_cur_time, 1); - compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); // ************************************************************************************* // Compute explicit viscous term // Note that for !advect_momentum, this actually computes divtau / rho // ************************************************************************************* - if (need_divtau() || use_tensor_correction ) - { + if (need_divtau() || use_tensor_correction ) { compute_divtau(get_divtau_old(),get_velocity_old_const(), get_density_old_const(),GetVecOfConstPtrs(vel_eta)); } - // ************************************************************************************* - // Compute explicit diffusive terms - // ************************************************************************************* - if (m_advect_tracer && need_divtau()) - { - compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(), - GetVecOfConstPtrs(tra_eta)); - } - // ********************************************************************************************** // Compute the MAC-projected velocities at all levels // ************************************************************************************* @@ -188,107 +174,12 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* // Update density first // ************************************************************************************* - if (m_constant_density) - { - for (int lev = 0; lev <= finest_level; lev++) - MultiFab::Copy(density_nph[lev], m_leveldata[lev]->density_o, 0, 0, 1, 1); - } else { - for (int lev = 0; lev <= finest_level; lev++) - { - auto& ld = *m_leveldata[lev]; -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& bx = mfi.tilebox(); - Array4 const& rho_o = ld.density_o.const_array(mfi); - Array4 const& rho_new = ld.density.array(mfi); - Array4 const& drdt = ld.conv_density_o.const_array(mfi); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - rho_new(i,j,k) = rho_o(i,j,k) + l_dt * drdt(i,j,k); - }); - } // mfi - - // Fill ghost cells of the new density field so that we can define density_nph - // on the valid region grown by 1 - int ng = 1; - fillpatch_density(lev, m_t_new[lev], ld.density, ng); - - for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& gbx = mfi.growntilebox(1); - Array4 const& rho_old = ld.density_o.const_array(mfi); - Array4 const& rho_new = ld.density.const_array(mfi); - Array4 const& rho_nph = density_nph[lev].array(mfi); - - amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - rho_nph(i,j,k) = m_half * (rho_old(i,j,k) + rho_new(i,j,k)); - }); - } // mfi - - } // lev - - // Average down solution - for (int lev = finest_level-1; lev >= 0; --lev) { -#ifdef AMREX_USE_EB - amrex::EB_average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density, - 0, 1, refRatio(lev)); -#else - amrex::average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density, - 0, 1, refRatio(lev)); -#endif - } - } // not constant density - - - // ************************************************************************************* - // Compute (or if Godunov, re-compute) the tracer forcing terms - // (forcing for (rho s) if conservative) - // ************************************************************************************* - if (m_advect_tracer) - compute_tra_forces(GetVecOfPtrs(tra_forces), GetVecOfConstPtrs(density_nph)); + update_density(StepType::Predictor); // ************************************************************************************* - // Update the tracer next + // Update tracer next // ************************************************************************************* - if (m_advect_tracer) - { - tracer_explicit_update(tra_forces); - } - - // ************************************************************************************* - // Solve diffusion equation for tracer - // ************************************************************************************* - if ( m_advect_tracer ) - { - if (m_diff_type == DiffusionType::Crank_Nicolson || m_diff_type == DiffusionType::Implicit) - { - const int ng_diffusion = 1; - for (int lev = 0; lev <= finest_level; ++lev) - fillphysbc_tracer(lev, new_time, m_leveldata[lev]->tracer, ng_diffusion); - - Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : m_half*m_dt; - diffuse_scalar(get_tracer_new(), get_density_new(), GetVecOfConstPtrs(tra_eta), dt_diff); - } - else - { - // Need to average down tracer since the diffusion solver didn't do it for us. - for (int lev = finest_level-1; lev >= 0; --lev) { -#ifdef AMREX_USE_EB - amrex::EB_average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer, - 0, m_ntrac, refRatio(lev)); -#else - amrex::average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer, - 0, m_ntrac, refRatio(lev)); -#endif - } - } - } // if (m_advect_tracer) - + update_tracer(StepType::Predictor, tra_forces); // ************************************************************************************* // Define (or if advection_type != "MOL", re-define) the forcing terms, without the viscous terms @@ -298,7 +189,6 @@ void incflo::ApplyPredictor (bool incremental_projection) GetVecOfConstPtrs(density_nph), get_tracer_old_const(), get_tracer_new_const()); - // ************************************************************************************* // Update the velocity // ************************************************************************************* diff --git a/src/incflo_update_density.cpp b/src/incflo_update_density.cpp new file mode 100644 index 00000000..63787640 --- /dev/null +++ b/src/incflo_update_density.cpp @@ -0,0 +1,80 @@ +#include + +using namespace amrex; + +void incflo::update_density (StepType step_type) +{ + BL_PROFILE("incflo::update_density"); + + int ng = (step_type == StepType::Corrector) ? 0 : 1; + + Real l_dt = m_dt; + + if (m_constant_density) + { + for (int lev = 0; lev <= finest_level; lev++) + MultiFab::Copy(m_leveldata[lev]->density_nph, m_leveldata[lev]->density_o, 0, 0, 1, ng); + + } else { + for (int lev = 0; lev <= finest_level; lev++) + { + auto& ld = *m_leveldata[lev]; +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + Array4 const& rho_old = ld.density_o.const_array(mfi); + Array4 const& rho_new = ld.density.array(mfi); + Array4 const& rho_nph = ld.density_nph.array(mfi); + Array4 const& drdt_o = ld.conv_density_o.const_array(mfi); + + if (step_type == StepType::Predictor) { + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + rho_new(i,j,k) = rho_old(i,j,k) + l_dt * drdt_o(i,j,k); + }); + + } else if (step_type == StepType::Corrector) { + Array4 const& drdt_n = ld.conv_density.const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + rho_new(i,j,k) = rho_old(i,j,k) + l_dt * Real(0.5) * (drdt_o(i,j,k) + drdt_n(i,j,k)); + rho_nph(i,j,k) = Real(0.5) * (rho_old(i,j,k) + rho_new(i,j,k)); + }); + } + } // mfi + + if (step_type == StepType::Predictor) { + // Fill ghost cells of the new density field so that we can define density_nph + // on the valid region grown by 1 + fillpatch_density(lev, m_t_new[lev], ld.density, ng); + + for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& gbx = mfi.growntilebox(1); + Array4 const& rho_new = ld.density.const_array(mfi); + Array4 const& rho_old = ld.density_o.const_array(mfi); + Array4 const& rho_nph = ld.density_nph.array(mfi); + + ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + rho_nph(i,j,k) = Real(0.5) * (rho_old(i,j,k) + rho_new(i,j,k)); + }); + } // mfi + } // Predictor + } // lev + + // Average down solution + for (int lev = finest_level-1; lev >= 0; --lev) { +#ifdef AMREX_USE_EB + amrex::EB_average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density, + 0, 1, refRatio(lev)); +#else + amrex::average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density, + 0, 1, refRatio(lev)); +#endif + } + } // not constant density +} diff --git a/src/incflo_update_tracer.cpp b/src/incflo_update_tracer.cpp new file mode 100644 index 00000000..a3a3e8b0 --- /dev/null +++ b/src/incflo_update_tracer.cpp @@ -0,0 +1,67 @@ +#include + +using namespace amrex; + +void incflo::update_tracer (StepType step_type, Vector& tra_forces) +{ + BL_PROFILE("incflo::update_tracer"); + + Vector tra_eta; + + Real new_time = m_cur_time + m_dt; + + if (m_advect_tracer) + { + // ************************************************************************************* + // Compute the tracer forcing terms ( forcing for (rho s) if conservative ) + // ************************************************************************************* + for (int lev = 0; lev <= finest_level; ++lev) { + tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev)); + } + + compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + + // ************************************************************************************* + // Compute explicit diffusive terms + // ************************************************************************************* + if (step_type == StepType::Predictor && need_divtau()) { + compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(), + GetVecOfConstPtrs(tra_eta)); + } else if (step_type == StepType::Corrector && m_diff_type == DiffusionType::Explicit) { + compute_laps(get_laps_new(), get_tracer_new_const(), get_density_new_const(), + GetVecOfConstPtrs(tra_eta)); + } + + if (step_type == StepType::Predictor) { + tracer_explicit_update(tra_forces); + } else if (step_type == StepType::Corrector) { + tracer_explicit_update_corrector(tra_forces); + } + + // ************************************************************************************* + // Solve diffusion equation for tracer + // ************************************************************************************* + if (m_diff_type == DiffusionType::Crank_Nicolson || m_diff_type == DiffusionType::Implicit) + { + const int ng_diffusion = 1; + for (int lev = 0; lev <= finest_level; ++lev) + 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; + diffuse_scalar(get_tracer_new(), get_density_new(), GetVecOfConstPtrs(tra_eta), dt_diff); + } + else + { + // Need to average down tracer since the diffusion solver didn't do it for us. + for (int lev = finest_level-1; lev >= 0; --lev) { +#ifdef AMREX_USE_EB + amrex::EB_average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer, + 0, m_ntrac, refRatio(lev)); +#else + amrex::average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer, + 0, m_ntrac, refRatio(lev)); +#endif + } + } + } // advect tracer +} diff --git a/src/setup/incflo_arrays.cpp b/src/setup/incflo_arrays.cpp index be3e4cf7..d903c7e8 100644 --- a/src/setup/incflo_arrays.cpp +++ b/src/setup/incflo_arrays.cpp @@ -8,19 +8,24 @@ incflo::LevelData::LevelData (amrex::BoxArray const& ba, int ntrac, int ng_state, const std::string& advection_type, bool implicit_diffusion, bool use_tensor_correction, bool advect_tracer) - : velocity (ba, dm, AMREX_SPACEDIM, ng_state, MFInfo(), fact), - velocity_o(ba, dm, AMREX_SPACEDIM, ng_state, MFInfo(), fact), - velocity_eb(ba, dm, AMREX_SPACEDIM, ng_state, MFInfo(), fact), - density (ba, dm, 1 , ng_state, MFInfo(), fact), - density_eb(ba, dm, 1 , ng_state, MFInfo(), fact), - density_o (ba, dm, 1 , ng_state, MFInfo(), fact), - tracer (ba, dm, ntrac , ng_state, MFInfo(), fact), - tracer_eb (ba, dm, ntrac , ng_state, MFInfo(), fact), - tracer_o (ba, dm, ntrac , ng_state, MFInfo(), fact), - mac_phi (ba, dm, 1 , 1 , MFInfo(), fact), - p_nd (amrex::convert(ba,IntVect::TheNodeVector()), + : velocity (ba, dm, AMREX_SPACEDIM, ng_state, MFInfo(), fact), + velocity_o (ba, dm, AMREX_SPACEDIM, ng_state, MFInfo(), fact), + velocity_eb (ba, dm, AMREX_SPACEDIM, ng_state, MFInfo(), fact), + + density (ba, dm, 1 , ng_state, MFInfo(), fact), + density_eb (ba, dm, 1 , ng_state, MFInfo(), fact), + density_o (ba, dm, 1 , ng_state, MFInfo(), fact), + density_nph (ba, dm, 1 , ng_state, MFInfo(), fact), + + tracer (ba, dm, ntrac , ng_state, MFInfo(), fact), + tracer_eb (ba, dm, ntrac , ng_state, MFInfo(), fact), + tracer_o (ba, dm, ntrac , ng_state, MFInfo(), fact), + + mac_phi (ba, dm, 1 , 1 , MFInfo(), fact), + p_nd (amrex::convert(ba,IntVect::TheNodeVector()), dm, 1 , 0 , MFInfo(), fact), - gp (ba, dm, AMREX_SPACEDIM, 0 , MFInfo(), fact), + gp (ba, dm, AMREX_SPACEDIM, 0 , MFInfo(), fact), + conv_velocity_o(ba, dm, AMREX_SPACEDIM, 0, MFInfo(), fact), conv_density_o (ba, dm, 1 , 0, MFInfo(), fact), conv_tracer_o (ba, dm, ntrac , 0, MFInfo(), fact)