diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a696aac4..f6602817 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -11,8 +11,12 @@ target_sources(incflo incflo_correct_small_cells.cpp incflo_redistribute.cpp incflo_explicit_update.cpp - incflo_tagging.cpp incflo_regrid.cpp + 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 e98e6a5a..135785eb 100644 --- a/src/Make.package +++ b/src/Make.package @@ -5,8 +5,12 @@ 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 += incflo_update_velocity.cpp +CEXE_sources += incflo_utils.cpp CEXE_sources += main.cpp ifeq ($(USE_EB), TRUE) diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 8f7fbcbf..0d64f979 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -107,21 +107,24 @@ incflo::compute_convective_term (Vector const& conv_u, for (int lev = 0; lev <= finest_level; ++lev) { AMREX_D_TERM( - face_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));, - face_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));, - face_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));); + face_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),u_mac[lev]->Factory());, + face_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),v_mac[lev]->Factory());, + face_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),w_mac[lev]->Factory());); AMREX_D_TERM( - flux_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));, - flux_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));, - flux_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));); + flux_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),u_mac[lev]->Factory());, + flux_y[lev].define(v_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),v_mac[lev]->Factory());, + flux_z[lev].define(w_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),w_mac[lev]->Factory());); - divu[lev].define(vel[lev]->boxArray(),dmap[lev],1,4,MFInfo(),Factory(lev)); - if (m_advect_momentum) + divu[lev].define(vel[lev]->boxArray(),dmap[lev],1,4,MFInfo(),vel[lev]->Factory()); + + if (m_advect_momentum) { rhovel[lev].define(vel[lev]->boxArray(),dmap[lev],AMREX_SPACEDIM, - vel[lev]->nGrow(),MFInfo(),Factory(lev)); - if (m_advect_tracer && m_ntrac > 0 && any_conserv_trac) - rhotrac[lev].define(vel[lev]->boxArray(),dmap[lev],tracer[lev]->nComp(), - tracer[lev]->nGrow(),MFInfo(),Factory(lev)); + vel[lev]->nGrow(),MFInfo(),vel[lev]->Factory()); + } + if (m_advect_tracer && m_ntrac > 0 && any_conserv_trac) { + rhotrac[lev].define(tracer[lev]->boxArray(),dmap[lev],tracer[lev]->nComp(), + tracer[lev]->nGrow(),MFInfo(),tracer[lev]->Factory()); + } AMREX_D_TERM(faces[lev][0] = &face_x[lev];, faces[lev][1] = &face_y[lev];, diff --git a/src/diffusion/DiffusionScalarOp.H b/src/diffusion/DiffusionScalarOp.H index a77ac786..df74be14 100644 --- a/src/diffusion/DiffusionScalarOp.H +++ b/src/diffusion/DiffusionScalarOp.H @@ -28,7 +28,6 @@ public: void compute_laps (amrex::Vector const& laps, amrex::Vector const& a_scalar, - amrex::Vector const& /*density*/, amrex::Vector const& eta); void compute_divtau (amrex::Vector const& a_divtau, diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp index d078c1de..032da0c1 100644 --- a/src/diffusion/DiffusionScalarOp.cpp +++ b/src/diffusion/DiffusionScalarOp.cpp @@ -474,7 +474,6 @@ DiffusionScalarOp::diffuse_vel_components (Vector const& vel, void DiffusionScalarOp::compute_laps (Vector const& a_laps, Vector const& a_scalar, - Vector const& /*a_density*/, Vector const& a_eta) { BL_PROFILE("DiffusionScalarOp::compute_laps"); diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp index d171d545..ae51d74f 100644 --- a/src/diffusion/incflo_diffusion.cpp +++ b/src/diffusion/incflo_diffusion.cpp @@ -57,10 +57,9 @@ incflo::compute_divtau(Vector const& divtau, void incflo::compute_laps(Vector const& laps, Vector const& scalar, - Vector const& density, Vector const& eta) { - get_diffusion_scalar_op()->compute_laps(laps, scalar, density, eta); + get_diffusion_scalar_op()->compute_laps(laps, scalar, eta); } void diff --git a/src/incflo.H b/src/incflo.H index 16061dba..3e229738 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,12 @@ 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_eta, + amrex::Vector& tra_forces); + void update_velocity (StepType step_type, amrex::Vector& vel_eta, + amrex::Vector& vel_forces); + /////////////////////////////////////////////////////////////////////////// // // derive @@ -201,11 +211,8 @@ public: [[nodiscard]] amrex::Array average_scalar_eta_to_faces (int lev, int comp, amrex::MultiFab const& cc_eta) const; - // FIXME? compute_laps never used density in any meaningful way, but - // then removing density would break anything downstream.... void compute_laps (amrex::Vector const& laps, amrex::Vector const& scalar, - amrex::Vector const& density, amrex::Vector const& eta); void diffuse_scalar (amrex::Vector const& scalar, @@ -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; @@ -745,6 +753,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; @@ -768,6 +777,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 da9d0300..9b0c6e3f 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 @@ -143,338 +133,35 @@ void incflo::ApplyCorrector() // ************************************************************************************* // Compute viscosity / diffusive coefficients // ************************************************************************************* - compute_viscosity(GetVecOfPtrs(vel_eta), - get_density_new(), get_velocity_new(), - new_time, 1); - compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 1); - // Here we create divtau of the (n+1,*) state that was computed in the predictor; - // we use this laps only if DiffusionType::Explicit + // Here we create divtau of the (n+1,*) state that was computed in the predictor if ( (m_diff_type == DiffusionType::Explicit) || use_tensor_correction ) { compute_divtau(get_divtau_new(), get_velocity_new_const(), 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. - // ************************************************************************************* - Real l_dt = m_dt; - - // ************************************************************************************* - // 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 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 - // ************************************************************************************* - if (m_advect_tracer) - { - tracer_explicit_update_corrector(tra_forces); - } - // ************************************************************************************* - // Solve diffusion equation for tracer + // Update density // ************************************************************************************* - 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_density(StepType::Corrector); // ************************************************************************************* - // Define the forcing terms to use in the final update (using half-time density) + // Update tracer // ************************************************************************************* - compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_new_const(), - GetVecOfConstPtrs(density_nph), - get_tracer_old_const(), get_tracer_new_const()); + update_tracer(StepType::Corrector, tra_eta, tra_forces); // ************************************************************************************* // 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); + update_velocity(StepType::Corrector, vel_eta, vel_forces); - 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); - } - - // ********************************************************************************************** - // // Project velocity field, update pressure + // ********************************************************************************************** bool incremental_projection = false; - ApplyProjection(GetVecOfConstPtrs(density_nph), new_time,m_dt,incremental_projection); - -#ifdef INCFLO_USE_PARTICLES - // ************************************************************************************** - // Update the particle positions - // ************************************************************************************** - evolveTracerParticles(AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac), - GetVecOfConstPtrs(w_mac))); -#endif + ApplyProjection(get_density_nph_const(), new_time,m_dt,incremental_projection); #ifdef AMREX_USE_EB // ********************************************************************************************** @@ -486,4 +173,12 @@ void incflo::ApplyCorrector() AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac), GetVecOfConstPtrs(w_mac))); #endif + +#ifdef INCFLO_USE_PARTICLES + // ************************************************************************************** + // Update the particle positions + // ************************************************************************************** + evolveTracerParticles(AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac), + GetVecOfConstPtrs(w_mac))); +#endif } diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index ca161c24..8039574c 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; @@ -99,10 +97,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; @@ -125,17 +119,12 @@ void incflo::ApplyPredictor (bool incremental_projection) } } - // ************************************************************************************* - // We now define the forcing terms to use in the Godunov prediction inside the predictor - // ************************************************************************************* - // ************************************************************************************* // Compute viscosity / diffusive coefficients // ************************************************************************************* 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 @@ -148,21 +137,27 @@ void incflo::ApplyPredictor (bool incremental_projection) } // ************************************************************************************* - // Compute explicit diffusive terms + // Compute explicit diffusive term -- note this is used inside compute_convective_term // ************************************************************************************* - if (m_advect_tracer && need_divtau()) + if (m_advect_tracer) { - compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(), - GetVecOfConstPtrs(tra_eta)); + compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + if (need_divtau()) { + compute_laps(get_laps_old(), get_tracer_old_const(), GetVecOfConstPtrs(tra_eta)); + } } // ********************************************************************************************** - // Compute the MAC-projected velocities at all levels + // Compute the forcing terms // ************************************************************************************* bool include_pressure_gradient = !(m_use_mac_phi_in_godunov); compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_old_const(), get_density_old_const(), get_tracer_old_const(), get_tracer_old_const(), include_pressure_gradient); + + // ********************************************************************************************** + // Compute the MAC-projected velocities at all levels + // ************************************************************************************* 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); @@ -172,6 +167,7 @@ void incflo::ApplyPredictor (bool incremental_projection) // Compute the explicit advective terms R_u^(n+1/2), R_s^(n+1/2) and R_t^(n+1/2) // if (advection_type == "MOL" ) // Compute the explicit advective terms R_u^n , R_s^n and R_t^n + // Note that if advection_type != "MOL" then we call compute_tra_forces inside this routine // ************************************************************************************* compute_convective_term(get_conv_velocity_old(), get_conv_density_old(), get_conv_tracer_old(), get_velocity_old_const(), get_density_old_const(), get_tracer_old_const(), @@ -181,269 +177,24 @@ void incflo::ApplyPredictor (bool incremental_projection) m_cur_time); // ************************************************************************************* - // Define local variables for lambda to capture. - // ************************************************************************************* - Real l_dt = m_dt; - + // Update density // ************************************************************************************* - // 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 + update_density(StepType::Predictor); - // 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 the 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) - - - // ************************************************************************************* - // 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 - // ************************************************************************************* - 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); - } + // ********************************************************************************************** + // Update tracer + // ********************************************************************************************** + update_tracer(StepType::Predictor, tra_eta, tra_forces); - 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 + // ********************************************************************************************** + 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 // ************************************************************************************** @@ -457,9 +208,7 @@ void incflo::ApplyPredictor (bool incremental_projection) #ifdef AMREX_USE_EB // ********************************************************************************************** - // // Over-write velocity in cells with vfrac < 1e-4 - // // ********************************************************************************************** if (m_advection_type == "MOL") incflo_correct_small_cells(get_velocity_new(), diff --git a/src/incflo_update_density.cpp b/src/incflo_update_density.cpp new file mode 100644 index 00000000..3ee93394 --- /dev/null +++ b/src/incflo_update_density.cpp @@ -0,0 +1,73 @@ +#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++) + { + 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& 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)); + }); + } + } // 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 + } + + for (int lev = 0; lev <= finest_level; lev++) + { + auto& ld = *m_leveldata[lev]; + + // Fill ghost cells of new-time density if needed (we assume ghost cells of old density are already filled) + if (ng > 0) { + fillpatch_density(lev, m_t_new[lev], ld.density, ng); + } + + // Define half-time density after the average down + MultiFab::LinComb(ld.density_nph, Real(0.5), ld.density, 0, Real(0.5), ld.density_o, 0, 0, 1, ng); + } + + } else { + for (int lev = 0; lev <= finest_level; lev++) { + MultiFab::Copy(m_leveldata[lev]->density_nph, m_leveldata[lev]->density_o, 0, 0, 1, ng); + } + } +} diff --git a/src/incflo_update_tracer.cpp b/src/incflo_update_tracer.cpp new file mode 100644 index 00000000..eb8dedad --- /dev/null +++ b/src/incflo_update_tracer.cpp @@ -0,0 +1,70 @@ +#include + +using namespace amrex; + +void incflo::update_tracer (StepType step_type, Vector& tra_eta, Vector& tra_forces) +{ + BL_PROFILE("incflo::update_tracer"); + + Real new_time = m_cur_time + m_dt; + + if (m_advect_tracer) + { + // ************************************************************************************* + // Compute diffusion coefficient for tracer + // ************************************************************************************* + + // ************************************************************************************* + // Compute the tracer forcing terms (forcing for (rho s), not for s) + // ************************************************************************************* + compute_tra_forces(GetVecOfPtrs(tra_forces), get_density_nph_const()); + + // ************************************************************************************* + // Compute explicit diffusive term (if corrector) + // ************************************************************************************* + if (step_type == StepType::Corrector) + { + compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1); + if (m_diff_type == DiffusionType::Explicit) { + compute_laps(get_laps_new(), get_tracer_new_const(), GetVecOfConstPtrs(tra_eta)); + } + } + + // ************************************************************************************* + // 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 + // ************************************************************************************* + 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/incflo_update_velocity.cpp b/src/incflo_update_velocity.cpp new file mode 100644 index 00000000..bd5a5375 --- /dev/null +++ b/src/incflo_update_velocity.cpp @@ -0,0 +1,341 @@ +#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& 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 = ld.density_nph.const_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 * ( + l_half*( dvdt_o(i,j,k,0)+ dvdt(i,j,k,0)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,1)+ dvdt(i,j,k,1)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,2)+ dvdt(i,j,k,2)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,0)+ dvdt(i,j,k,0)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,1)+ dvdt(i,j,k,1)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,2)+ dvdt(i,j,k,2)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + +l_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 * ( + l_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) + +l_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 * ( + l_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + +l_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 * ( + l_half*( dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,1)+dvdt(i,j,k,1)) + + l_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 * ( + l_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_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 * ( + l_half*( dvdt_o(i,j,k,2)+dvdt(i,j,k,2)) + vel_f(i,j,k,2) );); + }); + } + } + } + } + } // 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); + } +} 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)