From 9e82cf4c5950486f1ecab58433eb880b30ee1080 Mon Sep 17 00:00:00 2001 From: cgilet Date: Wed, 21 Feb 2024 21:27:27 -0500 Subject: [PATCH] More clean up. Improve a few comments. --- src/boundary_conditions/incflo_set_bcs.cpp | 149 +++++------------- ...ncflo_compute_MAC_projected_velocities.cpp | 68 ++------ .../incflo_compute_advection_term.cpp | 1 - src/diffusion/DiffusionScalarOp.cpp | 16 +- src/incflo.H | 3 +- src/prob/prob_bc.cpp | 4 +- .../incflo_apply_nodal_projection.cpp | 44 ------ 7 files changed, 59 insertions(+), 226 deletions(-) diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp index 237a66a4..78b4d222 100644 --- a/src/boundary_conditions/incflo_set_bcs.cpp +++ b/src/boundary_conditions/incflo_set_bcs.cpp @@ -6,42 +6,6 @@ using namespace amrex; -// iMultiFab -// incflo::make_ccBC_mask(int lev) -// { -// // BC info stored in ghost cell. -// iMultiFab new_mask(grids[lev], dmap[lev], 1, 1); -// // No expectation in MLMG for there to be valid data where there is no Robin BC -// //new_mask = 0; - -// Geometry const& gm = Geom(lev); -// Box const& domain = gm.Domain(); -// for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { -// Orientation olo(dir,Orientation::low); -// Orientation ohi(dir,Orientation::high); -// if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) { -// Box dlo = (m_bc_type[olo] == BC::mixed) ? adjCellLo(domain,dir) : Box(); -// Box dhi = (m_bc_type[ohi] == BC::mixed) ? adjCellHi(domain,dir) : Box(); -// #ifdef _OPENMP -// #pragma omp parallel if (Gpu::notInLaunchRegion()) -// #endif -// for (MFIter mfi(new_mask); mfi.isValid(); ++mfi) { -// Box blo = mfi.growntilebox() & dlo; -// Box bhi = mfi.growntilebox() & dhi; -// Array4 const& mask_arr = new_mask.array(mfi); -// if (blo.ok()) { -// prob_set_BC_MF(olo, blo, mask_arr, lev); -// } -// if (bhi.ok()) { -// prob_set_BC_MF(ohi, bhi, mask_arr, lev); -// } -// } -// } -// } - -// return new_mask; -// } - iMultiFab incflo::make_nodalBC_mask(int lev) { @@ -81,6 +45,7 @@ incflo::make_nodalBC_mask(int lev) return new_mask; } +// This function makes a position dependent BC MF for the advection routines std::unique_ptr incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector& bcs, std::string field) @@ -133,38 +98,38 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector& bcs, return BC_MF; } -Vector> -incflo::make_diffusion_robinBC_MFs(int lev, MultiFab& phi) +// TODO There's probably a way to combine this with above function rather than repeating +// much of the same code... +Vector +incflo::make_robinBC_MFs(int lev, MultiFab* phi) { int nghost = 1; - Vector> robin(3); + Vector robin(3); for (int n = 0; n < 3; n++) { - robin[n].reset(new MultiFab(grids[lev], dmap[lev], 1, 1)); + robin[n].define(grids[lev], dmap[lev], 1, nghost); } - MultiFab& robin_a = *robin[0]; - MultiFab& robin_b = *robin[1]; - MultiFab& robin_f = *robin[2]; -// fixme - for vis, init robin MFs, although i don't know that this is generally needed - robin_a = 0; - robin_b = 0; - robin_f = 0; + MultiFab& robin_a = robin[0]; + MultiFab& robin_b = robin[1]; + MultiFab& robin_f = robin[2]; - // only ghost cells of robin BC arrays are used in MLMG - // bc in ghost cells that are outside the domain. + // define and fill Robin BC a, b, and f + // BC values are stored in the ghost cells, although the BC is considered on face. + // Only the ghost cells of robin BC arrays are used in MLMG. Box const& domain = Geom(lev).Domain(); for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { Orientation olo(dir,Orientation::low); Orientation ohi(dir,Orientation::high); + // Only need to fill Robin BC sides, MLMG will check for Robin BC first if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) { Box dlo = (m_bc_type[olo] == BC::mixed) ? adjCellLo(domain,dir) : Box(); Box dhi = (m_bc_type[ohi] == BC::mixed) ? adjCellHi(domain,dir) : Box(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - // FIXME - don't think we want tiling here... + // FIXME - do we want tiling here... for (MFIter mfi(robin_a,TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& gbx = amrex::grow(mfi.validbox(),nghost); Box blo = gbx & dlo; @@ -172,18 +137,29 @@ incflo::make_diffusion_robinBC_MFs(int lev, MultiFab& phi) Array4 const& a_arr = robin_a.array(mfi); Array4 const& b_arr = robin_b.array(mfi); Array4 const& f_arr = robin_f.array(mfi); - Array4 const& bcv = phi.const_array(mfi); - // Robin BC: a u + b du/dn = f -- inflow, Dirichlet a=1, b=0, f=bcv - // -- outflow, Neumann a=0, b=1, f=0 - // Only need to fill Robin BC sides, MLMG will check for Robin BC first - if (blo.ok()) { - // fixme -- here i think we need to go to physbcfunct here -// phi has already had boundry filled, so can just use that here - prob_set_diffusion_robinBCs(olo, blo, a_arr, b_arr, f_arr, bcv, lev); - } - if (bhi.ok()) { - prob_set_diffusion_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, bcv, lev); + if (phi) { + // phi has filled boundry, so can use its ghost cells for the BC values + Array4 const& bcv = phi->const_array(mfi); + // Diffusion + // Robin BC: a u + b du/dn = f -- inflow, Dirichlet a=1, b=0, f=bcv + // -- outflow, Neumann a=0, b=1, f=0 + if (blo.ok()) { + prob_set_diffusion_robinBCs(olo, blo, a_arr, b_arr, f_arr, bcv, lev); + } + if (bhi.ok()) { + prob_set_diffusion_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, bcv, lev); + } + } else { + // MAC + // Robin BC: a u + b du/dn = f -- inflow, Neumann a=0, b=1, f=0 + // -- outflow, Dirichlet a=1, b=0, f=0 + if (blo.ok()) { + prob_set_MAC_robinBCs(olo, blo, a_arr, b_arr, f_arr, lev); + } + if (bhi.ok()) { + prob_set_MAC_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, lev); + } } } } @@ -192,57 +168,6 @@ incflo::make_diffusion_robinBC_MFs(int lev, MultiFab& phi) return robin; } -// void -// incflo::make_nodalBC_mask(int lev, -// const BoxArray& ba, // do we really need to pass these, use member vars? - i think we need to pass to be safe when calling from RemakeLevel... -// const DistributionMapping& dm) -// { -// // FIXME do we really want this check here? -// bool has_mixedBC = false; -// for (OrientationIter oit; oit; ++oit) { -// if (m_bc_type[oit()] == BC::mixed ) -// { -// has_mixedBC = true; -// break; -// } -// } - -// if (!has_mixedBC) { return; } - - -// // MLNodeLap does not require any ghost cells... -// std::unique_ptr new_mask(new iMultiFab(amrex::convert(ba,IntVect::TheNodeVector()), -// dm, 1, 0)); -// *new_mask = 1; - -// Geometry const& gm = Geom(lev); -// Box const& domain = gm.Domain(); -// for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { -// Orientation olo(dir,Orientation::low); -// Orientation ohi(dir,Orientation::high); -// if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) { -// Box dlo = (m_bc_type[olo] == BC::mixed) ? surroundingNodes(bdryLo(domain,dir)) : Box(); -// Box dhi = (m_bc_type[ohi] == BC::mixed) ? surroundingNodes(bdryHi(domain,dir)) : Box(); -// #ifdef _OPENMP -// #pragma omp parallel if (Gpu::notInLaunchRegion()) -// #endif -// for (MFIter mfi(*new_mask); mfi.isValid(); ++mfi) { -// Box blo = mfi.validbox() & dlo; -// Box bhi = mfi.validbox() & dhi; -// Array4 const& mask_arr = new_mask->array(mfi); -// if (blo.ok()) { -// prob_set_BC_MF(olo, blo, mask_arr, lev); -// } -// if (bhi.ok()) { -// prob_set_BC_MF(ohi, bhi, mask_arr, lev); -// } -// } -// } -// } - -// m_BC_MF[lev] = std::move(new_mask); -// } - #ifdef AMREX_USE_EB void incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost) diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp index fc9c74f9..63804ec7 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -91,7 +91,6 @@ incflo::compute_MAC_projected_velocities ( // // Initialize (or redefine the beta in) the MacProjector // - Print()<<"Mac proj needs init? "<needInitialization()<needInitialization()) { LPInfo lp_info; @@ -112,61 +111,12 @@ incflo::compute_MAC_projected_velocities ( } macproj->setDomainBC(get_mac_projection_bc(Orientation::low), get_mac_projection_bc(Orientation::high)); - // Not sure exactly where this goes. Does it only need initialization? set every time??? - // initProj only take as many levels as 1/rho, so if a level is added, would have to - // re-init proj anyway. - // robin_a, etc needs to exist to finest_level - // Values need to be in the ghost cells, although the BC is considered on face. if ( m_has_mixedBC ) { - int nghost = 1; - for (int lev = 0; lev <= finest_level; ++lev) { - // define and fill Robin BC a, b, and f - MultiFab robin_a (grids[lev],dmap[lev],1,nghost,MFInfo(),Factory(lev)); - MultiFab robin_b (grids[lev],dmap[lev],1,nghost,MFInfo(),Factory(lev)); - MultiFab robin_f (grids[lev],dmap[lev],1,nghost,MFInfo(),Factory(lev)); - - // fixme? - for vis, init robin MFs, although i don't know that this is generally needed - robin_a = 0; - robin_b = 0; - robin_f = 0; - - // only ghost cells of robin BC arrays are used in MLMG - // bc in ghost cells that are outside the domain. - Box const& domain = Geom(lev).Domain(); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - Orientation olo(dir,Orientation::low); - Orientation ohi(dir,Orientation::high); - - if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) { - Box dlo = (m_bc_type[olo] == BC::mixed) ? adjCellLo(domain,dir) : Box(); - Box dhi = (m_bc_type[ohi] == BC::mixed) ? adjCellHi(domain,dir) : Box(); -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - // FIXME - do we want tiling here... - for (MFIter mfi(robin_a,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - Box const& gbx = amrex::grow(mfi.validbox(),nghost); - Box blo = gbx & dlo; - Box bhi = gbx & dhi; - Array4 const& a_arr = robin_a.array(mfi); - Array4 const& b_arr = robin_b.array(mfi); - Array4 const& f_arr = robin_f.array(mfi); - - // Robin BC: a u + b du/dn = f -- inflow, Neumann a=0, b=1, f=0 - // -- outflow, Dirichlet a=1, b=0, f=0 - // Only need to fill Robin BC sides, MLMG will check for Robin BC first - if (blo.ok()) { - prob_set_MAC_robinBCs(olo, blo, a_arr, b_arr, f_arr, lev); - } - if (bhi.ok()) { - prob_set_MAC_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, lev); - } - } - } - } - macproj->setLevelBC(lev, nullptr, &robin_a, &robin_b, &robin_f); + auto const robin = make_robinBC_MFs(lev); + macproj->setLevelBC(lev, nullptr, + &robin[0], &robin[1], &robin[2]); } } } else { @@ -197,12 +147,16 @@ incflo::compute_MAC_projected_velocities ( l_advection_type = "Godunov"; } -// FIXME? - need to create the BC MF here. This is a iMF vs a Real MF, so can't use - // the robin BC info from MAC. This is because we need this to hold the BCType, - // which are enums, so ints - // std::unique_ptr BC_MF; if (m_has_mixedBC) { + // Create a MF to hold the BCType info. Note that this is different than the + // bcs for the MAC projection because the MAC operates on phi, this is velocity. + // + // TODO? Could consider creating an incflo member variable to save the BC_MF + // amrex::Vector > m_BC_MF; + // Could stack components as vel, density, tracer, and then could use for + // scalars' advective step as well. But not sure it really matters one way or + // the other. BC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity"); } diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 0f49c0c9..8915ba17 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -372,7 +372,6 @@ incflo::compute_convective_term (Vector const& conv_u, is_velocity, fluxes_are_area_weighted, m_advection_type, PPM::default_limiter, velbc_arr); - // ************************************************************************ // Density // ************************************************************************ diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp index 4e3bfee4..d8139a28 100644 --- a/src/diffusion/DiffusionScalarOp.cpp +++ b/src/diffusion/DiffusionScalarOp.cpp @@ -221,10 +221,10 @@ DiffusionScalarOp::diffuse_scalar (Vector const& tracer, #ifdef AMREX_USE_EB if (m_eb_scal_solve_op) { if ( m_incflo->m_has_mixedBC ) { - auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, phi[lev]); + auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]); m_eb_scal_solve_op->setLevelBC(lev, &phi[lev], - robin[0].get(), robin[1].get(), robin[2].get()); + &robin[0], &robin[1], &robin[2]); } else { m_eb_scal_solve_op->setLevelBC(lev, &phi[lev]); @@ -376,10 +376,10 @@ DiffusionScalarOp::diffuse_vel_components (Vector const& vel, // m_eb_vel_solve_op->setPhiOnCentroid(); if ( m_incflo->m_has_mixedBC ) { - auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, phi[lev]); + auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]); m_eb_vel_solve_op->setLevelBC(lev, &phi[lev], - robin[0].get(), robin[1].get(), robin[2].get()); + &robin[0], &robin[1], &robin[2]); } else { m_eb_vel_solve_op->setLevelBC(lev, &phi[lev]); @@ -483,10 +483,10 @@ void DiffusionScalarOp::compute_laps (Vector const& a_laps, m_eb_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid); if ( m_incflo->m_has_mixedBC ) { - auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, scalar_comp[lev]); + auto const robin = m_incflo->make_robinBC_MFs(lev, &scalar_comp[lev]); m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev], - robin[0].get(), robin[1].get(), robin[2].get()); + &robin[0], &robin[1], &robin[2]); } else { m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev]); @@ -605,10 +605,10 @@ void DiffusionScalarOp::compute_divtau (Vector const& a_divtau, vel_single.emplace_back( vel[lev],amrex::make_alias,comp,1); if ( m_incflo->m_has_mixedBC ) { - auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, vel_single[lev]); + auto const robin = m_incflo->make_robinBC_MFs(lev, &vel_single[lev]); m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev], - robin[0].get(), robin[1].get(), robin[2].get()); + &robin[0], &robin[1], &robin[2]); } else { m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev]); diff --git a/src/incflo.H b/src/incflo.H index b0553e03..3ce060b1 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -114,8 +114,7 @@ public: amrex::Gpu::DeviceVector& bcs, std::string field); amrex::iMultiFab make_nodalBC_mask (int lev); - amrex::Vector> make_diffusion_robinBC_MFs(int lev, - amrex::MultiFab& phi); + amrex::Vector make_robinBC_MFs(int lev, amrex::MultiFab* phi = nullptr); // void make_ccBC_mask (int lev, const amrex::BoxArray& ba, // const amrex::DistributionMapping& dm); // void make_nodalBC_mask (int lev, const amrex::BoxArray& ba, diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp index 1aa89989..48d3bf9d 100644 --- a/src/prob/prob_bc.cpp +++ b/src/prob/prob_bc.cpp @@ -88,13 +88,13 @@ void incflo::prob_set_BC_MF (Orientation ori, Box const& bx, } } -// For MAC, the BC is on phi, so always 0 or 1 void incflo::prob_set_MAC_robinBCs (Orientation ori, Box const& bx, Array4 const& robin_a, Array4 const& robin_b, Array4 const& robin_f, int lev) { + // For MAC, the BC is on phi, so always 0 or 1 // Robin BC: a u + b du/dn = f -- inflow, Neumann a=0, b=1, f=0 // -- outflow, Dirichlet a=1, b=0, f=0 @@ -186,7 +186,6 @@ void incflo::prob_set_MAC_robinBCs (Orientation ori, Box const& bx, } } -// For diffusion, we also pass in the dirichlet bc void incflo::prob_set_diffusion_robinBCs (Orientation ori, Box const& bx, Array4 const& robin_a, Array4 const& robin_b, @@ -194,6 +193,7 @@ void incflo::prob_set_diffusion_robinBCs (Orientation ori, Box const& bx, Array4 const& bcval, int lev) { + // For diffusion, we also pass in the dirichlet bc (bcval) // Robin BC: a u + b du/dn = f -- inflow, Dirichlet a=1, b=0, f=bcval // -- outflow, Neumann a=0, b=1, f=0 diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index 1ea2e5d3..8847a275 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -161,50 +161,6 @@ void incflo::ApplyNodalProjection (Vector density, for(int lev = 0; lev <= finest_level; ++lev) { const auto mask = make_nodalBC_mask(lev); -// // MLNodeLap does not require any ghost cells... -// std::unique_ptr nodal_mask(new iMultiFab(amrex::convert(ba,IntVect::TheNodeVector()), -// dm, 1, 0)); -// *nodal_mask = 1; - -// Box const& domain = Geom(lev).Domain(); -// for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { -// Orientation olo(dir,Orientation::low); -// Orientation ohi(dir,Orientation::high); -// if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) { -// Box dlo = (m_bc_type[olo] == BC::mixed) ? surroundingNodes(bdryLo(domain,dir)) : Box(); -// Box dhi = (m_bc_type[ohi] == BC::mixed) ? surroundingNodes(bdryHi(domain,dir)) : Box(); -// #ifdef _OPENMP -// #pragma omp parallel if (Gpu::notInLaunchRegion()) -// #endif -// for (MFIter mfi(*nodal_mask); mfi.isValid(); ++mfi) { -// Box blo = mfi.validbox() & dlo; -// Box bhi = mfi.validbox() & dhi; -// Array4 const& nodal_arr = nodal_mask->array(mfi); -// Array4 const& bcs_arr = m_mixedBC_mask[lev]->const_array(mfi); -// if (blo.ok()) { -// // nodal_mask is the same as the mixedBC_mask, except with vals on -// // on BC face rather than in cell-centered ghosts -// // So, lo side i_nd->i_cc-1, hi side i_nd=i_cc - -// //FIXME, theres a corner cell issue here :( -// // maybe answer is to just go to make_mask fn in prob... -// Dim3 shift = IntVect::TheDimensionVector(dir).dim3(); - -// amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -// { -// nodal_arr(i,j,k) = bcs_arr(i-shift.x, j-shift.y, k-shift.z); -// }); -// } -// if (bhi.ok()) { -// amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -// { -// nodal_arr(i,j,k) = bcs_arr(i, j, k); -// }); -// } -// } -// } -// } - nodal_projector->getLinOp().setOversetMask(lev, mask); } }