From 3f9c0f5045c587fda424b0a847c7bbb314ed0d0f Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 24 Sep 2024 08:12:45 -0700 Subject: [PATCH] put in basic data structures / operations for cc-projection as alternative to nodal projection --- src/boundary_conditions/CMakeLists.txt | 1 + src/boundary_conditions/Make.package | 1 + .../incflo_set_velocity_bcs.cpp | 39 ++ src/incflo.H | 23 +- src/incflo.cpp | 11 +- src/incflo_apply_corrector.cpp | 4 +- src/incflo_apply_predictor.cpp | 4 +- src/incflo_regrid.cpp | 4 + src/prob/prob_bc.cpp | 96 ++++ src/prob/prob_init_fluid.cpp | 2 + src/projection/CMakeLists.txt | 1 + src/projection/Make.package | 1 + src/projection/incflo_apply_cc_projection.cpp | 467 ++++++++++++++++++ .../incflo_apply_nodal_projection.cpp | 4 +- src/projection/incflo_projection_bc.cpp | 4 +- src/setup/init.cpp | 59 ++- src/utilities/io.cpp | 38 +- 17 files changed, 732 insertions(+), 27 deletions(-) create mode 100644 src/boundary_conditions/incflo_set_velocity_bcs.cpp create mode 100644 src/projection/incflo_apply_cc_projection.cpp diff --git a/src/boundary_conditions/CMakeLists.txt b/src/boundary_conditions/CMakeLists.txt index ce9bd6251..7bf0c1725 100644 --- a/src/boundary_conditions/CMakeLists.txt +++ b/src/boundary_conditions/CMakeLists.txt @@ -4,4 +4,5 @@ target_sources(incflo incflo_fillpatch.cpp incflo_fillphysbc.cpp incflo_set_bcs.cpp + incflo_set_velocity_bcs.cpp ) diff --git a/src/boundary_conditions/Make.package b/src/boundary_conditions/Make.package index 9de2ae895..8e382a0ff 100644 --- a/src/boundary_conditions/Make.package +++ b/src/boundary_conditions/Make.package @@ -2,3 +2,4 @@ CEXE_sources += boundary_conditions.cpp CEXE_sources += incflo_fillpatch.cpp incflo_fillphysbc.cpp CEXE_sources += incflo_set_bcs.cpp +CEXE_sources += incflo_set_velocity_bcs.cpp diff --git a/src/boundary_conditions/incflo_set_velocity_bcs.cpp b/src/boundary_conditions/incflo_set_velocity_bcs.cpp new file mode 100644 index 000000000..31cb74bd5 --- /dev/null +++ b/src/boundary_conditions/incflo_set_velocity_bcs.cpp @@ -0,0 +1,39 @@ +#ifdef AMREX_USE_EB +#include +#endif + +#include + +using namespace amrex; + +void +incflo::set_inflow_velocity (int lev, amrex::Real time, MultiFab& vel, int nghost) +{ + Geometry const& gm = Geom(lev); + Box const& domain = gm.growPeriodicDomain(nghost); + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + Orientation olo(dir,Orientation::low); + Orientation ohi(dir,Orientation::high); + if (m_bc_type[olo] == BC::mass_inflow || m_bc_type[ohi] == BC::mass_inflow) { + Box dlo = (m_bc_type[olo] == BC::mass_inflow) ? amrex::adjCellLo(domain,dir,nghost) : Box(); + Box dhi = (m_bc_type[ohi] == BC::mass_inflow) ? amrex::adjCellHi(domain,dir,nghost) : Box(); +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(vel); mfi.isValid(); ++mfi) { + Box const& gbx = amrex::grow(mfi.validbox(),nghost); + Box blo = gbx & dlo; + Box bhi = gbx & dhi; + Array4 const& v = vel[mfi].array(); + int gid = mfi.index(); + if (blo.ok()) { + prob_set_inflow_velocity(gid, olo, blo, v, lev, time); + } + if (bhi.ok()) { + prob_set_inflow_velocity(gid, ohi, bhi, v, lev, time); + } + } + } + } + vel.EnforcePeriodicity(gm.periodicity()); +} diff --git a/src/incflo.H b/src/incflo.H index 08670e23a..d419d7cd4 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -123,6 +123,8 @@ public: amrex::iMultiFab make_nodalBC_mask (int lev); amrex::Vector make_robinBC_MFs(int lev, amrex::MultiFab* state = nullptr); + void set_inflow_velocity (int lev, amrex::Real time, amrex::MultiFab& vel, int nghost); + #ifdef AMREX_USE_EB void set_eb_velocity (int lev, amrex::Real time, amrex::MultiFab& eb_vel, int nghost); void set_eb_density (int lev, amrex::Real time, amrex::MultiFab& eb_density, int nghost); @@ -249,6 +251,9 @@ public: amrex::Array4 const& bcval, int lev); + void prob_set_inflow_velocity (int grid_id, amrex::Orientation ori, amrex::Box const& bx, + amrex::Array4 const& v, int lev, amrex::Real time); + #include "incflo_prob_I.H" #include "incflo_prob_usr_I.H" @@ -259,6 +264,9 @@ public: /////////////////////////////////////////////////////////////////////////// void ApplyProjection(amrex::Vector density, + AMREX_D_DECL(amrex::Vector const& u_mac, + amrex::Vector const& v_mac, + amrex::Vector const& w_mac), amrex::Real time, amrex::Real scaling_factor, bool incremental); void ApplyNodalProjection(amrex::Vector density, amrex::Real time, amrex::Real scaling_factor, bool incremental); @@ -267,6 +275,11 @@ public: amrex::Vector const& divu_Source, amrex::Real time, amrex::Real scaling_factor, bool incremental, bool set_inflow_bc); + void ApplyCCProjection(amrex::Vector density, + AMREX_D_DECL(amrex::Vector const& u_mac, + amrex::Vector const& v_mac, + amrex::Vector const& w_mac), + amrex::Real time, amrex::Real scaling_factor, bool incremental); /////////////////////////////////////////////////////////////////////////// // @@ -508,6 +521,9 @@ private: // use gradient of mac phi which contains the full pressure bool m_use_mac_phi_in_godunov = false; + // If true use CC projection; if false use nodal projection + bool m_use_cc_proj = false; + enum struct DiffusionType { Invalid, Explicit, Crank_Nicolson, Implicit }; @@ -554,7 +570,7 @@ private: int m_plt_gpz = 1; int m_plt_rho = 1; int m_plt_tracer = 1; - int m_plt_p_nd = 0; + int m_plt_p = 0; int m_plt_macphi = 0; int m_plt_eta = 0; int m_plt_magvel = 1; @@ -615,7 +631,8 @@ private: amrex::MultiFab mac_phi; // cell-centered pressure used in MAC projection - // nodal pressure multifab + // Pressure MultiFabs (only one will actually be used) + amrex::MultiFab p_cc; amrex::MultiFab p_nd; // cell-centered pressure gradient @@ -819,7 +836,7 @@ private: return m_bcrec_force_d.data(); } [[nodiscard]] amrex::Array - get_nodal_projection_bc (amrex::Orientation::Side side) const noexcept; + get_projection_bc (amrex::Orientation::Side side) const noexcept; [[nodiscard]] amrex::Array get_mac_projection_bc (amrex::Orientation::Side side) const noexcept; diff --git a/src/incflo.cpp b/src/incflo.cpp index 33823686a..9a07446cd 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -197,10 +197,19 @@ void incflo::Evolve() void incflo::ApplyProjection (Vector density, + AMREX_D_DECL(Vector const& u_mac, + Vector const& v_mac, + Vector const& w_mac), Real time, Real scaling_factor, bool incremental) { BL_PROFILE("incflo::ApplyProjection"); - ApplyNodalProjection(std::move(density),time,scaling_factor,incremental); + if (m_use_cc_proj) + { + ApplyCCProjection(density,AMREX_D_DECL(u_mac,v_mac,w_mac), + time,scaling_factor,incremental); + } + else + ApplyNodalProjection(density,time,scaling_factor,incremental); } // Make a new level from scratch using provided BoxArray and DistributionMapping. diff --git a/src/incflo_apply_corrector.cpp b/src/incflo_apply_corrector.cpp index 9b0c6e3f6..4d00cae43 100644 --- a/src/incflo_apply_corrector.cpp +++ b/src/incflo_apply_corrector.cpp @@ -161,7 +161,9 @@ void incflo::ApplyCorrector() // Project velocity field, update pressure // ********************************************************************************************** bool incremental_projection = false; - ApplyProjection(get_density_nph_const(), new_time,m_dt,incremental_projection); + ApplyProjection(get_density_nph_const(), + AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), + GetVecOfPtrs(w_mac)),new_time,m_dt,incremental_projection); #ifdef AMREX_USE_EB // ********************************************************************************************** diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index 8039574cc..58296b862 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -194,7 +194,9 @@ void incflo::ApplyPredictor (bool incremental_projection) // ********************************************************************************************** // Project velocity field, update pressure // ********************************************************************************************** - ApplyProjection(get_density_nph_const(),new_time,m_dt,incremental_projection); + ApplyProjection(get_density_nph_const(), + AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac), + GetVecOfPtrs(w_mac)),new_time,m_dt,incremental_projection); #ifdef INCFLO_USE_PARTICLES // ************************************************************************************** diff --git a/src/incflo_regrid.cpp b/src/incflo_regrid.cpp index afaf5f863..618f1b8b1 100644 --- a/src/incflo_regrid.cpp +++ b/src/incflo_regrid.cpp @@ -34,6 +34,8 @@ void incflo::MakeNewLevelFromCoarse (int lev, fillcoarsepatch_tracer(lev, time, new_leveldata->tracer, 0); } fillcoarsepatch_gradp(lev, time, new_leveldata->gp, 0); + + new_leveldata->p_cc.setVal(0.0); new_leveldata->p_nd.setVal(0.0); m_leveldata[lev] = std::move(new_leveldata); @@ -82,6 +84,8 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba, fillpatch_tracer(lev, time, new_leveldata->tracer, 0); } fillpatch_gradp(lev, time, new_leveldata->gp, 0); + + new_leveldata->p_cc.setVal(0.0); new_leveldata->p_nd.setVal(0.0); m_leveldata[lev] = std::move(new_leveldata); diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp index 2545f635f..a2fdd275e 100644 --- a/src/prob/prob_bc.cpp +++ b/src/prob/prob_bc.cpp @@ -296,3 +296,99 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx, +std::to_string(m_probtype)); } } + +void incflo::prob_set_inflow_velocity (int /*grid_id*/, Orientation ori, Box const& bx, + Array4 const& vel, int lev, Real /*time*/) +{ + if (6 == m_probtype) + { + AMREX_D_TERM(Real u = m_ic_u;, + Real v = m_ic_v;, + Real w = m_ic_w;); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(vel(i,j,k,0) = u;, + vel(i,j,k,1) = v;, + vel(i,j,k,2) = w;); + }); + } + else if (31 == m_probtype) + { + Real dyinv = 1.0 / Geom(lev).Domain().length(1); + Real u = m_ic_u; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real y = (j+0.5)*dyinv; + vel(i,j,k,0) = 6. * u * y * (1.-y); + }); + } + else if (311 == m_probtype) + { + Real dzinv = 1.0 / Geom(lev).Domain().length(2); + Real u = m_ic_u; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real z = (k+0.5)*dzinv; + vel(i,j,k,0) = 6. * u * z * (1.-z); + }); + } + else if (41 == m_probtype) + { + Real dzinv = 1.0 / Geom(lev).Domain().length(2); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real z = (k+0.5)*dzinv; + vel(i,j,k,0) = 0.5*z; + }); + } + else if (32 == m_probtype) + { + Real dzinv = 1.0 / Geom(lev).Domain().length(2); + Real v = m_ic_v; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real z = (k+0.5)*dzinv; + vel(i,j,k,1) = 6. * v * z * (1.-z); + }); + } + else if (322 == m_probtype) + { + Real dxinv = 1.0 / Geom(lev).Domain().length(0); + Real v = m_ic_v; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real x = (i+0.5)*dxinv; + vel(i,j,k,1) = 6. * v * x * (1.-x); + }); + } + else if (33 == m_probtype) + { + Real dxinv = 1.0 / Geom(lev).Domain().length(0); + Real w = m_ic_w; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real x = (i+0.5)*dxinv; + vel(i,j,k,2) = 6. * w * x * (1.-x); + }); + } + else if (333 == m_probtype) + { + Real dyinv = 1.0 / Geom(lev).Domain().length(1); + Real w = m_ic_w; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real y = (j+0.5)*dyinv; + vel(i,j,k,2) = 6. * w * y * (1.-y); + }); + } + else + { + const int dir = ori.coordDir(); + const Real bcv = m_bc_velocity[ori][dir]; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + vel(i,j,k,dir) = bcv; + }); + }; +} diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index 0e50205f5..7c2f5e092 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -10,7 +10,9 @@ void incflo::prob_init_fluid (int lev) auto const& problo = geom[lev].ProbLoArray(); auto const& probhi = geom[lev].ProbHiArray(); + ld.p_cc.setVal(0.0); ld.p_nd.setVal(0.0); + ld.gp.setVal(0.0); ld.density.setVal(m_ro_0); diff --git a/src/projection/CMakeLists.txt b/src/projection/CMakeLists.txt index e5139323f..9484857a6 100644 --- a/src/projection/CMakeLists.txt +++ b/src/projection/CMakeLists.txt @@ -3,5 +3,6 @@ target_include_directories(incflo PRIVATE ${CMAKE_CURRENT_LIST_DIR}) target_sources(incflo PRIVATE incflo_projection_bc.cpp + incflo_apply_cc_projection.cpp incflo_apply_nodal_projection.cpp ) diff --git a/src/projection/Make.package b/src/projection/Make.package index fe856079b..2c57d6bab 100644 --- a/src/projection/Make.package +++ b/src/projection/Make.package @@ -1,2 +1,3 @@ CEXE_sources += incflo_projection_bc.cpp +CEXE_sources += incflo_apply_cc_projection.cpp CEXE_sources += incflo_apply_nodal_projection.cpp diff --git a/src/projection/incflo_apply_cc_projection.cpp b/src/projection/incflo_apply_cc_projection.cpp new file mode 100644 index 000000000..8dff4f09a --- /dev/null +++ b/src/projection/incflo_apply_cc_projection.cpp @@ -0,0 +1,467 @@ +#include +#include +#include + +using namespace amrex; + +void +average_ccvel_to_mac (const Array& fc, const MultiFab& cc) +{ + AMREX_ASSERT(cc.nComp() == AMREX_SPACEDIM); + AMREX_ASSERT(cc.nGrow() >= 1); + AMREX_ASSERT(fc[0]->nComp() == 1); +#if (AMREX_SPACEDIM >= 2) + AMREX_ASSERT(fc[1]->nComp() == 1); +#endif +#if (AMREX_SPACEDIM == 3) + AMREX_ASSERT(fc[2]->nComp() == 1); +#endif + + int ncomp = AMREX_SPACEDIM; + +#ifdef AMREX_USE_GPU + if (Gpu::inLaunchRegion() && cc.isFusingCandidate()) { + auto const& ccma = cc.const_arrays(); + AMREX_D_TERM(auto const& fxma = fc[0]->arrays();, + auto const& fyma = fc[1]->arrays();, + auto const& fzma = fc[2]->arrays();); + MultiFab foo(amrex::convert(cc.boxArray(),IntVect(1)), cc.DistributionMap(), 1, 0, + MFInfo().SetAlloc(false)); + IntVect ng = -cc.nGrowVect(); + ParallelFor(foo, IntVect(0), ncomp, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept + { + Box ccbx(ccma[box_no]); + ccbx.grow(ng); + AMREX_D_TERM(Box const& xbx = amrex::surroundingNodes(ccbx,0);, + Box const& ybx = amrex::surroundingNodes(ccbx,1);, + Box const& zbx = amrex::surroundingNodes(ccbx,2);); + if (xbx.contains(i,j,k) and n == 0) { + fxma[box_no](i,j,k) = Real(0.5)*(ccma[box_no](i-1,j,k,n) + ccma[box_no](i,j,k,n)); + } + if (ybx.contains(i,j,k) and n == 1) { + fyma[box_no](i,j,k) = Real(0.5)*(ccma[box_no](i,j-1,k,n) + ccma[box_no](i,j,k,n)); + } +#if (AMREX_SPACEDIM == 3) + if (zbx.contains(i,j,k) and n == 2) { + fzma[box_no](i,j,k) = Real(0.5)*(ccma[box_no](i,j,k-1,n) + ccma[box_no](i,j,k,n)); + } +#endif + }); + Gpu::streamSynchronize(); + } else +#endif + { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(cc,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + AMREX_D_TERM(const Box& xbx = mfi.nodaltilebox(0);, + const Box& ybx = mfi.nodaltilebox(1);, + const Box& zbx = mfi.nodaltilebox(2);); + const auto& index_bounds = amrex::getIndexBounds(AMREX_D_DECL(xbx,ybx,zbx)); + + AMREX_D_TERM(Array4 const& fxarr = fc[0]->array(mfi);, + Array4 const& fyarr = fc[1]->array(mfi);, + Array4 const& fzarr = fc[2]->array(mfi);); + Array4 const& ccarr = cc.const_array(mfi); + + AMREX_HOST_DEVICE_PARALLEL_FOR_4D(index_bounds, ncomp, i, j, k, n, + { + if (xbx.contains(i,j,k) and n == 0) { + fxarr(i,j,k) = Real(0.5)*(ccarr(i-1,j,k,n) + ccarr(i,j,k,n)); + } + if (ybx.contains(i,j,k) and n == 1) { + fyarr(i,j,k) = Real(0.5)*(ccarr(i,j-1,k,n) + ccarr(i,j,k,n)); + } +#if (AMREX_SPACEDIM == 3) + if (zbx.contains(i,j,k) and n == 2) { + fzarr(i,j,k) = Real(0.5)*(ccarr(i,j,k-1,n) + ccarr(i,j,k,n)); + } +#endif + }); + } + } +} + +void +average_mac_to_ccvel (const Array& fc, MultiFab& cc) +{ + AMREX_ASSERT(cc.nComp() == AMREX_SPACEDIM); + AMREX_ASSERT(fc[0]->nComp() == 1); +#if (AMREX_SPACEDIM >= 2) + AMREX_ASSERT(fc[1]->nComp() == 1); +#endif +#if (AMREX_SPACEDIM == 3) + AMREX_ASSERT(fc[2]->nComp() == 1); +#endif + + int ncomp = AMREX_SPACEDIM; + +#ifdef AMREX_USE_GPU + if (Gpu::inLaunchRegion() && cc.isFusingCandidate()) { + auto const& ccma = cc.const_arrays(); + AMREX_D_TERM(auto const& fxma = fc[0]->arrays();, + auto const& fyma = fc[1]->arrays();, + auto const& fzma = fc[2]->arrays();); + MultiFab foo(amrex::convert(cc.boxArray(),IntVect(1)), cc.DistributionMap(), 1, 0, + MFInfo().SetAlloc(false)); + ParallelFor(foo, IntVect(0), ncomp, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept + { + Box ccbx(ccma[box_no]); + AMREX_D_TERM(Box const& xbx = amrex::surroundingNodes(ccbx,0);, + Box const& ybx = amrex::surroundingNodes(ccbx,1);, + Box const& zbx = amrex::surroundingNodes(ccbx,2);); + if (xbx.contains(i,j,k) and n == 0) { + ccma[box_no](i,j,k,n) = Real(0.5) * (fxma[box_no](i,j,k) + fxma[box_no](i+1,j,k)); + } + if (ybx.contains(i,j,k) and n == 1) { + ccma[box_no](i,j,k,n) = Real(0.5) * (fyma[box_no](i,j,k) + fyma[box_no](i,j+1,k)); + } +#if (AMREX_SPACEDIM == 3) + if (zbx.contains(i,j,k) and n == 2) { + ccma[box_no](i,j,k,n) = Real(0.5) * (fzma[box_no](i,j,k) + fzma[box_no](i,j,k+1)); + } +#endif + }); + Gpu::streamSynchronize(); + } else +#endif + { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(cc,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + AMREX_D_TERM(Array4 const& fxarr = fc[0]->array(mfi);, + Array4 const& fyarr = fc[1]->array(mfi);, + Array4 const& fzarr = fc[2]->array(mfi);); + Array4 const& ccarr = cc.array(mfi); + + AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n, + { + if (n == 0) + ccarr(i,j,k,0) = Real(0.5) * (fxarr(i,j,k) + fxarr(i+1,j,k)); + + if (n == 1) + ccarr(i,j,k,1) = Real(0.5) * (fyarr(i,j,k) + fyarr(i,j+1,k)); + +#if (AMREX_SPACEDIM == 3) + if (n == 2) + ccarr(i,j,k,2) = Real(0.5) * (fzarr(i,j,k) + fzarr(i,j,k+1)); +#endif + }); + } + } +} + +// +// Computes the following decomposition: +// +// u + dt grad(phi) / ro = u*, where div(u) = 0 +// +// where u* is a non-div-free velocity field, stored +// by components in u, v, and w. The resulting div-free +// velocity field, u, overwrites the value of u* in u, v, and w. +// +// phi is an auxiliary function related to the pressure p by the relation: +// +// new p = phi +// +// except in the initial projection when +// +// new p = old p + phi (nstep has its initial value -1) +// +// Note: scaling_factor equals dt except when called during initial projection, when it is 1.0 +// +void incflo::ApplyCCProjection (Vector density, + AMREX_D_DECL(Vector const& u_mac, + Vector const& v_mac, + Vector const& w_mac), + Real time, Real scaling_factor, bool incremental) +{ + BL_PROFILE("incflo::ApplyCCProjection"); + + // If we have dropped the dt substantially for whatever reason, + // use a different form of the approximate projection that + // projects (U^*-U^n + dt Gp) rather than (U^* + dt Gp) + + bool proj_for_small_dt = (time > 0.0 && m_dt < 0.1 * m_prev_dt); + + // Add ( dt grad p /ro ) to u* + if (!incremental) + { + 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& u = ld.velocity.array(mfi); + Array4 const& gp = ld.gp.const_array(mfi); + Array4 const& rho = density[lev]->const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real soverrho = scaling_factor / rho(i,j,k); + AMREX_D_TERM(u(i,j,k,0) += gp(i,j,k,0) * soverrho;, + u(i,j,k,1) += gp(i,j,k,1) * soverrho;, + u(i,j,k,2) += gp(i,j,k,2) * soverrho;); + }); + } + } + } + + // Define "vel" to be U^* - U^n rather than U^* + if (proj_for_small_dt || incremental) + { + for (int lev = 0; lev <= finest_level; ++lev) { + MultiFab::Subtract(m_leveldata[lev]->velocity, + m_leveldata[lev]->velocity_o, 0, 0, AMREX_SPACEDIM, 0); + } + } + + auto bclo = get_projection_bc(Orientation::low); + auto bchi = get_projection_bc(Orientation::high); + + Vector vel; + for (int lev = 0; lev <= finest_level; ++lev) { + vel.push_back(&(m_leveldata[lev]->velocity)); + vel[lev]->setBndry(0.0); + if (!proj_for_small_dt && !incremental) { + set_inflow_velocity(lev, time, *vel[lev], 1); + } + } + + // *************************************************************************************** + // START OF MAC STUFF + // *************************************************************************************** + + // This will hold (dt/rho) on faces + Vector > inv_rho(finest_level+1); + + for (int lev = 0; lev <= finest_level; ++lev) { + AMREX_D_TERM( + inv_rho[lev][0].define(u_mac[lev]->boxArray(),dmap[lev],1,0,MFInfo(),Factory(lev));, + inv_rho[lev][1].define(v_mac[lev]->boxArray(),dmap[lev],1,0,MFInfo(),Factory(lev));, + inv_rho[lev][2].define(w_mac[lev]->boxArray(),dmap[lev],1,0,MFInfo(),Factory(lev))); + AMREX_D_TERM( + inv_rho[lev][0].setVal(0.);, + inv_rho[lev][1].setVal(0.);, + inv_rho[lev][2].setVal(0.);); + } + + // Compute (1/rho) on faces + for (int lev=0; lev <= finest_level; ++lev) + { +#ifdef AMREX_USE_EB + EB_interp_CellCentroid_to_FaceCentroid (*density[lev], GetArrOfPtrs(inv_rho[lev]), + 0, 0, 1, geom[lev], get_density_bcrec()); +#else + amrex::average_cellcenter_to_face(GetArrOfPtrs(inv_rho[lev]), *density[lev], geom[lev]); +#endif + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + inv_rho[lev][idim].invert(scaling_factor, 0); + } + } + + // + // Initialize (or redefine the beta in) the MacProjector + if (macproj->needInitialization()) + { + LPInfo lp_info; + lp_info.setMaxCoarseningLevel(m_mac_mg_max_coarsening_level); +#ifndef AMREX_USE_EB + if (m_constant_density) { + Vector ba; + Vector dm; + for (auto const& ir : inv_rho) { + ba.push_back(ir[0].boxArray()); + dm.push_back(ir[0].DistributionMap()); + } + macproj->initProjector(ba, dm, lp_info, scaling_factor/m_ro_0); + } else +#endif + { + macproj->initProjector(lp_info, GetVecOfArrOfConstPtrs(inv_rho)); + } + macproj->setDomainBC(bclo, bchi); + } else { +#ifndef AMREX_USE_EB + if (m_constant_density) { + macproj->updateBeta(scaling_factor/m_ro_0); // unnecessary unless m_ro_0 changes. + } else +#endif + { + macproj->updateBeta(GetVecOfArrOfConstPtrs(inv_rho)); + } + } + + Vector cc_phi(finest_level+1); + Vector cc_gphi(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev ) + { + cc_phi[lev] = new MultiFab(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev)); + cc_gphi[lev] = new MultiFab(grids[lev], dmap[lev], AMREX_SPACEDIM, 0, MFInfo(), Factory(lev)); + cc_phi[lev]->setVal(0.); + cc_gphi[lev]->setVal(0.); + } + + Vector > mac_vec(finest_level+1); + for (int lev=0; lev <= finest_level; ++lev) + { + AMREX_D_TERM(mac_vec[lev][0] = u_mac[lev];, + mac_vec[lev][1] = v_mac[lev];, + mac_vec[lev][2] = w_mac[lev];); + } + + // Compute velocity on faces + for (int lev = 0; lev <= finest_level; ++lev) + { + // Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac} + // returned from this call are on face CENTROIDS + vel[lev]->FillBoundary(geom[lev].periodicity()); +#if 1 + MOL::ExtrapVelToFaces(*vel[lev], + AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]), + geom[lev], + get_velocity_bcrec(), get_velocity_bcrec_device_ptr()); + +#else + average_ccvel_to_mac( mac_vec[lev], *vel[lev]); +#endif + } + + macproj->setUMAC(mac_vec); + + if (m_verbose > 2) amrex::Print() << "CC Projection:\n"; + // + // Perform MAC projection: - del dot (dt/rho) grad phi = div(U) + // + macproj->project(cc_phi,m_mac_mg_rtol,m_mac_mg_atol); + + // + // After the projection we grab the dt/rho (grad phi) used in the projection + // + Vector > m_fluxes; + m_fluxes.resize(finest_level+1); + for (int lev=0; lev <= finest_level; ++lev) + { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + m_fluxes[lev][idim].define( + amrex::convert(grids[lev], IntVect::TheDimensionVector(idim)), + dmap[lev], 1, 0, MFInfo(), Factory(lev)); + } + } + + // + // Note that "fluxes" comes back as MINUS (dt/rho) Gphi + // +#ifdef AMREX_USE_EB + macproj->getFluxes(amrex::GetVecOfArrOfPtrs(m_fluxes), cc_phi, MLMG::Location::FaceCentroid); +#else + macproj->getFluxes(amrex::GetVecOfArrOfPtrs(m_fluxes), cc_phi, MLMG::Location::FaceCenter); +#endif + + for (int lev=0; lev <= finest_level; ++lev) + { +#ifdef AMREX_USE_EB + amrex::Abort("Haven't written mac_to_ccvel for EB"); +#else + average_mac_to_ccvel(GetArrOfPtrs(m_fluxes[lev]),*cc_gphi[lev]); +#endif + } + + 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.gp,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box const& tbx = mfi.tilebox(); + Array4 const& gp_cc = ld.gp.array(mfi); + Array4 const& p_cc = ld.p_cc.array(mfi); + Array4 const& gphi = cc_gphi[lev]->const_array(mfi); + Array4 const& phi = cc_phi[lev]->const_array(mfi); + + Array4 const& u = ld.velocity.array(mfi); + Array4 const& rho = density[lev]->const_array(mfi); + + amrex::ParallelFor(tbx, [u,gphi,p_cc,phi,incremental] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(u(i,j,k,0) += gphi(i,j,k,0);, + u(i,j,k,1) += gphi(i,j,k,1);, + u(i,j,k,2) += gphi(i,j,k,2);); + if (incremental) + p_cc (i,j,k) += phi(i,j,k); + else + p_cc (i,j,k) = phi(i,j,k); + }); + + if (incremental && m_constant_density) { + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(gp_cc(i,j,k,0) -= gphi(i,j,k,0) * m_ro_0 / scaling_factor;, + gp_cc(i,j,k,1) -= gphi(i,j,k,1) * m_ro_0 / scaling_factor;, + gp_cc(i,j,k,2) -= gphi(i,j,k,2) * m_ro_0 / scaling_factor;); + + }); + } else if (incremental && !m_constant_density) { + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(gp_cc(i,j,k,0) -= gphi(i,j,k,0) * rho(i,j,k) / scaling_factor;, + gp_cc(i,j,k,1) -= gphi(i,j,k,1) * rho(i,j,k) / scaling_factor;, + gp_cc(i,j,k,2) -= gphi(i,j,k,2) * rho(i,j,k) / scaling_factor;); + }); + } else if (!incremental && m_constant_density) { + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(gp_cc(i,j,k,0) = -gphi(i,j,k,0) * m_ro_0 / scaling_factor;, + gp_cc(i,j,k,1) = -gphi(i,j,k,1) * m_ro_0 / scaling_factor;, + gp_cc(i,j,k,2) = -gphi(i,j,k,2) * m_ro_0 / scaling_factor;); + }); + } else if (!incremental && !m_constant_density) { + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + AMREX_D_TERM(gp_cc(i,j,k,0) = -gphi(i,j,k,0) * rho(i,j,k) / scaling_factor;, + gp_cc(i,j,k,1) = -gphi(i,j,k,1) * rho(i,j,k) / scaling_factor;, + gp_cc(i,j,k,2) = -gphi(i,j,k,2) * rho(i,j,k) / scaling_factor;); + }); + } + } + ld.gp.FillBoundary(geom[lev].periodicity()); + ld.p_cc.FillBoundary(geom[lev].periodicity()); + } + + // *************************************************************************************** + // END OF MAC STUFF + // *************************************************************************************** + + // Define "vel" to be U^{n+1} rather than (U^{n+1}-U^n) + if (proj_for_small_dt || incremental) + { + for (int lev = 0; lev <= finest_level; ++lev) { + MultiFab::Add(m_leveldata[lev]->velocity, + m_leveldata[lev]->velocity_o, 0, 0, AMREX_SPACEDIM, 0); + } + } + + for (int lev = finest_level-1; lev >= 0; --lev) { +#ifdef AMREX_USE_EB + amrex::EB_average_down(m_leveldata[lev+1]->gp, m_leveldata[lev]->gp, + 0, AMREX_SPACEDIM, refRatio(lev)); +#else + amrex::average_down(m_leveldata[lev+1]->gp, m_leveldata[lev]->gp, + 0, AMREX_SPACEDIM, refRatio(lev)); +#endif + } +} diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index 57a18043b..df1da3c46 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -123,8 +123,8 @@ void incflo::ApplyNodalProjection (Vector const& density, // Perform projection std::unique_ptr nodal_projector; - auto bclo = get_nodal_projection_bc(Orientation::low); - auto bchi = get_nodal_projection_bc(Orientation::high); + auto bclo = get_projection_bc(Orientation::low); + auto bchi = get_projection_bc(Orientation::high); for (int lev = 0; lev <= finest_level; ++lev) { #ifdef AMREX_USE_EB diff --git a/src/projection/incflo_projection_bc.cpp b/src/projection/incflo_projection_bc.cpp index 2b038fdc1..c52642209 100644 --- a/src/projection/incflo_projection_bc.cpp +++ b/src/projection/incflo_projection_bc.cpp @@ -3,7 +3,7 @@ using namespace amrex; Array -incflo::get_nodal_projection_bc (Orientation::Side side) const noexcept +incflo::get_projection_bc (Orientation::Side side) const noexcept { amrex::Array r; for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { @@ -33,7 +33,7 @@ incflo::get_nodal_projection_bc (Orientation::Side side) const noexcept break; } default: - amrex::Abort("get_nodal_projection_bc: undefined BC type"); + amrex::Abort("get_projection_bc: undefined BC type"); }; } } diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 2a785bef3..95f6f3177 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -246,7 +246,7 @@ void incflo::ReadIOParameters() m_plt_gpz = 1; m_plt_rho = 1; m_plt_tracer = 1; - m_plt_p_nd = 0; + m_plt_p = 0; m_plt_macphi = 0; m_plt_eta = 0; m_plt_vort = 0; @@ -271,7 +271,7 @@ void incflo::ReadIOParameters() pp.query("plt_rho", m_plt_rho ); pp.query("plt_tracer", m_plt_tracer); - pp.query("plt_p_nd", m_plt_p_nd ); + pp.query("plt_p ", m_plt_p ); pp.query("plt_macphi", m_plt_macphi); pp.query("plt_eta", m_plt_eta ); pp.query("plt_magvel", m_plt_magvel); @@ -348,18 +348,43 @@ void incflo::InitialProjection() { BL_PROFILE("incflo::InitialProjection()"); + // ************************************************************************************* + // Allocate space for the temporary MAC velocities + // ************************************************************************************* + Vector u_mac_tmp(finest_level+1), v_mac_tmp(finest_level+1), w_mac_tmp(finest_level+1); + int ngmac = nghost_mac(); + + for (int lev = 0; lev <= finest_level; ++lev) { + AMREX_D_TERM(u_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(0)), dmap[lev], + 1, ngmac, MFInfo(), Factory(lev));, + v_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(1)), dmap[lev], + 1, ngmac, MFInfo(), Factory(lev));, + w_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(2)), dmap[lev], + 1, ngmac, MFInfo(), Factory(lev));); + if (ngmac > 0) { + AMREX_D_TERM(u_mac_tmp[lev].setBndry(0.0);, + v_mac_tmp[lev].setBndry(0.0);, + w_mac_tmp[lev].setBndry(0.0);); + } + } + Real dummy_dt = 1.0; - bool incremental = false; + bool incremental_projection = false; for (int lev = 0; lev <= finest_level; lev++) { m_leveldata[lev]->density.FillBoundary(geom[lev].periodicity()); } - ApplyProjection(get_density_new_const(),m_cur_time,dummy_dt,incremental); + + ApplyProjection(get_density_new_const(), + AMREX_D_DECL(GetVecOfPtrs(u_mac_tmp), GetVecOfPtrs(v_mac_tmp), + GetVecOfPtrs(w_mac_tmp)),m_cur_time,dummy_dt,incremental_projection); + // We set p and gp back to zero (p0 may still be still non-zero) for (int lev = 0; lev <= finest_level; lev++) { m_leveldata[lev]->p_nd.setVal(0.0); + m_leveldata[lev]->p_cc.setVal(0.0); m_leveldata[lev]->gp.setVal(0.0); } } @@ -381,6 +406,26 @@ void incflo::InitialPressureProjection() m_leveldata[lev]->density.FillBoundary(geom[lev].periodicity()); } + // ************************************************************************************* + // Allocate space for the temporary MAC velocities + // ************************************************************************************* + Vector u_mac_tmp(finest_level+1), v_mac_tmp(finest_level+1), w_mac_tmp(finest_level+1); + int ngmac = nghost_mac(); + + for (int lev = 0; lev <= finest_level; ++lev) { + AMREX_D_TERM(u_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(0)), dmap[lev], + 1, ngmac, MFInfo(), Factory(lev));, + v_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(1)), dmap[lev], + 1, ngmac, MFInfo(), Factory(lev));, + w_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(2)), dmap[lev], + 1, ngmac, MFInfo(), Factory(lev));); + if (ngmac > 0) { + AMREX_D_TERM(u_mac_tmp[lev].setBndry(0.0);, + v_mac_tmp[lev].setBndry(0.0);, + w_mac_tmp[lev].setBndry(0.0);); + } + } + // Set the velocity to the gravity field Vector vel(finest_level + 1); for (int lev = 0; lev <= finest_level; ++lev) { @@ -416,9 +461,9 @@ void incflo::InitialPressureProjection() // Always zero this here Vector Source(finest_level+1, nullptr); - ApplyNodalProjection(get_density_new_const(), GetVecOfPtrs(vel), Source, - m_cur_time, dummy_dt, false /*incremental*/, - true /*set_inflow_bc*/); + ApplyProjection(get_density_new_const(), + AMREX_D_DECL(GetVecOfPtrs(u_mac_tmp), GetVecOfPtrs(v_mac_tmp), + GetVecOfPtrs(w_mac_tmp)),m_cur_time,dummy_dt,false); } #ifdef AMREX_USE_EB diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index b7946cc5b..d3de63406 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -96,8 +96,13 @@ void incflo::WriteCheckPointFile() const VisMF::Write(m_leveldata[lev]->gp, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "gradp")); - VisMF::Write(m_leveldata[lev]->p_nd, - amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "p_nd")); + if (m_use_cc_proj) { + VisMF::Write(m_leveldata[lev]->p_nd, + amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "p_nd")); + } else { + VisMF::Write(m_leveldata[lev]->p_cc, + amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "p_cc")); + } } #ifdef INCFLO_USE_PARTICLES @@ -226,8 +231,13 @@ void incflo::ReadCheckpointFile() VisMF::Read(m_leveldata[lev]->gp, amrex::MultiFabFileFullPrefix(lev, m_restart_file, level_prefix, "gradp")); - VisMF::Read(m_leveldata[lev]->p_nd, - amrex::MultiFabFileFullPrefix(lev, m_restart_file, level_prefix, "p_nd")); + if (m_use_cc_proj) { + VisMF::Read(m_leveldata[lev]->p_nd, + amrex::MultiFabFileFullPrefix(lev, m_restart_file, level_prefix, "p_nd")); + } else { + VisMF::Read(m_leveldata[lev]->p_cc, + amrex::MultiFabFileFullPrefix(lev, m_restart_file, level_prefix, "p_cc")); + } } #ifdef INCFLO_USE_PARTICLES @@ -368,7 +378,7 @@ void incflo::WritePlotFile() if (m_plt_tracer) ncomp += m_ntrac; // Pressure - if(m_plt_p_nd) ++ncomp; + if(m_plt_p) ++ncomp; // MAC phi if(m_plt_macphi) ncomp += 1; @@ -488,16 +498,24 @@ void incflo::WritePlotFile() } icomp += m_ntrac; } - if (m_plt_p_nd) { - for (int lev = 0; lev <= finest_level; ++lev) - amrex::average_node_to_cellcenter(mf[lev], icomp, m_leveldata[lev]->p_nd, 0, 1); - pltscaVarsName.push_back("p_nd"); + if (m_plt_p) { + if (m_use_cc_proj) { + for (int lev = 0; lev <= finest_level; ++lev) { + MultiFab::Copy(mf[lev], m_leveldata[lev]->p_cc, 0, icomp, 1, 0); + } + } else { + for (int lev = 0; lev <= finest_level; ++lev) { + amrex::average_node_to_cellcenter(mf[lev], icomp, m_leveldata[lev]->p_nd, 0, 1); + } + } + pltscaVarsName.push_back("p"); ++icomp; } if (m_plt_macphi) { - for (int lev = 0; lev <= finest_level; ++lev) + for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy(mf[lev], m_leveldata[lev]->mac_phi, 0, icomp, 1, 0); + } pltscaVarsName.push_back("mac_phi"); ++icomp; }