From 07f9240b16e58cad63f09e460f104ac4bf94f1d1 Mon Sep 17 00:00:00 2001 From: cgilet Date: Wed, 7 Feb 2024 15:27:41 -0500 Subject: [PATCH] WIP - the mac and nodal proj work as expected for mixed BCs (given the appropriate setup in inputs.split). advection not working (hacked to get umac to project) not in any way elegantly coded Also requires minor update in AMReX-Hydro to pass through robin BC for MAC projector --- .../boundary_conditions.cpp | 1 + src/boundary_conditions/incflo_set_bcs.cpp | 162 ++++++++-- ...ncflo_compute_MAC_projected_velocities.cpp | 293 ++++++++++++------ .../incflo_compute_advection_term.cpp | 64 ++++ src/incflo.H | 23 +- src/incflo.cpp | 2 +- src/incflo_apply_predictor.cpp | 11 +- src/incflo_regrid.cpp | 6 +- src/prob/prob_bc.cpp | 118 ++++++- .../incflo_apply_nodal_projection.cpp | 50 ++- src/setup/incflo_arrays.cpp | 2 +- 11 files changed, 600 insertions(+), 132 deletions(-) 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..42ce6881 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..18c8818b 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); }