From ff6374bea8948070219bf90fe87249bb1fd002cf Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 15 Jun 2024 10:23:23 -0700 Subject: [PATCH] pull velocity update into separate file --- src/CMakeLists.txt | 2 + src/Make.package | 2 + src/incflo.H | 8 +- src/incflo.cpp | 342 ------------------------------ src/incflo_apply_corrector.cpp | 217 +------------------- src/incflo_apply_predictor.cpp | 166 +-------------- src/incflo_update_velocity.cpp | 294 ++++++++++++++++++++++++++ src/incflo_utils.cpp | 365 +++++++++++++++++++++++++++++++++ 8 files changed, 676 insertions(+), 720 deletions(-) create mode 100644 src/incflo_update_velocity.cpp create mode 100644 src/incflo_utils.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3e8d6f07..f6602817 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -15,6 +15,8 @@ target_sources(incflo incflo_tagging.cpp incflo_update_density.cpp incflo_update_tracer.cpp + incflo_update_velocity.cpp + incflo_utils.cpp main.cpp ) diff --git a/src/Make.package b/src/Make.package index 3c2b87e0..135785eb 100644 --- a/src/Make.package +++ b/src/Make.package @@ -9,6 +9,8 @@ CEXE_sources += incflo_regrid.cpp CEXE_sources += incflo_tagging.cpp CEXE_sources += incflo_update_density.cpp CEXE_sources += incflo_update_tracer.cpp +CEXE_sources += incflo_update_velocity.cpp +CEXE_sources += incflo_utils.cpp CEXE_sources += main.cpp ifeq ($(USE_EB), TRUE) diff --git a/src/incflo.H b/src/incflo.H index 1006118d..b9accecc 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -162,8 +162,10 @@ 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); + void update_density (StepType step_type); + void update_tracer (StepType step_type, amrex::Vector& tra_forces); + void update_velocity (StepType step_type, amrex::Vector& vel_eta, + amrex::Vector& vel_forces); /////////////////////////////////////////////////////////////////////////// // @@ -753,6 +755,7 @@ private: amrex::Vector get_tracer_eb () noexcept; amrex::Vector get_density_old () noexcept; amrex::Vector get_density_new () noexcept; + amrex::Vector get_density_nph () noexcept; amrex::Vector get_tracer_old () noexcept; amrex::Vector get_tracer_new () noexcept; amrex::Vector get_mac_phi () noexcept; @@ -776,6 +779,7 @@ private: [[nodiscard]] amrex::Vector get_tracer_eb () const noexcept; [[nodiscard]] amrex::Vector get_density_old_const () const noexcept; [[nodiscard]] amrex::Vector get_density_new_const () const noexcept; + [[nodiscard]] amrex::Vector get_density_nph_const () const noexcept; [[nodiscard]] amrex::Vector get_tracer_old_const () const noexcept; [[nodiscard]] amrex::Vector get_tracer_new_const () const noexcept; [[nodiscard]] amrex::Vector get_vel_forces_const () const noexcept; diff --git a/src/incflo.cpp b/src/incflo.cpp index fd6ba731..8f4c5906 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -295,345 +295,3 @@ incflo::writeNow() return write_now; } - -Vector incflo::get_velocity_eb () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->velocity_eb)); - } - return r; -} - -Vector incflo::get_density_eb () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->density_eb)); - } - return r; -} - -Vector incflo::get_tracer_eb () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->tracer_eb)); - } - return r; -} - -Vector incflo::get_velocity_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->velocity_o)); - } - return r; -} - -Vector incflo::get_velocity_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->velocity)); - } - return r; -} - -Vector incflo::get_density_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->density_o)); - } - return r; -} - -Vector incflo::get_density_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->density)); - } - return r; -} - -Vector incflo::get_tracer_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->tracer_o)); - } - return r; -} - -Vector incflo::get_tracer_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->tracer)); - } - return r; -} - -Vector incflo::get_mac_phi() noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->mac_phi)); - } - return r; -} - -Vector incflo::get_conv_velocity_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->conv_velocity_o)); - } - return r; -} - -Vector incflo::get_conv_velocity_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->conv_velocity)); - } - return r; -} - -Vector incflo::get_conv_density_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->conv_density_o)); - } - return r; -} - -Vector incflo::get_conv_density_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->conv_density)); - } - return r; -} - -Vector incflo::get_conv_tracer_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->conv_tracer_o)); - } - return r; -} - -Vector incflo::get_conv_tracer_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->conv_tracer)); - } - return r; -} - -Vector incflo::get_divtau_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->divtau_o)); - } - return r; -} - -Vector incflo::get_divtau_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->divtau)); - } - return r; -} - -Vector incflo::get_laps_old () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->laps_o)); - } - return r; -} - -Vector incflo::get_laps_new () noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->laps)); - } - return r; -} - -Vector incflo::get_velocity_old_const () const noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->velocity_o)); - } - return r; -} - -Vector incflo::get_velocity_new_const () const noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->velocity)); - } - return r; -} - -Vector incflo::get_density_old_const () const noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->density_o)); - } - return r; -} - -Vector incflo::get_density_new_const () const noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->density)); - } - return r; -} - -Vector incflo::get_tracer_old_const () const noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->tracer_o)); - } - return r; -} - -Vector incflo::get_tracer_new_const () const noexcept -{ - Vector r; - r.reserve(finest_level+1); - for (int lev = 0; lev <= finest_level; ++lev) { - r.push_back(&(m_leveldata[lev]->tracer)); - } - return r; -} - -void incflo::copy_from_new_to_old_velocity (IntVect const& ng) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - copy_from_new_to_old_velocity(lev, ng); - } -} - -void incflo::copy_from_new_to_old_velocity (int lev, IntVect const& ng) -{ - MultiFab::Copy(m_leveldata[lev]->velocity_o, - m_leveldata[lev]->velocity, 0, 0, AMREX_SPACEDIM, ng); -} - -void incflo::copy_from_old_to_new_velocity (IntVect const& ng) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - copy_from_old_to_new_velocity(lev, ng); - } -} - -void incflo::copy_from_old_to_new_velocity (int lev, IntVect const& ng) -{ - MultiFab::Copy(m_leveldata[lev]->velocity, - m_leveldata[lev]->velocity_o, 0, 0, AMREX_SPACEDIM, ng); -} - -void incflo::copy_from_new_to_old_density (IntVect const& ng) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - copy_from_new_to_old_density(lev, ng); - } -} - -void incflo::copy_from_new_to_old_density (int lev, IntVect const& ng) -{ - MultiFab::Copy(m_leveldata[lev]->density_o, - m_leveldata[lev]->density, 0, 0, 1, ng); -} - -void incflo::copy_from_old_to_new_density (IntVect const& ng) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - copy_from_old_to_new_density(lev, ng); - } -} - -void incflo::copy_from_old_to_new_density (int lev, IntVect const& ng) -{ - MultiFab::Copy(m_leveldata[lev]->density, - m_leveldata[lev]->density_o, 0, 0, 1, ng); -} - -void incflo::copy_from_new_to_old_tracer (IntVect const& ng) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - copy_from_new_to_old_tracer(lev, ng); - } -} - -void incflo::copy_from_new_to_old_tracer (int lev, IntVect const& ng) -{ - if (m_ntrac > 0) { - MultiFab::Copy(m_leveldata[lev]->tracer_o, - m_leveldata[lev]->tracer, 0, 0, m_ntrac, ng); - } -} - -void incflo::copy_from_old_to_new_tracer (IntVect const& ng) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - copy_from_old_to_new_tracer(lev, ng); - } -} - -void incflo::copy_from_old_to_new_tracer (int lev, IntVect const& ng) -{ - if (m_ntrac > 0) { - MultiFab::Copy(m_leveldata[lev]->tracer, - m_leveldata[lev]->tracer_o, 0, 0, m_ntrac, ng); - } -} diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index aa388010..27df84ef 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -69,8 +69,6 @@ void incflo::ApplyCorrector() { BL_PROFILE("incflo::ApplyCorrector"); - constexpr Real m_half = Real(0.5); - // We use the new time value for things computed on the "*" state Real new_time = m_cur_time + m_dt; @@ -94,14 +92,6 @@ void incflo::ApplyCorrector() } } - // ************************************************************************************* - // Allocate space for half-time density - // ************************************************************************************* - Vector density_nph; - for (int lev = 0; lev <= finest_level; ++lev) { - density_nph.emplace_back(grids[lev], dmap[lev], 1, 0, MFInfo(), Factory(lev)); - } - // ********************************************************************************************** // We only reach the corrector if advection_type == MOL which means we don't use the forces // in constructing the advection term @@ -152,11 +142,6 @@ void incflo::ApplyCorrector() get_density_new_const(), GetVecOfConstPtrs(vel_eta)); } - // ************************************************************************************* - // Define local variables for lambda to capture. - // ************************************************************************************* - Real l_dt = m_dt; - // ************************************************************************************* // Update density first // ************************************************************************************* @@ -167,214 +152,16 @@ void incflo::ApplyCorrector() // ************************************************************************************* update_tracer(StepType::Corrector, tra_forces); - // ************************************************************************************* - // Define the forcing terms to use in the final update (using half-time density) - // ************************************************************************************* - compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_new_const(), - GetVecOfConstPtrs(density_nph), - get_tracer_old_const(), get_tracer_new_const()); - // ************************************************************************************* // Update velocity // ************************************************************************************* - 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& vel = ld.velocity.array(mfi); - Array4 const& vel_o = ld.velocity_o.const_array(mfi); - Array4 const& dvdt = ld.conv_velocity.const_array(mfi); - Array4 const& dvdt_o = ld.conv_velocity_o.const_array(mfi); - Array4 const& vel_f = vel_forces[lev].const_array(mfi); - - 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); - - if (m_diff_type == DiffusionType::Explicit) - { - Array4 const& divtau_o = ld.divtau_o.const_array(mfi); - Array4 const& divtau = ld.divtau.const_array(mfi); - - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( - m_half*( dvdt_o(i,j,k,0)+ dvdt(i,j,k,0)) - + m_half*(divtau_o(i,j,k,0)+divtau(i,j,k,0)) - + rho_nph(i,j,k) * vel_f(i,j,k,0) );, - vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( - m_half*( dvdt_o(i,j,k,1)+ dvdt(i,j,k,1)) - + m_half*(divtau_o(i,j,k,1)+divtau(i,j,k,1)) - + rho_nph(i,j,k) * vel_f(i,j,k,1) );, - vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( - m_half*( dvdt_o(i,j,k,2)+ dvdt(i,j,k,2)) - + m_half*(divtau_o(i,j,k,2)+divtau(i,j,k,2)) - + rho_nph(i,j,k) * vel_f(i,j,k,2) );); - - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( - m_half*( dvdt_o(i,j,k,0)+ dvdt(i,j,k,0)) - + m_half*(divtau_o(i,j,k,0)+divtau(i,j,k,0)) - + vel_f(i,j,k,0) );, - vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( - m_half*( dvdt_o(i,j,k,1)+ dvdt(i,j,k,1)) - + m_half*(divtau_o(i,j,k,1)+divtau(i,j,k,1)) - + vel_f(i,j,k,1) );, - vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( - m_half*( dvdt_o(i,j,k,2)+ dvdt(i,j,k,2)) - + m_half*(divtau_o(i,j,k,2)+divtau(i,j,k,2)) - + vel_f(i,j,k,2) );); - }); - } - } - else if (m_diff_type == DiffusionType::Crank_Nicolson) - { - Array4 const& divtau_o = ld.divtau_o.const_array(mfi); - - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( - m_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) - +m_half*(divtau_o(i,j,k,0) ) + rho_nph(i,j,k)*vel_f(i,j,k,0) );, - vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( - m_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) - +m_half*(divtau_o(i,j,k,1) ) + rho_nph(i,j,k)*vel_f(i,j,k,1) );, - vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( - m_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) - +m_half*(divtau_o(i,j,k,2) ) + rho_nph(i,j,k)*vel_f(i,j,k,2) );); - - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( - m_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) - + m_half* divtau_o(i,j,k,0) + vel_f(i,j,k,0) );, - vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( - m_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) - + m_half* divtau_o(i,j,k,1) + vel_f(i,j,k,1) );, - vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( - m_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) - + m_half* divtau_o(i,j,k,2) + vel_f(i,j,k,2) );); - }); - } - } - else if (m_diff_type == DiffusionType::Implicit) - { - if (use_tensor_correction) - { - Array4 const& divtau = ld.divtau.const_array(mfi); - if (m_advect_momentum) - { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( - m_half*(dvdt_o(i,j,k,0) + dvdt(i,j,k,0)) - + rho_nph(i,j,k) * vel_f(i,j,k,0) + divtau(i,j,k,0));, -// - vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( - m_half*(dvdt_o(i,j,k,1) + dvdt(i,j,k,1)) - + rho_nph(i,j,k) * vel_f(i,j,k,1) + divtau(i,j,k,1));, -// - vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( - m_half*(dvdt_o(i,j,k,2) + dvdt(i,j,k,2)) - + rho_nph(i,j,k) * vel_f(i,j,k,2) + divtau(i,j,k,2));); - - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( - m_half*( dvdt_o(i,j,k,0) + dvdt(i,j,k,0)) - + vel_f(i,j,k,0) + divtau(i,j,k,0));, -// - vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( - m_half*(dvdt_o(i,j,k,1) + dvdt(i,j,k,1)) - + vel_f(i,j,k,1) + divtau(i,j,k,1) );, -// - vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( - m_half*( dvdt_o(i,j,k,2) + dvdt(i,j,k,2)) - + vel_f(i,j,k,2) + divtau(i,j,k,2) );); - }); - } - } else { - if (m_advect_momentum) - { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * ( - m_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) - + rho_nph(i,j,k) * vel_f(i,j,k,0) );, - vel(i,j,k,1) = rho_old(i,j,k) * vel_o(i,j,k,1) + l_dt * ( - m_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) - + rho_nph(i,j,k) * vel_f(i,j,k,1) );, - vel(i,j,k,2) = rho_old(i,j,k) * vel_o(i,j,k,2) + l_dt * ( - m_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) - + rho_nph(i,j,k) * vel_f(i,j,k,2) );); - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * ( - m_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + vel_f(i,j,k,0) );, - vel(i,j,k,1) = vel_o(i,j,k,1) + l_dt * ( - m_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) + vel_f(i,j,k,1) );, - vel(i,j,k,2) = vel_o(i,j,k,2) + l_dt * ( - m_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + vel_f(i,j,k,2) );); - }); - } - } - } - } - } - - // ********************************************************************************************** - // - // Solve diffusion equation for u* at t^{n+1} but using eta at predicted new time - // - // ********************************************************************************************** - - 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_velocity(lev, new_time, m_leveldata[lev]->velocity, ng_diffusion); - fillphysbc_density (lev, new_time, m_leveldata[lev]->density , ng_diffusion); - } - - Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : m_half*m_dt; - diffuse_velocity(get_velocity_new(), get_density_new(), GetVecOfConstPtrs(vel_eta), dt_diff); - } + update_velocity(StepType::Corrector, vel_eta, vel_forces); // ********************************************************************************************** // // Project velocity field, update pressure bool incremental_projection = false; - ApplyProjection(GetVecOfConstPtrs(density_nph), new_time,m_dt,incremental_projection); + ApplyProjection(get_density_nph_const(), new_time,m_dt,incremental_projection); #ifdef INCFLO_USE_PARTICLES // ************************************************************************************** diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index 205bc08a..ebbcc734 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -70,8 +70,6 @@ void incflo::ApplyPredictor (bool incremental_projection) { BL_PROFILE("incflo::ApplyPredictor"); - constexpr Real m_half = Real(0.5); - // We use the new time value for things computed on the "*" state Real new_time = m_cur_time + m_dt; @@ -96,13 +94,6 @@ void incflo::ApplyPredictor (bool incremental_projection) } } - // ************************************************************************************* - // Allocate space for half-time density - // ************************************************************************************* - Vector density_nph; - for (int lev = 0; lev <= finest_level; ++lev) - density_nph.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); - // Forcing terms Vector vel_forces, tra_forces; @@ -167,173 +158,26 @@ void incflo::ApplyPredictor (bool incremental_projection) m_cur_time); // ************************************************************************************* - // Define local variables for lambda to capture. - // ************************************************************************************* - Real l_dt = m_dt; - - // ************************************************************************************* - // Update density first + // Update density // ************************************************************************************* update_density(StepType::Predictor); // ************************************************************************************* - // Update tracer next + // Update tracer // ************************************************************************************* update_tracer(StepType::Predictor, tra_forces); // ************************************************************************************* - // Define (or if advection_type != "MOL", re-define) the forcing terms, without the viscous terms - // and using the half-time density - // ************************************************************************************* - compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_old_const(), - GetVecOfConstPtrs(density_nph), - get_tracer_old_const(), get_tracer_new_const()); - - // ************************************************************************************* - // Update the velocity + // Update velocity // ************************************************************************************* - 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& vel = ld.velocity.array(mfi); - Array4 const& dvdt = ld.conv_velocity_o.const_array(mfi); - Array4 const& vel_f = vel_forces[lev].const_array(mfi); - 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); - - if (m_diff_type == DiffusionType::Implicit) { - - if (use_tensor_correction) - { - Array4 const& divtau_o = ld.divtau_o.const_array(mfi); - // Here divtau_o is the difference of tensor and scalar divtau_o! - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); - }); - } - } else { - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2));); - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2));); - }); - } - } - } - else if (m_diff_type == DiffusionType::Crank_Nicolson) - { - - Array4 const& divtau_o = ld.divtau_o.const_array(mfi); - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+m_half*divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+m_half*divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+m_half*divtau_o(i,j,k,2));); - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+m_half*divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+m_half*divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+m_half*divtau_o(i,j,k,2));); - }); - } - } - else if (m_diff_type == DiffusionType::Explicit) - { - Array4 const& divtau_o = ld.divtau_o.const_array(mfi); - if (m_advect_momentum) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, - vel(i,j,k,1) *= rho_old(i,j,k);, - vel(i,j,k,2) *= rho_old(i,j,k);); - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); - AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, - vel(i,j,k,1) /= rho_new(i,j,k);, - vel(i,j,k,2) /= rho_new(i,j,k);); - }); - } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, - vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, - vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); - }); - } - } - } // mfi - } // lev - - // ************************************************************************************* - // Solve diffusion equation for u* but using eta_old at old time - // ************************************************************************************* - 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_velocity(lev, new_time, m_leveldata[lev]->velocity, ng_diffusion); - fillphysbc_density (lev, new_time, m_leveldata[lev]->density , ng_diffusion); - } - - Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : m_half*m_dt; - diffuse_velocity(get_velocity_new(), get_density_new(), GetVecOfConstPtrs(vel_eta), dt_diff); - } + update_velocity(StepType::Predictor, vel_eta, vel_forces); // ********************************************************************************************** // // Project velocity field, update pressure // // ********************************************************************************************** - ApplyProjection(GetVecOfConstPtrs(density_nph),new_time,m_dt,incremental_projection); + ApplyProjection(get_density_nph_const(),new_time,m_dt,incremental_projection); #ifdef INCFLO_USE_PARTICLES // ************************************************************************************** diff --git a/src/incflo_update_velocity.cpp b/src/incflo_update_velocity.cpp new file mode 100644 index 00000000..48fbd183 --- /dev/null +++ b/src/incflo_update_velocity.cpp @@ -0,0 +1,294 @@ +#include + +using namespace amrex; + +void incflo::update_velocity (StepType step_type, Vector& vel_eta, Vector& vel_forces) +{ + BL_PROFILE("incflo::update_velocity"); + + Real new_time = m_cur_time + m_dt; + + Real l_dt = m_dt; + Real l_half = Real(0.5); + + if (step_type == StepType::Predictor) { + + // ************************************************************************************* + // Define (or if advection_type != "MOL", re-define) the forcing terms, without the viscous terms + // and using the half-time density + // ************************************************************************************* + compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_old_const(), + get_density_nph_const(), get_tracer_old_const(), get_tracer_new_const()); + + // ************************************************************************************* + // Update the velocity + // ************************************************************************************* + 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& vel = ld.velocity.array(mfi); + Array4 const& dvdt = ld.conv_velocity_o.const_array(mfi); + Array4 const& vel_f = vel_forces[lev].const_array(mfi); + Array4 const& rho_old = ld.density_o.const_array(mfi); + Array4 const& rho_new = ld.density.const_array(mfi); + Array4 const& rho_nph = ld.density_nph.const_array(mfi); + + if (m_diff_type == DiffusionType::Implicit) { + + if (use_tensor_correction) + { + Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + // Here divtau_o is the difference of tensor and scalar divtau_o! + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + }); + } + } else { + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2));); + }); + } + } + } + else if (m_diff_type == DiffusionType::Crank_Nicolson) + { + + Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+l_half*divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+l_half*divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+l_half*divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+l_half*divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+l_half*divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+l_half*divtau_o(i,j,k,2));); + }); + } + } + else if (m_diff_type == DiffusionType::Explicit) + { + Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + }); + } + } + } // mfi + } // lev + + } else if (step_type == StepType::Corrector) { + + // ************************************************************************************* + // Define the forcing terms to use in the final update (using half-time density) + // ************************************************************************************* + compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_new_const(), + get_density_nph_const(), get_tracer_old_const(), get_tracer_new_const()); + + 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& vel = ld.velocity.array(mfi); + Array4 const& dvdt = ld.conv_velocity_o.const_array(mfi); + Array4 const& vel_f = vel_forces[lev].const_array(mfi); + Array4 const& rho_old = ld.density_o.const_array(mfi); + Array4 const& rho_new = ld.density.const_array(mfi); + Array4 const& rho_nph = ld.density_nph.const_array(mfi); + + if (m_diff_type == DiffusionType::Implicit) { + + if (use_tensor_correction) + { + Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + // Here divtau_o is the difference of tensor and scalar divtau_o! + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + }); + } + } else { + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2));); + }); + } + } + } + else if (m_diff_type == DiffusionType::Crank_Nicolson) + { + + Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+l_half*divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+l_half*divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+l_half*divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+l_half*divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+l_half*divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+l_half*divtau_o(i,j,k,2));); + }); + } + } + else if (m_diff_type == DiffusionType::Explicit) + { + Array4 const& divtau_o = ld.divtau_o.const_array(mfi); + if (m_advect_momentum) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);, + vel(i,j,k,1) *= rho_old(i,j,k);, + vel(i,j,k,2) *= rho_old(i,j,k);); + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+rho_nph(i,j,k)*vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+rho_nph(i,j,k)*vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+rho_nph(i,j,k)*vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + AMREX_D_TERM(vel(i,j,k,0) /= rho_new(i,j,k);, + vel(i,j,k,1) /= rho_new(i,j,k);, + vel(i,j,k,2) /= rho_new(i,j,k);); + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));, + vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));, + vel(i,j,k,2) += l_dt*(dvdt(i,j,k,2)+vel_f(i,j,k,2)+divtau_o(i,j,k,2));); + }); + } + } + } // mfi + } // lev + + } // Corrector + + // ************************************************************************************* + // Solve diffusion equation for u* but using eta_old at old time + // ************************************************************************************* + 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_velocity(lev, new_time, m_leveldata[lev]->velocity, ng_diffusion); + 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; + diffuse_velocity(get_velocity_new(), get_density_new(), GetVecOfConstPtrs(vel_eta), dt_diff); + } +} + diff --git a/src/incflo_utils.cpp b/src/incflo_utils.cpp new file mode 100644 index 00000000..3e155e1a --- /dev/null +++ b/src/incflo_utils.cpp @@ -0,0 +1,365 @@ +#include + +using namespace amrex; + +Vector incflo::get_velocity_eb () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->velocity_eb)); + } + return r; +} + +Vector incflo::get_density_eb () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->density_eb)); + } + return r; +} + +Vector incflo::get_tracer_eb () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->tracer_eb)); + } + return r; +} + +Vector incflo::get_velocity_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->velocity_o)); + } + return r; +} + +Vector incflo::get_velocity_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->velocity)); + } + return r; +} + +Vector incflo::get_density_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->density_o)); + } + return r; +} + +Vector incflo::get_density_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->density)); + } + return r; +} + +Vector incflo::get_density_nph () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->density_nph)); + } + return r; +} + +Vector incflo::get_tracer_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->tracer_o)); + } + return r; +} + +Vector incflo::get_tracer_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->tracer)); + } + return r; +} + +Vector incflo::get_mac_phi() noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->mac_phi)); + } + return r; +} + +Vector incflo::get_conv_velocity_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_velocity_o)); + } + return r; +} + +Vector incflo::get_conv_velocity_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_velocity)); + } + return r; +} + +Vector incflo::get_conv_density_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_density_o)); + } + return r; +} + +Vector incflo::get_conv_density_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_density)); + } + return r; +} + +Vector incflo::get_conv_tracer_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_tracer_o)); + } + return r; +} + +Vector incflo::get_conv_tracer_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->conv_tracer)); + } + return r; +} + +Vector incflo::get_divtau_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->divtau_o)); + } + return r; +} + +Vector incflo::get_divtau_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->divtau)); + } + return r; +} + +Vector incflo::get_laps_old () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->laps_o)); + } + return r; +} + +Vector incflo::get_laps_new () noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->laps)); + } + return r; +} + +Vector incflo::get_velocity_old_const () const noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->velocity_o)); + } + return r; +} + +Vector incflo::get_velocity_new_const () const noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->velocity)); + } + return r; +} + +Vector incflo::get_density_old_const () const noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->density_o)); + } + return r; +} + +Vector incflo::get_density_new_const () const noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->density)); + } + return r; +} + +Vector incflo::get_density_nph_const () const noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->density_nph)); + } + return r; +} + +Vector incflo::get_tracer_old_const () const noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->tracer_o)); + } + return r; +} + +Vector incflo::get_tracer_new_const () const noexcept +{ + Vector r; + r.reserve(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + r.push_back(&(m_leveldata[lev]->tracer)); + } + return r; +} + +void incflo::copy_from_new_to_old_velocity (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + copy_from_new_to_old_velocity(lev, ng); + } +} + +void incflo::copy_from_new_to_old_velocity (int lev, IntVect const& ng) +{ + MultiFab::Copy(m_leveldata[lev]->velocity_o, + m_leveldata[lev]->velocity, 0, 0, AMREX_SPACEDIM, ng); +} + +void incflo::copy_from_old_to_new_velocity (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + copy_from_old_to_new_velocity(lev, ng); + } +} + +void incflo::copy_from_old_to_new_velocity (int lev, IntVect const& ng) +{ + MultiFab::Copy(m_leveldata[lev]->velocity, + m_leveldata[lev]->velocity_o, 0, 0, AMREX_SPACEDIM, ng); +} + +void incflo::copy_from_new_to_old_density (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + copy_from_new_to_old_density(lev, ng); + } +} + +void incflo::copy_from_new_to_old_density (int lev, IntVect const& ng) +{ + MultiFab::Copy(m_leveldata[lev]->density_o, + m_leveldata[lev]->density, 0, 0, 1, ng); +} + +void incflo::copy_from_old_to_new_density (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + copy_from_old_to_new_density(lev, ng); + } +} + +void incflo::copy_from_old_to_new_density (int lev, IntVect const& ng) +{ + MultiFab::Copy(m_leveldata[lev]->density, + m_leveldata[lev]->density_o, 0, 0, 1, ng); +} + +void incflo::copy_from_new_to_old_tracer (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + copy_from_new_to_old_tracer(lev, ng); + } +} + +void incflo::copy_from_new_to_old_tracer (int lev, IntVect const& ng) +{ + if (m_ntrac > 0) { + MultiFab::Copy(m_leveldata[lev]->tracer_o, + m_leveldata[lev]->tracer, 0, 0, m_ntrac, ng); + } +} + +void incflo::copy_from_old_to_new_tracer (IntVect const& ng) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + copy_from_old_to_new_tracer(lev, ng); + } +} + +void incflo::copy_from_old_to_new_tracer (int lev, IntVect const& ng) +{ + if (m_ntrac > 0) { + MultiFab::Copy(m_leveldata[lev]->tracer, + m_leveldata[lev]->tracer_o, 0, 0, m_ntrac, ng); + } +}