diff --git a/src/boundary_conditions/boundary_conditions.cpp b/src/boundary_conditions/boundary_conditions.cpp index d34b16d6..1781be73 100644 --- a/src/boundary_conditions/boundary_conditions.cpp +++ b/src/boundary_conditions/boundary_conditions.cpp @@ -94,6 +94,7 @@ void incflo::init_bcs () else if (bc_type == "mixed" ) { amrex::Print() << bcid << " set to mixed inflow outflow.\n"; + m_has_mixedBC = true; #ifdef AMREX_USE_EB ParmParse ipp("incflo"); diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp index c5d734d4..303ae4a0 100644 --- a/src/boundary_conditions/incflo_set_bcs.cpp +++ b/src/boundary_conditions/incflo_set_bcs.cpp @@ -6,27 +6,49 @@ using namespace amrex; -void -incflo::make_mixedBC_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) +iMultiFab +incflo::make_ccBC_mask(int lev) { - bool has_mixedBC = false; - for (OrientationIter oit; oit; ++oit) { - if (m_bc_type[oit()] == BC::mixed ) - { - has_mixedBC = true; - break; + // 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); + } + } } } - if (!has_mixedBC) { return; } - + return new_mask; +} +iMultiFab +incflo::make_nodalBC_mask(int lev) +{ // 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; + iMultiFab new_mask(amrex::convert(grids[lev], IntVect::TheNodeVector()), dmap[lev], 1, 0); + // Overset mask needs to be defined everywhere + new_mask = 1; Geometry const& gm = Geom(lev); Box const& domain = gm.Domain(); @@ -39,23 +61,123 @@ incflo::make_mixedBC_mask(int lev, #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*new_mask); mfi.isValid(); ++mfi) { + 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); + Array4 const& mask_arr = new_mask.array(mfi); if (blo.ok()) { - prob_set_mixedBC_mask(olo, blo, mask_arr, lev); + prob_set_BC_MF(olo, blo, mask_arr, lev); } if (bhi.ok()) { - prob_set_mixedBC_mask(ohi, bhi, mask_arr, lev); + prob_set_BC_MF(ohi, bhi, mask_arr, lev); } } } } - m_mixedBC_mask[lev] = std::move(new_mask); + 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; + +// 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); +// } +// } +// } +// } + +// m_BC_MF[lev] = std::move(new_mask); +// } + +// 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 41d56457..05c32e21 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -91,6 +91,7 @@ incflo::compute_MAC_projected_velocities ( // // Initialize (or redefine the beta in) the MacProjector // + Print()<<"Mac proj needs init? "<needInitialization()<needInitialization()) { LPInfo lp_info; @@ -116,7 +117,7 @@ incflo::compute_MAC_projected_velocities ( // re-init proj anyway. // robin_a, etc needs to exist to finest_level and that's it // Values need to be in the ghost cells, although the BC is considered on face. - if ( m_mixedBC_mask[0] ) { + if ( m_has_mixedBC ) { int nghost = 1; for (int lev = 0; lev <= finest_level; ++lev) @@ -126,6 +127,11 @@ incflo::compute_MAC_projected_velocities ( 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; + // I don't believe interior values get used, // only ghost cells of robin BC arrays are used in MLMG // bc in ghost cells that are outside the domain. @@ -138,47 +144,59 @@ incflo::compute_MAC_projected_velocities ( Orientation olo(dir,Orientation::low); Orientation ohi(dir,Orientation::high); - Box dlo = amrex::adjCellLo(domain,dir,nghost); - Box dhi = amrex::adjCellHi(domain,dir,nghost); + 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(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); - 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.; - }); + // FIXME - don't think 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); + //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 +// Only need to fill Robin BC sides, MLMG will check for Robin BC first + if (blo.ok()) { + // robin_b is the same as the nodalBC_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(); + + prob_set_MAC_robinBCs(olo, blo, a_arr, b_arr, f_arr, lev); + + // 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()) { + + prob_set_MAC_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, lev); + + // 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.; + // }); + } } } } + VisMF::Write(robin_a, "ra"); + VisMF::Write(robin_b, "rb"); + VisMF::Write(robin_f, "rf"); macproj->setLevelBC(lev, nullptr, &robin_a, &robin_b, &robin_f); } } @@ -210,6 +228,62 @@ 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 if we want this to be extensible to + // 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.; +// }); +// } +// } +// } + // Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac} // returned from this call are on face CENTROIDS HydroUtils::ExtrapVelToFaces(*vel[lev], *vel_forces[lev], @@ -224,69 +298,98 @@ incflo::compute_MAC_projected_velocities ( l_advection_type); } - 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) + //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 { - 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); - } - }); - } + 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.; + } + } + }); } +// 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 + 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) { @@ -325,7 +428,9 @@ incflo::compute_MAC_projected_velocities ( } else { macproj->project(m_mac_mg_rtol,m_mac_mg_atol); } - +//fixme + VisMF::Write(*u_mac[0],"umacp"); + VisMF::Write(*v_mac[0],"vmacp"); // Note that the macproj->project call above ensures that the MAC velocities are averaged down -- // we don't need to do that again here } diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp index 1e52b2f6..54bf310e 100644 --- a/src/convection/incflo_compute_advection_term.cpp +++ b/src/convection/incflo_compute_advection_term.cpp @@ -443,6 +443,70 @@ incflo::compute_convective_term (Vector const& conv_u, } // mfi } // lev +// FIXME - probably want to make this a fn, since it's the same as get used in compute_MAC +// 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 making edge states +// 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); +// } +// }); +// } +// } +// } +// } +// } + // In order to enforce conservation across coarse-fine boundaries we must be sure to average down the fluxes // before we use them. Note we also need to average down the face states if we are going to do // convective differencing diff --git a/src/incflo.H b/src/incflo.H index 5bc9f969..83f504ea 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -110,8 +110,13 @@ public: // /////////////////////////////////////////////////////////////////////////// - void make_mixedBC_mask (int lev, const amrex::BoxArray& ba, - const amrex::DistributionMapping& dm); + amrex::iMultiFab make_ccBC_mask (int lev); + amrex::iMultiFab make_nodalBC_mask (int lev); + // void make_ccBC_mask (int lev, const amrex::BoxArray& ba, + // const amrex::DistributionMapping& dm); + // void make_nodalBC_mask (int lev, const amrex::BoxArray& ba, + // const amrex::DistributionMapping& dm); + #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); @@ -212,8 +217,13 @@ public: /////////////////////////////////////////////////////////////////////////// void prob_init_fluid (int lev); - void prob_set_mixedBC_mask (amrex::Orientation ori, amrex::Box const& bx, - amrex::Array4 const& mask, int lev); + void prob_set_BC_MF (amrex::Orientation ori, amrex::Box const& bx, + amrex::Array4 const& mask, int lev); + void prob_set_MAC_robinBCs (amrex::Orientation ori, amrex::Box const& bx, + amrex::Array4 const& robin_a, + amrex::Array4 const& robin_b, + amrex::Array4 const& robin_f, + int lev); #include "incflo_prob_I.H" #include "incflo_prob_usr_I.H" @@ -570,9 +580,12 @@ private: periodic, mixed, undefined }; +// FIXME - don;t think this is the way to go, just use flag and make when needed since different + // code routines need the info in different containers... // mask for mixed BCs: 1 = inflow (Neumann for NodalProj, Dirichlet for advection), // 0 = outflow (Homogeneous Dirichlet for NodalProj, FOEXTRAP for advection) - amrex::Vector > m_mixedBC_mask; + amrex::Vector > m_BC_MF; + bool m_has_mixedBC = false; amrex::GpuArray m_bc_type; amrex::GpuArray m_bc_pressure; diff --git a/src/incflo.cpp b/src/incflo.cpp index 13e43e59..70398a15 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -223,7 +223,7 @@ void incflo::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_gr if (m_restart_file.empty()) { prob_init_fluid(lev); } - make_mixedBC_mask(lev, grids[lev], dmap[lev]); + //make_mixedBC_mask(lev, grids[lev], dmap[lev]); #ifdef AMREX_USE_EB macproj = std::make_unique(Geom(0,finest_level), diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index ee57f33a..67a41c05 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -137,10 +137,15 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************* // Compute explicit viscous term // ************************************************************************************* - if (need_divtau() || use_tensor_correction) + if (need_divtau() || use_tensor_correction ) { - compute_divtau(get_divtau_old(),get_velocity_old_const(), - get_density_old_const(),GetVecOfConstPtrs(vel_eta)); + // FIXME - deal with viscosity later... + auto divtau = get_divtau_old(); + for (int lev = 0; lev <= finest_level; ++lev) { + *divtau[lev] = 0.; + } + // compute_divtau(get_divtau_old(),get_velocity_old_const(), + // get_density_old_const(),GetVecOfConstPtrs(vel_eta)); } // ************************************************************************************* diff --git a/src/incflo_regrid.cpp b/src/incflo_regrid.cpp index 2832b172..33c4a301 100644 --- a/src/incflo_regrid.cpp +++ b/src/incflo_regrid.cpp @@ -43,7 +43,7 @@ void incflo::MakeNewLevelFromCoarse (int lev, m_leveldata[lev] = std::move(new_leveldata); m_factory[lev] = std::move(new_fact); - make_mixedBC_mask(lev, ba, dm); + //make_mixedBC_mask(lev, ba, dm); m_diffusion_tensor_op.reset(); m_diffusion_scalar_op.reset(); @@ -96,7 +96,7 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba, m_leveldata[lev] = std::move(new_leveldata); m_factory[lev] = std::move(new_fact); - make_mixedBC_mask(lev, ba, dm); + //make_mixedBC_mask(lev, ba, dm); m_diffusion_tensor_op.reset(); m_diffusion_scalar_op.reset(); @@ -118,7 +118,7 @@ void incflo::ClearLevel (int lev) BL_PROFILE("incflo::ClearLevel()"); m_leveldata[lev].reset(); m_factory[lev].reset(); - m_mixedBC_mask[lev].reset(); + m_BC_MF[lev].reset(); m_diffusion_tensor_op.reset(); m_diffusion_scalar_op.reset(); macproj.reset(); diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp index cf7f82c1..f4c924b3 100644 --- a/src/prob/prob_bc.cpp +++ b/src/prob/prob_bc.cpp @@ -2,8 +2,8 @@ using namespace amrex; -void incflo::prob_set_mixedBC_mask (Orientation ori, Box const& bx, - Array4 const& mask, int lev) +void incflo::prob_set_BC_MF (Orientation ori, Box const& bx, + Array4 const& mask, int lev) { if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype) { @@ -24,16 +24,22 @@ void incflo::prob_set_mixedBC_mask (Orientation ori, Box const& bx, if (direction == 0) { if (i <= half_num_cells) { mask(i,j,k,0) = 0; // 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; + } else { + mask(i,j,k,0) = BCType::ext_dir; } } else if (direction == 2) { if (k <= half_num_cells) { mask(i,j,k,0) = 0; + } else { + mask(i,j,k,0) = BCType::ext_dir; } } }); @@ -43,16 +49,22 @@ void incflo::prob_set_mixedBC_mask (Orientation ori, Box const& bx, if (direction == 0) { if (i > half_num_cells) { mask(i,j,k,0) = 0; // 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; + } else { + mask(i,j,k,0) = BCType::ext_dir; } } else if (direction == 2) { if (k > half_num_cells) { mask(i,j,k,0) = 0; + } else { + mask(i,j,k,0) = BCType::ext_dir; } } }); @@ -60,7 +72,107 @@ void incflo::prob_set_mixedBC_mask (Orientation ori, Box const& bx, } else { - Abort("incflo::prob_set_mixedBC_mask: No masking function for probtype " + Abort("incflo::prob_set_BC_MF: No masking function for probtype " + +std::to_string(m_probtype)); + } +} + +// FIXME - need one like this for diffusion, but then it really needs to be going to the PhysBCFunct +// Because we need the Dirichelet val, right? +// For MAC, the BC is on phi, so always 0 or 1, I think +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) +{ + // 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 (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype) + { + int direction = 0; + if (1101 == m_probtype) { + direction = 1; + } + else if (1102 == m_probtype) { + direction = 2; + } + Box const& domain = geom[lev].Domain(); + int half_num_cells = domain.length(direction) / 2; + + Orientation::Side side = ori.faceDir(); + if (side == Orientation::low) { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + robin_f(i,j,k) = 0.; + + if (direction == 0) { + if (i <= half_num_cells) { + robin_a(i,j,k,0) = 1.; + robin_b(i,j,k,0) = 0; // outflow on bottom + } else { + robin_a(i,j,k,0) = 0.; + robin_b(i,j,k,0) = 1.; // inflow on top + } + } + else if (direction == 1) { + if (j <= half_num_cells) { + robin_a(i,j,k,0) = 1.; + robin_b(i,j,k,0) = 0; + } else { + robin_a(i,j,k,0) = 0.; + robin_b(i,j,k,0) = 1.; + } + } + else if (direction == 2) { + if (k <= half_num_cells) { + robin_a(i,j,k,0) = 1.; + robin_b(i,j,k,0) = 0; + } else { + robin_a(i,j,k,0) = 0.; + robin_b(i,j,k,0) = 1.; + } + } + }); + } else { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + robin_f(i,j,k,0) = 0.; + + if (direction == 0) { + if (i > half_num_cells) { + robin_a(i,j,k,0) = 1.; + robin_b(i,j,k,0) = 0; // outflow on top + } else { + robin_a(i,j,k,0) = 0.; + robin_b(i,j,k,0) = 1.; // inflow on bottom + } + } + else if (direction == 1) { + if (j > half_num_cells) { + robin_a(i,j,k,0) = 1.; + robin_b(i,j,k,0) = 0; + } else { + robin_a(i,j,k,0) = 0.; + robin_b(i,j,k,0) = 1.; + } + } + else if (direction == 2) { + if (k > half_num_cells) { + robin_a(i,j,k,0) = 1.; + robin_b(i,j,k,0) = 0; + } else { + robin_a(i,j,k,0) = 0.; + robin_b(i,j,k,0) = 1.; + } + } + }); + } + } + else + { + Abort("incflo::prob_set_MAC_robinBCs: No masking function for probtype " +std::to_string(m_probtype)); } } diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index 63a8a34b..1ea2e5d3 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -155,11 +155,57 @@ void incflo::ApplyNodalProjection (Vector density, } } - if(m_mixedBC_mask[0]) + if(m_has_mixedBC) { + // We use an overset mask to effectively apply the mixed BC for(int lev = 0; lev <= finest_level; ++lev) { - nodal_projector->getLinOp().setOversetMask(lev, *m_mixedBC_mask[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); } } diff --git a/src/setup/incflo_arrays.cpp b/src/setup/incflo_arrays.cpp index ee2cfb18..67c79ba6 100644 --- a/src/setup/incflo_arrays.cpp +++ b/src/setup/incflo_arrays.cpp @@ -60,5 +60,5 @@ void incflo::ResizeArrays () m_factory.resize(max_level+1); - m_mixedBC_mask.resize(max_level+1); + //m_BC_MF.resize(max_level+1); }