diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp index 303ae4a0..b201ae80 100644 --- a/src/boundary_conditions/incflo_set_bcs.cpp +++ b/src/boundary_conditions/incflo_set_bcs.cpp @@ -6,41 +6,41 @@ 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; +// 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); - } - } - } - } +// 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; -} +// return new_mask; +// } iMultiFab incflo::make_nodalBC_mask(int lev) @@ -49,6 +49,9 @@ incflo::make_nodalBC_mask(int lev) iMultiFab new_mask(amrex::convert(grids[lev], IntVect::TheNodeVector()), dmap[lev], 1, 0); // Overset mask needs to be defined everywhere new_mask = 1; + // NodalProj expects outflow regions to be indcatied with a 0. This is treated as a boolean + // mask, so any other regions can have any other int value + int outflow = 0; Geometry const& gm = Geom(lev); Box const& domain = gm.Domain(); @@ -66,10 +69,10 @@ incflo::make_nodalBC_mask(int lev) 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); + prob_set_BC_MF(olo, blo, mask_arr, lev, outflow); } if (bhi.ok()) { - prob_set_BC_MF(ohi, bhi, mask_arr, lev); + prob_set_BC_MF(ohi, bhi, mask_arr, lev, outflow); } } } @@ -78,54 +81,56 @@ incflo::make_nodalBC_mask(int lev) return new_mask; } -// void -// incflo::make_ccBC_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) -// { -// bool has_mixedBC = false; -// for (OrientationIter oit; oit; ++oit) { -// if (m_bc_type[oit()] == BC::mixed ) -// { -// has_mixedBC = true; -// break; -// } -// } - -// if (!has_mixedBC) { return; } - - -// // BC info stored in ghost cell. -// std::unique_ptr new_mask(new iMultiFab(ba, dm, 1, 1)); -// *new_mask = 0; +std::unique_ptr +incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector bcs) //, int scomp, int ncomp) +{ + int ncomp = bcs.size(); + // BC info stored in ghost cells to avoid any ambiguity + std::unique_ptr BC_MF(new iMultiFab(grids[lev], dmap[lev], ncomp, 1)); + // initialize to 0 == BCType::int_dir, not sure this is needed really + *BC_MF = 0; + // For advection, we use FOEXTRAP for outflow regions per incflo::init_bcs() + int outflow = BCType::foextrap; + + Box const& domain = geom[lev].Domain(); + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + Orientation olo(dir,Orientation::low); + Orientation ohi(dir,Orientation::high); -// 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); -// } -// } -// } -// } + Box dlo = adjCellLo(domain,dir); + Box dhi = adjCellHi(domain,dir); +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*BC_MF); mfi.isValid(); ++mfi) { + Box blo = mfi.growntilebox() & dlo; + Box bhi = mfi.growntilebox() & dhi; + Array4 const& mask_arr = BC_MF->array(mfi); + if (m_bc_type[olo] == BC::mixed) { + prob_set_BC_MF(olo, blo, mask_arr, lev, outflow); + } else { + amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int n = 0; n < ncomp; n++) { + mask_arr(i,j,k,n) = bcs[n].lo(dir); + } + }); + } + if (m_bc_type[ohi] == BC::mixed) { + prob_set_BC_MF(ohi, bhi, mask_arr, lev, outflow); + } else { + amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + for (int n = 0; n < ncomp; n++) { + mask_arr(i,j,k,n) = bcs[n].hi(dir); + } + }); + } + } + } -// m_BC_MF[lev] = std::move(new_mask); -// } + return BC_MF; +} // void // incflo::make_nodalBC_mask(int lev, diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp index 42ce6881..49c231ea 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -233,56 +233,10 @@ incflo::compute_MAC_projected_velocities ( // full position dependent BCs, would need the MF to hold the BCType, which are enums, // so ints... // -// if (m_mixedBC_mask[lev]) { - -// // Create MF to hold position-dependent BC info in ghost cells: -// // Low x-face is represented by indices ((i_lo-1, j_lo, k_lo)(i_lo-1, j_hi, k_hi)) -// // High x-face is represented by indices ((i_hi+1, j_lo, k_lo)(i_hi+1, j_hi, k_hi)) -// // Interior values don't get used -// iMultiFab bcs (grids[lev],dmap[lev],1,1,MFInfo(),Factory(lev)); - -// for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { -// Orientation olo(dir,Orientation::low); -// Orientation ohi(dir,Orientation::high); - -// Box dlo = amrex::adjCellLo(domain,dir,nghost); -// Box dhi = amrex::adjCellHi(domain,dir,nghost); -// #ifdef _OPENMP -// #pragma omp parallel if (Gpu::notInLaunchRegion()) -// #endif -// for (MFIter mfi(bcs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { -// Box const& gbx = mfi.growntilebox(); -// Box blo = gbx & dlo; -// Box bhi = gbx & dhi; -// Array4 const& bcs_arr = bcs.array(mfi); -// Array4 const& mask = m_mixedBC_mask[lev]->const_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 -// if (blo.ok()) { -// // robin_b is the same as the mixedBC_mask, except with vals in -// // cell-centered ghosts rather than on BC face. -// // So, lo side i_cc->i_nd+1, hi side i_cc=i_nd -// Dim3 shift = IntVect::TheDimensionVector(dir).dim3(); - -// amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -// { -// b_arr(i,j,k) = mask(i+shift.x, j+shift.y, k+shift.z); -// //robin_a is the "opposite" of b; 0->1 and vice versa -// a_arr(i,j,k) = (mask(i+shift.x, j+shift.y, k+shift.z) == 1) ? 0. : 1.; -// f_arr(i,j,k) = 0.; -// }); -// } -// if (bhi.ok()) { -// amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -// { -// b_arr(i,j,k) = mask(i,j,k); -// a_arr(i,j,k) = (mask(i,j,k) == 1) ? 0. : 1.; -// f_arr(i,j,k) = 0.; -// }); -// } -// } -// } + std::unique_ptr BC_MF; + if (m_has_mixedBC) { + BC_MF = make_BC_MF(lev, m_bcrec_velocity_d); + } // Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac} // returned from this call are on face CENTROIDS @@ -295,101 +249,13 @@ incflo::compute_MAC_projected_velocities ( m_eb_flow.enabled ? get_velocity_eb()[lev] : nullptr, #endif m_godunov_ppm, m_godunov_use_forces_in_trans, - l_advection_type); - } - - //fixme - fix up umac to confirm projection working as expected... - Box const& domain_box = geom[0].Domain(); - int direction = 1; - int half_num_cells = domain_box.length(direction) / 2; - for (MFIter mfi(*u_mac[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) { - Box const& bx = mfi.validbox(); - Array4 const& um = u_mac[0]->array(mfi); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (i == domain_box.smallEnd(0)) - { - if (j > half_num_cells) { - // Here we take the inflow BC specified in inputs file - um(i,j,k) = 1.; - } - } - else if (i == domain_box.bigEnd(0)+1) - { - if (j <= half_num_cells) { - // Here we take minus the inflow BC specified in inputs file - um(i,j,k) = -1.; - } - } - }); + l_advection_type, PPM::default_limiter, BC_MF.get()); } -// if (m_mixedBC_mask[0]) { -// // Fix up the mixedBC faces which went through advection with FOEXTRAP (outflow) BCs -// // Here we overwrite the inflow portion with the Dirichlet BC -// // NOTE - there's a subtle difference in this versus what AMReX-Hydro does. See comment -// // in Utils/hydro_bcs_K.H . If weird gradients start arising at the inflow, we might -// // want to revisit this... -// // vel has filled ghost cells, this is needed for ExtrapVel -// for (int lev = 0; lev <= finest_level; ++lev) -// { -// int nghost = 1; -// Box const& domain = Geom(lev).growPeriodicDomain(nghost); -// for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { -// Orientation olo(dir,Orientation::low); -// Orientation ohi(dir,Orientation::high); - -// Box dlo = amrex::adjCellLo(domain,dir,nghost); -// Box dhi = amrex::adjCellHi(domain,dir,nghost); -// #ifdef _OPENMP -// #pragma omp parallel if (Gpu::notInLaunchRegion()) -// #endif -// for (MFIter mfi(*vel[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { -// Box const& gbx = amrex::grow(mfi.validbox(),1); -// Box blo = gbx & dlo; -// Box bhi = gbx & dhi; - -// MultiFab* mac; -// if (dir == 0) { -// mac = u_mac[lev]; } -// else if (dir == 1) { -// mac = v_mac[lev]; } -// #if AMREX_SPACEDIM > 2 -// else { -// mac = w_mac[lev]; -// } -// #endif -// Array4 const& mac_arr = mac->array(mfi); -// Array4 const& vel_arr = vel[lev]->array(mfi, dir); -// Array4 const& mask = m_mixedBC_mask[lev]->const_array(mfi); - -// if (blo.ok()) { -// Dim3 shift = IntVect::TheDimensionVector(dir).dim3(); - -// amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -// { -// if (mask(i+shift.x, j+shift.y, k+shift.z) == 1) // inflow -// { -// mac_arr(i+shift.x, j+shift.y, k+shift.z) = vel_arr(i,j,k); -// } -// }); -// } -// if (bhi.ok()) { -// amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -// { -// if (mask(i,j,k) == 1) { -// mac_arr(i,j,k) = vel_arr(i,j,k); -// } -// }); -// } -// } -// } -// } -// } -//fixme +//fixmeamrex:: VisMF::Write(*u_mac[0],"umac"); VisMF::Write(*v_mac[0],"vmac"); - +// Vector > mac_vec(finest_level+1); for (int lev=0; lev <= finest_level; ++lev) { diff --git a/src/incflo.H b/src/incflo.H index 83f504ea..4f3dd271 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -110,7 +110,7 @@ public: // /////////////////////////////////////////////////////////////////////////// - amrex::iMultiFab make_ccBC_mask (int lev); + std::unique_ptr make_BC_MF (int lev, amrex::Gpu::DeviceVector bcs); amrex::iMultiFab make_nodalBC_mask (int lev); // void make_ccBC_mask (int lev, const amrex::BoxArray& ba, // const amrex::DistributionMapping& dm); @@ -218,7 +218,8 @@ public: void prob_init_fluid (int lev); void prob_set_BC_MF (amrex::Orientation ori, amrex::Box const& bx, - amrex::Array4 const& mask, int lev); + amrex::Array4 const& mask, int lev, + int outflow_val); void prob_set_MAC_robinBCs (amrex::Orientation ori, amrex::Box const& bx, amrex::Array4 const& robin_a, amrex::Array4 const& robin_b, diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp index 18c8818b..8d975c01 100644 --- a/src/prob/prob_bc.cpp +++ b/src/prob/prob_bc.cpp @@ -2,8 +2,12 @@ using namespace amrex; +// This function is used to prepare both the nodal projection and advection for +// mixed BCs (i.e. position dependent BCs). The necessarily value to indicate an +// outflow BC differs between advection and the NodalProj, so we pass it in void incflo::prob_set_BC_MF (Orientation ori, Box const& bx, - Array4 const& mask, int lev) + Array4 const& mask, int lev, + int outflow_val) { if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype) { @@ -23,21 +27,21 @@ void incflo::prob_set_BC_MF (Orientation ori, Box const& bx, { if (direction == 0) { if (i <= half_num_cells) { - mask(i,j,k,0) = 0; // outflow on bottom + mask(i,j,k,0) = outflow_val; // outflow on bottom } else { mask(i,j,k,0) = BCType::ext_dir; // inflow on top } } else if (direction == 1) { if (j <= half_num_cells) { - mask(i,j,k,0) = 0; + mask(i,j,k,0) = outflow_val; } else { mask(i,j,k,0) = BCType::ext_dir; } } else if (direction == 2) { if (k <= half_num_cells) { - mask(i,j,k,0) = 0; + mask(i,j,k,0) = outflow_val; } else { mask(i,j,k,0) = BCType::ext_dir; } @@ -48,21 +52,21 @@ void incflo::prob_set_BC_MF (Orientation ori, Box const& bx, { if (direction == 0) { if (i > half_num_cells) { - mask(i,j,k,0) = 0; // outflow on top + mask(i,j,k,0) = outflow_val; // outflow on top } else { mask(i,j,k,0) = BCType::ext_dir; // inflow on bottom } } else if (direction == 1) { if (j > half_num_cells) { - mask(i,j,k,0) = 0; + mask(i,j,k,0) = outflow_val; } else { mask(i,j,k,0) = BCType::ext_dir; } } else if (direction == 2) { if (k > half_num_cells) { - mask(i,j,k,0) = 0; + mask(i,j,k,0) = outflow_val; } else { mask(i,j,k,0) = BCType::ext_dir; }