diff --git a/Src/EB/AMReX_EB2.cpp b/Src/EB/AMReX_EB2.cpp index 8786da3c124..ba4c1c3ec07 100644 --- a/Src/EB/AMReX_EB2.cpp +++ b/Src/EB/AMReX_EB2.cpp @@ -267,10 +267,10 @@ int comp_max_crse_level (Box cdomain, const Box& domain) { int ilev; for (ilev = 0; ilev < 30; ++ilev) { - if (cdomain.contains(domain)) break; + if (cdomain.contains(domain)) { break; } cdomain.refine(2); } - if (cdomain != domain) ilev = -1; + if (cdomain != domain) { ilev = -1; } return ilev; } } diff --git a/Src/EB/AMReX_EB2_2D_C.H b/Src/EB/AMReX_EB2_2D_C.H index 873467a9143..529d513163b 100644 --- a/Src/EB/AMReX_EB2_2D_C.H +++ b/Src/EB/AMReX_EB2_2D_C.H @@ -271,10 +271,10 @@ void build_cellflag_from_ap (int i, int j, Array4 const& cflag, // By default, all neighbors are already set. auto flg = cflag(i,j,k); - if (apx(i ,j ,k) == 0.0_rt) flg.setDisconnected(-1, 0, 0); - if (apx(i+1,j ,k) == 0.0_rt) flg.setDisconnected( 1, 0, 0); - if (apy(i ,j ,k) == 0.0_rt) flg.setDisconnected( 0,-1, 0); - if (apy(i ,j+1,k) == 0.0_rt) flg.setDisconnected( 0, 1, 0); + if (apx(i ,j ,k) == 0.0_rt) { flg.setDisconnected(-1, 0, 0); } + if (apx(i+1,j ,k) == 0.0_rt) { flg.setDisconnected( 1, 0, 0); } + if (apy(i ,j ,k) == 0.0_rt) { flg.setDisconnected( 0,-1, 0); } + if (apy(i ,j+1,k) == 0.0_rt) { flg.setDisconnected( 0, 1, 0); } if ((apx(i,j ,k) == 0.0_rt || apy(i-1,j,k) == 0.0_rt) && (apx(i,j-1,k) == 0.0_rt || apy(i ,j,k) == 0.0_rt)) diff --git a/Src/EB/AMReX_EB2_2D_C.cpp b/Src/EB/AMReX_EB2_2D_C.cpp index 5144ac09a42..b77b2ebd00e 100644 --- a/Src/EB/AMReX_EB2_2D_C.cpp +++ b/Src/EB/AMReX_EB2_2D_C.cpp @@ -278,10 +278,10 @@ int build_faces (Box const& bx, Array4 const& cell, else { int ncuts = 0; - if (fx(i ,j ,0) == Type::irregular) ++ncuts; - if (fx(i+1,j ,0) == Type::irregular) ++ncuts; - if (fy(i ,j ,0) == Type::irregular) ++ncuts; - if (fy(i ,j+1,0) == Type::irregular) ++ncuts; + if (fx(i ,j ,0) == Type::irregular) { ++ncuts; } + if (fx(i+1,j ,0) == Type::irregular) { ++ncuts; } + if (fy(i ,j ,0) == Type::irregular) { ++ncuts; } + if (fy(i ,j+1,0) == Type::irregular) { ++ncuts; } if (ncuts > 2) { Gpu::Atomic::Add(dp,1); } @@ -437,10 +437,10 @@ void set_connection_flags (Box const& bxg1, auto flg = cell(i,j,0); - if (fx(i ,j ,0) == Type::covered) flg.setDisconnected(IntVect(-1, 0)); - if (fx(i+1,j ,0) == Type::covered) flg.setDisconnected(IntVect( 1, 0)); - if (fy(i ,j ,0) == Type::covered) flg.setDisconnected(IntVect( 0,-1)); - if (fy(i ,j+1,0) == Type::covered) flg.setDisconnected(IntVect( 0, 1)); + if (fx(i ,j ,0) == Type::covered) { flg.setDisconnected(IntVect(-1, 0)); } + if (fx(i+1,j ,0) == Type::covered) { flg.setDisconnected(IntVect( 1, 0)); } + if (fy(i ,j ,0) == Type::covered) { flg.setDisconnected(IntVect( 0,-1)); } + if (fy(i ,j+1,0) == Type::covered) { flg.setDisconnected(IntVect( 0, 1)); } if (((fx(i,j,0) == Type::covered) || fy(i-1,j,0) == Type::covered) && ((fx(i,j-1,0) == Type::covered) || fy(i,j,0) == Type::covered)) diff --git a/Src/EB/AMReX_EB2_3D_C.H b/Src/EB/AMReX_EB2_3D_C.H index 394b114cae4..10dbdb342d2 100644 --- a/Src/EB/AMReX_EB2_3D_C.H +++ b/Src/EB/AMReX_EB2_3D_C.H @@ -636,107 +636,107 @@ void build_cellflag_from_ap (int i, int j, int k, Array4 const& cfla { flg.setConnected(0,0,0); - if (apx(i ,j,k) != 0.0_rt) flg.setConnected(-1, 0, 0); - if (apx(i+1,j,k) != 0.0_rt) flg.setConnected( 1, 0, 0); - if (apy(i,j ,k) != 0.0_rt) flg.setConnected( 0, -1, 0); - if (apy(i,j+1,k) != 0.0_rt) flg.setConnected( 0, 1, 0); - if (apz(i,j,k ) != 0.0_rt) flg.setConnected( 0, 0, -1); - if (apz(i,j,k+1) != 0.0_rt) flg.setConnected( 0, 0, 1); + if (apx(i ,j,k) != 0.0_rt) { flg.setConnected(-1, 0, 0); } + if (apx(i+1,j,k) != 0.0_rt) { flg.setConnected( 1, 0, 0); } + if (apy(i,j ,k) != 0.0_rt) { flg.setConnected( 0, -1, 0); } + if (apy(i,j+1,k) != 0.0_rt) { flg.setConnected( 0, 1, 0); } + if (apz(i,j,k ) != 0.0_rt) { flg.setConnected( 0, 0, -1); } + if (apz(i,j,k+1) != 0.0_rt) { flg.setConnected( 0, 0, 1); } if ( (apx(i,j,k) != 0.0_rt && apy(i-1,j,k) != 0.0_rt) || (apy(i,j,k) != 0.0_rt && apx(i,j-1,k) != 0.0_rt) ) { flg.setConnected(-1, -1, 0); - if (apz(i-1,j-1,k ) != 0.0_rt) flg.setConnected(-1,-1,-1); - if (apz(i-1,j-1,k+1) != 0.0_rt) flg.setConnected(-1,-1, 1); + if (apz(i-1,j-1,k ) != 0.0_rt) { flg.setConnected(-1,-1,-1); } + if (apz(i-1,j-1,k+1) != 0.0_rt) { flg.setConnected(-1,-1, 1); } } if ( (apx(i+1,j,k) != 0.0_rt && apy(i+1,j ,k) != 0.0_rt) || (apy(i ,j,k) != 0.0_rt && apx(i+1,j-1,k) != 0.0_rt) ) { flg.setConnected(1, -1, 0); - if (apz(i+1,j-1,k ) != 0.0_rt) flg.setConnected(1,-1,-1); - if (apz(i+1,j-1,k+1) != 0.0_rt) flg.setConnected(1,-1, 1); + if (apz(i+1,j-1,k ) != 0.0_rt) { flg.setConnected(1,-1,-1); } + if (apz(i+1,j-1,k+1) != 0.0_rt) { flg.setConnected(1,-1, 1); } } if ( (apx(i,j ,k) != 0.0_rt && apy(i-1,j+1,k) != 0.0_rt) || (apy(i,j+1,k) != 0.0_rt && apx(i ,j+1,k) != 0.0_rt) ) { flg.setConnected(-1, 1, 0); - if (apz(i-1,j+1,k ) != 0.0_rt) flg.setConnected(-1, 1,-1); - if (apz(i-1,j+1,k+1) != 0.0_rt) flg.setConnected(-1, 1, 1); + if (apz(i-1,j+1,k ) != 0.0_rt) { flg.setConnected(-1, 1,-1); } + if (apz(i-1,j+1,k+1) != 0.0_rt) { flg.setConnected(-1, 1, 1); } } if ( (apx(i+1,j ,k) != 0.0_rt && apy(i+1,j+1,k) != 0.0_rt) || (apy(i ,j+1,k) != 0.0_rt && apx(i+1,j+1,k) != 0.0_rt) ) { flg.setConnected(1, 1, 0); - if (apz(i+1,j+1,k ) != 0.0_rt) flg.setConnected(1, 1,-1); - if (apz(i+1,j+1,k+1) != 0.0_rt) flg.setConnected(1, 1, 1); + if (apz(i+1,j+1,k ) != 0.0_rt) { flg.setConnected(1, 1,-1); } + if (apz(i+1,j+1,k+1) != 0.0_rt) { flg.setConnected(1, 1, 1); } } if ( (apx(i,j,k) != 0.0_rt && apz(i-1,j,k ) != 0.0_rt) || (apz(i,j,k) != 0.0_rt && apx(i ,j,k-1) != 0.0_rt) ) { flg.setConnected(-1, 0, -1); - if (apy(i-1,j ,k-1) != 0.0_rt) flg.setConnected(-1,-1,-1); - if (apy(i-1,j+1,k-1) != 0.0_rt) flg.setConnected(-1, 1,-1); + if (apy(i-1,j ,k-1) != 0.0_rt) { flg.setConnected(-1,-1,-1); } + if (apy(i-1,j+1,k-1) != 0.0_rt) { flg.setConnected(-1, 1,-1); } } if ( (apx(i+1,j,k) != 0.0_rt && apz(i+1,j,k ) != 0.0_rt) || (apz(i ,j,k) != 0.0_rt && apx(i+1,j,k-1) != 0.0_rt) ) { flg.setConnected(1, 0, -1); - if (apy(i+1,j ,k-1) != 0.0_rt) flg.setConnected(1,-1,-1); - if (apy(i+1,j+1,k-1) != 0.0_rt) flg.setConnected(1, 1,-1); + if (apy(i+1,j ,k-1) != 0.0_rt) { flg.setConnected(1,-1,-1); } + if (apy(i+1,j+1,k-1) != 0.0_rt) { flg.setConnected(1, 1,-1); } } if ( (apx(i,j,k ) != 0.0_rt && apz(i-1,j,k+1) != 0.0_rt) || (apz(i,j,k+1) != 0.0_rt && apx(i ,j,k+1) != 0.0_rt) ) { flg.setConnected(-1, 0, 1); - if (apy(i-1,j ,k+1) != 0.0_rt) flg.setConnected(-1,-1, 1); - if (apy(i-1,j+1,k+1) != 0.0_rt) flg.setConnected(-1, 1, 1); + if (apy(i-1,j ,k+1) != 0.0_rt) { flg.setConnected(-1,-1, 1); } + if (apy(i-1,j+1,k+1) != 0.0_rt) { flg.setConnected(-1, 1, 1); } } if ( (apx(i+1,j,k ) != 0.0_rt && apz(i+1,j,k+1) != 0.0_rt) || (apz(i ,j,k+1) != 0.0_rt && apx(i+1,j,k+1) != 0.0_rt) ) { flg.setConnected(1, 0, 1); - if (apy(i+1,j ,k+1) != 0.0_rt) flg.setConnected(1,-1, 1); - if (apy(i+1,j+1,k+1) != 0.0_rt) flg.setConnected(1, 1, 1); + if (apy(i+1,j ,k+1) != 0.0_rt) { flg.setConnected(1,-1, 1); } + if (apy(i+1,j+1,k+1) != 0.0_rt) { flg.setConnected(1, 1, 1); } } if ( (apy(i,j,k) != 0.0_rt && apz(i,j-1,k ) != 0.0_rt) || (apz(i,j,k) != 0.0_rt && apy(i,j ,k-1) != 0.0_rt) ) { flg.setConnected(0, -1, -1); - if (apx(i ,j-1,k-1) != 0.0_rt) flg.setConnected(-1,-1,-1); - if (apx(i+1,j-1,k-1) != 0.0_rt) flg.setConnected( 1,-1,-1); + if (apx(i ,j-1,k-1) != 0.0_rt) { flg.setConnected(-1,-1,-1); } + if (apx(i+1,j-1,k-1) != 0.0_rt) { flg.setConnected( 1,-1,-1); } } if ( (apy(i,j+1,k) != 0.0_rt && apz(i,j+1,k ) != 0.0_rt) || (apz(i,j ,k) != 0.0_rt && apy(i,j+1,k-1) != 0.0_rt) ) { flg.setConnected(0, 1, -1); - if (apx(i ,j+1,k-1) != 0.0_rt) flg.setConnected(-1, 1,-1); - if (apx(i+1,j+1,k-1) != 0.0_rt) flg.setConnected( 1, 1,-1); + if (apx(i ,j+1,k-1) != 0.0_rt) { flg.setConnected(-1, 1,-1); } + if (apx(i+1,j+1,k-1) != 0.0_rt) { flg.setConnected( 1, 1,-1); } } if ( (apy(i,j,k ) != 0.0_rt && apz(i,j-1,k+1) != 0.0_rt) || (apz(i,j,k+1) != 0.0_rt && apy(i,j ,k+1) != 0.0_rt) ) { flg.setConnected(0, -1, 1); - if (apx(i ,j-1,k+1) != 0.0_rt) flg.setConnected(-1,-1, 1); - if (apx(i+1,j-1,k+1) != 0.0_rt) flg.setConnected( 1,-1, 1); + if (apx(i ,j-1,k+1) != 0.0_rt) { flg.setConnected(-1,-1, 1); } + if (apx(i+1,j-1,k+1) != 0.0_rt) { flg.setConnected( 1,-1, 1); } } if ( (apy(i,j+1,k ) != 0.0_rt && apz(i,j+1,k+1) != 0.0_rt) || (apz(i,j ,k+1) != 0.0_rt && apy(i,j+1,k+1) != 0.0_rt) ) { flg.setConnected(0, 1, 1); - if (apx(i ,j+1,k+1) != 0.0_rt) flg.setConnected(-1, 1, 1); - if (apx(i+1,j+1,k+1) != 0.0_rt) flg.setConnected( 1, 1, 1); + if (apx(i ,j+1,k+1) != 0.0_rt) { flg.setConnected(-1, 1, 1); } + if (apx(i+1,j+1,k+1) != 0.0_rt) { flg.setConnected( 1, 1, 1); } } } diff --git a/Src/EB/AMReX_EB2_GeometryShop.H b/Src/EB/AMReX_EB2_GeometryShop.H index 6b5e255286e..974738e94f1 100644 --- a/Src/EB/AMReX_EB2_GeometryShop.H +++ b/Src/EB/AMReX_EB2_GeometryShop.H @@ -123,7 +123,7 @@ BrentRootFinder (GpuArray const& lo, } // Check whether in bounds - if (p > 0) q = -q; + if (p > 0) { q = -q; } p = std::abs(p); @@ -158,8 +158,8 @@ BrentRootFinder (GpuArray const& lo, } else { - if (xm < 0) bPt[rangedir] = bPt[rangedir] - tol1; - else bPt[rangedir] = bPt[rangedir] + tol1; + if (xm < 0) { bPt[rangedir] = bPt[rangedir] - tol1; } + else { bPt[rangedir] = bPt[rangedir] + tol1; } } fb = IF_f(f, bPt); @@ -220,7 +220,7 @@ public: } else { ++nfluid; } - if (nbody > 0 && nfluid > 0) return mixedcells; + if (nbody > 0 && nfluid > 0) { return mixedcells; } } } } diff --git a/Src/EB/AMReX_EB2_Level.cpp b/Src/EB/AMReX_EB2_Level.cpp index 9289f5d77f9..795338db135 100644 --- a/Src/EB/AMReX_EB2_Level.cpp +++ b/Src/EB/AMReX_EB2_Level.cpp @@ -162,7 +162,7 @@ Level::coarsenFromFine (Level& fineLevel, bool fill_boundary) ParallelDescriptor::ReduceBoolOr(b); mvmc_error = b; } - if (mvmc_error) return mvmc_error; + if (mvmc_error) { return mvmc_error; } const int ng = 2; m_cellflag.define(m_grids, m_dmap, 1, ng); @@ -241,9 +241,9 @@ Level::coarsenFromFine (Level& fineLevel, bool fill_boundary) vfrac(i,j,k) = 0.0; cflag(i,j,k) = EBCellFlag::TheCoveredCell(); } - AMREX_D_TERM(if (xbx.contains(cell)) apx(i,j,k) = 0.0;, - if (ybx.contains(cell)) apy(i,j,k) = 0.0;, - if (zbx.contains(cell)) apz(i,j,k) = 0.0;); + AMREX_D_TERM(if (xbx.contains(cell)) { apx(i,j,k) = 0.0; }, + if (ybx.contains(cell)) { apy(i,j,k) = 0.0; }, + if (zbx.contains(cell)) { apz(i,j,k) = 0.0; }) }); } } @@ -483,7 +483,7 @@ void Level::fillVolFrac (MultiFab& vfrac, const Geometry& geom) const { vfrac.setVal(1.0); - if (isAllRegular()) return; + if (isAllRegular()) { return; } vfrac.ParallelCopy(m_volfrac,0,0,1,0,vfrac.nGrow(),geom.periodicity()); @@ -715,7 +715,7 @@ Level::fillAreaFrac (Array const& a_areafrac, const Ge a_areafrac[idim]->setVal(1.0); } - if (isAllRegular()) return; + if (isAllRegular()) { return; } for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { diff --git a/Src/EB/AMReX_EBFluxRegister.cpp b/Src/EB/AMReX_EBFluxRegister.cpp index ed2c42f758c..7ba67a562c7 100644 --- a/Src/EB/AMReX_EBFluxRegister.cpp +++ b/Src/EB/AMReX_EBFluxRegister.cpp @@ -166,7 +166,7 @@ EBFluxRegister::FineAdd (const MFIter& mfi, const int li = mfi.LocalIndex(); Vector& cfp_fabs = m_cfp_fab[li]; - if (cfp_fabs.empty()) return; + if (cfp_fabs.empty()) { return; } const Box& tbx = mfi.tilebox(); BL_ASSERT(tbx.cellCentered()); @@ -340,7 +340,7 @@ EBFluxRegister::Reflux (MultiFab& crse_state, const amrex::MultiFab& crse_vfrac, const Box& gdomain = m_crse_geom.growPeriodicDomain(1); MFItInfo info; - if (Gpu::notInLaunchRegion()) info.EnableTiling().SetDynamic(true); + if (Gpu::notInLaunchRegion()) { info.EnableTiling().SetDynamic(true); } #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif diff --git a/Src/EB/AMReX_EBMultiFabUtil.cpp b/Src/EB/AMReX_EBMultiFabUtil.cpp index a5202c7d805..532653b6025 100644 --- a/Src/EB/AMReX_EBMultiFabUtil.cpp +++ b/Src/EB/AMReX_EBMultiFabUtil.cpp @@ -27,7 +27,7 @@ void EB_set_covered (MultiFab& mf, int icomp, int ncomp, int ngrow, Real val) { const auto *const factory = dynamic_cast(&(mf.Factory())); - if (factory == nullptr) return; + if (factory == nullptr) { return; } const auto& flags = factory->getMultiEBCellFlagFab(); AMREX_ALWAYS_ASSERT(mf.ixType().cellCentered() || mf.ixType().nodeCentered()); @@ -69,7 +69,7 @@ void EB_set_covered (MultiFab& mf, int icomp, int ncomp, int ngrow, const Vector& a_vals) { const auto *const factory = dynamic_cast(&(mf.Factory())); - if (factory == nullptr) return; + if (factory == nullptr) { return; } const auto& flags = factory->getMultiEBCellFlagFab(); AMREX_ALWAYS_ASSERT(mf.ixType().cellCentered() || mf.ixType().nodeCentered()); @@ -109,7 +109,7 @@ void EB_set_covered_faces (const Array& umac, Real val) { const auto *const factory = dynamic_cast(&(umac[0]->Factory())); - if (factory == nullptr) return; + if (factory == nullptr) { return; } const auto& area = factory->getAreaFrac(); const auto& flags = factory->getMultiEBCellFlagFab(); @@ -220,7 +220,7 @@ void EB_set_covered_faces (const Array& umac, const int scomp, const int ncomp, const Vector& a_vals ) { const auto *const factory = dynamic_cast(&(umac[0]->Factory())); - if (factory == nullptr) return; + if (factory == nullptr) { return; } const auto& area = factory->getAreaFrac(); const auto& flags = factory->getMultiEBCellFlagFab(); @@ -299,7 +299,7 @@ EB_set_covered_faces (const Array& umac, const int sco for (int k = lo.z; k <= hi.z; ++k) { for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { - if (ax(i,j,k) == 0.0) u(i,j,k,n) = vals[n]; + if (ax(i,j,k) == Real(0.0)) { u(i,j,k,n) = vals[n]; } }}} } } @@ -311,7 +311,7 @@ EB_set_covered_faces (const Array& umac, const int sco for (int k = lo.z; k <= hi.z; ++k) { for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { - if (ay(i,j,k) == 0.0) v(i,j,k,n) = vals[n]; + if (ay(i,j,k) == Real(0.0)) { v(i,j,k,n) = vals[n]; } }}} } } @@ -323,7 +323,7 @@ EB_set_covered_faces (const Array& umac, const int sco for (int k = lo.z; k <= hi.z; ++k) { for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { - if (az(i,j,k) == 0.0) w(i,j,k,n) = vals[n]; + if (az(i,j,k) == Real(0.0)) { w(i,j,k,n) = vals[n]; } }}} } } @@ -642,7 +642,7 @@ void EB_average_down_boundaries (const MultiFab& fine, MultiFab& crse, if (isMFIterSafe(fine, crse)) { MFItInfo info; - if (Gpu::notInLaunchRegion()) info.EnableTiling().SetDynamic(true); + if (Gpu::notInLaunchRegion()) { info.EnableTiling().SetDynamic(true); } #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -712,7 +712,7 @@ void EB_computeDivergence (MultiFab& divu, const Array dxinv = geom.InvCellSizeArray(); MFItInfo info; - if (Gpu::notInLaunchRegion()) info.EnableTiling().SetDynamic(true); + if (Gpu::notInLaunchRegion()) { info.EnableTiling().SetDynamic(true); } #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -774,7 +774,7 @@ void EB_computeDivergence (MultiFab& divu, const Array Real(1.e-8)) + if (std::abs(fcx(i,j,k,0)) > Real(1.e-8)) { jj = (fcx(i,j,k,0) < Real(0.0)) ? j - 1 : j + 1; - else if (apx(i,j-1,k) > Real(0.)) + } else if (apx(i,j-1,k) > Real(0.)) { jj = j-1; - else + } else { jj = j+1; + } - if (std::abs(fcx(i,j,k,1)) > Real(1.e-8)) + if (std::abs(fcx(i,j,k,1)) > Real(1.e-8)) { kk = (fcx(i,j,k,1) < Real(0.0)) ? k - 1 : k + 1; - else if (apx(i,j,k-1) > Real(0.)) + } else if (apx(i,j,k-1) > Real(0.)) { kk = k-1; - else + } else { kk = k+1; + } // If any of these cells has zero volume we don't want to use this stencil Real test_zero = cvol(i-1,jj,k ) * cvol(i-1,j ,kk) * cvol(i-1,jj,kk) * @@ -702,12 +704,11 @@ void eb_interp_centroid2facecent_x (Box const& ubx, Real test_zero_jkalt = cvol(i-1,jalt,k) * cvol(i-1,j,kalt) * cvol(i-1,jalt,kalt) * cvol(i ,jalt,k) * cvol(i ,j,kalt) * cvol(i ,jalt,kalt); - if (test_zero_jalt > Real(0.)) + if (test_zero_jalt > Real(0.)) { jj = jalt; - else if (test_zero_kalt > Real(0.)) + } else if (test_zero_kalt > Real(0.)) { kk = kalt; - else if (test_zero_jkalt > Real(0.)) - { + } else if (test_zero_jkalt > Real(0.)) { jj = jalt; kk = kalt; } @@ -766,15 +767,17 @@ void eb_interp_centroid2facecent_x (Box const& ubx, // This is the location of the face centroid relative to the central node // Recall fcx holds (y,z) of the x-face centroid as components ( /0/1) - if (j < jj) + if (j < jj) { y = Real(-0.5) + fcx(i,j,k,0); // (j,k) is in lower half of stencil so y < 0 - else + } else { y = Real(0.5) + fcx(i,j,k,0); // (j,k) is in upper half of stencil so y > 0 + } - if (k < kk) + if (k < kk) { z = Real(-0.5) + fcx(i,j,k,1); // (j,k) is in lower half of stencil so z < 0 - else + } else { z = Real(0.5) + fcx(i,j,k,1); // (j,k) is in upper half of stencil so z > 0 + } if (j < jj && k < kk) // (j,k) is lower left, (j+1,k+1) is upper right { @@ -889,19 +892,21 @@ void eb_interp_centroid2facecent_y (Box const& vbx, // We must add additional tests to avoid the case where fcy is very close to zero, but i-1/i+1 or k-1/k+1 // might be covered cells -- this can happen when the EB is exactly aligned with the grid planes int ii,kk; - if (std::abs(fcy(i,j,k,0)) > Real(1.e-8)) + if (std::abs(fcy(i,j,k,0)) > Real(1.e-8)) { ii = (fcy(i,j,k,0) < Real(0.0)) ? i - 1 : i + 1; - else if (apy(i-1,j,k) > Real(0.)) + } else if (apy(i-1,j,k) > Real(0.)) { ii = i-1; - else + } else { ii = i+1; + } - if (std::abs(fcy(i,j,k,1)) > Real(1.e-8)) + if (std::abs(fcy(i,j,k,1)) > Real(1.e-8)) { kk = (fcy(i,j,k,1) < Real(0.0)) ? k - 1 : k + 1; - else if (apy(i,j,k-1) > Real(0.)) + } else if (apy(i,j,k-1) > Real(0.)) { kk = k-1; - else + } else { kk = k+1; + } // If any of these cells has zero volume we don't want to use this stencil Real test_zero = cvol(ii,j-1,k) * cvol(i,j-1,kk) * cvol(ii,j-1,kk) * @@ -920,12 +925,11 @@ void eb_interp_centroid2facecent_y (Box const& vbx, Real test_zero_ikalt = cvol(ialt,j-1,k) * cvol(i,j-1,kalt) * cvol(ialt,j-1,kalt) * cvol(ialt,j ,k) * cvol(i,j ,kalt) * cvol(ialt,j ,kalt); - if (test_zero_ialt > Real(0.)) + if (test_zero_ialt > Real(0.)) { ii = ialt; - else if (test_zero_kalt > Real(0.)) + } else if (test_zero_kalt > Real(0.)) { kk = kalt; - else if (test_zero_ikalt > Real(0.)) - { + } else if (test_zero_ikalt > Real(0.)) { ii = ialt; kk = kalt; } @@ -981,15 +985,17 @@ void eb_interp_centroid2facecent_y (Box const& vbx, // This is the location of the face centroid relative to the central node // Recall fcy holds (x,z) of the x-face centroid as components (0/ /1) - if (i < ii) + if (i < ii) { x = Real(-0.5) + fcy(i,j,k,0); // (i,k) is in lower half of stencil so x < 0 - else + } else { x = Real(0.5) + fcy(i,j,k,0); // (i,k) is in upper half of stencil so x > 0 + } - if (k < kk) + if (k < kk) { z = Real(-0.5) + fcy(i,j,k,1); // (i,k) is in lower half of stencil so z < 0 - else + } else { z = Real(0.5) + fcy(i,j,k,1); // (i,k) is in upper half of stencil so z > 0 + } if (i < ii && k < kk) // (i,k) is lower left, (i+1,k+1) is upper right { @@ -1103,19 +1109,21 @@ void eb_interp_centroid2facecent_z (Box const& wbx, // We must add additional tests to avoid the case where fcz is very close to zero, but i-1/i+1 or j-1/j+1 // might be covered cells -- this can happen when the EB is exactly aligned with the grid planes int ii,jj; - if (std::abs(fcz(i,j,k,0)) > Real(1.e-8)) + if (std::abs(fcz(i,j,k,0)) > Real(1.e-8)) { ii = (fcz(i,j,k,0) < Real(0.0)) ? i - 1 : i + 1; - else if (apz(i-1,j,k) > Real(0.)) + } else if (apz(i-1,j,k) > Real(0.)) { ii = i-1; - else + } else { ii = i+1; + } - if (std::abs(fcz(i,j,k,1)) > Real(1.e-8)) + if (std::abs(fcz(i,j,k,1)) > Real(1.e-8)) { jj = (fcz(i,j,k,1) < Real(0.0)) ? j - 1 : j + 1; - else if (apz(i,j-1,k) > Real(0.)) + } else if (apz(i,j-1,k) > Real(0.)) { jj = j-1; - else + } else { jj = j+1; + } // If any of these cells has zero volume we don't want to use this stencil Real test_zero = cvol(ii,j,k-1) * cvol(i,jj,k-1) * cvol(ii,jj,k-1) * @@ -1134,12 +1142,11 @@ void eb_interp_centroid2facecent_z (Box const& wbx, Real test_zero_ijalt = cvol(ialt,j,k-1) * cvol(i,jalt,k-1) * cvol(ialt,jalt,k-1) * cvol(ialt,j,k ) * cvol(i,jalt,k ) * cvol(ialt,jalt,k ); - if (test_zero_ialt > Real(0.)) + if (test_zero_ialt > Real(0.)) { ii = ialt; - else if (test_zero_jalt > Real(0.)) + } else if (test_zero_jalt > Real(0.)) { jj = jalt; - else if (test_zero_ijalt > Real(0.)) - { + } else if (test_zero_ijalt > Real(0.)) { ii = ialt; jj = jalt; } @@ -1195,15 +1202,17 @@ void eb_interp_centroid2facecent_z (Box const& wbx, // This is the location of the face centroid relative to the central node // Recall fcz holds (x,y) of the x-face centroid as components (0/1/ ) - if (i < ii) + if (i < ii) { x = Real(-0.5) + fcz(i,j,k,0); // (i,k) is in lower half of stencil so x < 0 - else + } else { x = Real(0.5) + fcz(i,j,k,0); // (i,k) is in upper half of stencil so x > 0 + } - if (j < jj) + if (j < jj) { y = Real(-0.5) + fcz(i,j,k,1); // (i,k) is in lower half of stencil so z < 0 - else + } else { y = Real(0.5) + fcz(i,j,k,1); // (i,k) is in upper half of stencil so z > 0 + } if (i < ii && j < jj) // (i,j) is lower left, (i+1,j+1) is upper right { diff --git a/Src/EB/AMReX_EBToPVD.cpp b/Src/EB/AMReX_EBToPVD.cpp index dc6e93d5fb2..3835af4d3df 100644 --- a/Src/EB/AMReX_EBToPVD.cpp +++ b/Src/EB/AMReX_EBToPVD.cpp @@ -246,12 +246,14 @@ void EBToPVD::reorder_polygon(const std::vector>& lpoints, int longest = 2; if(std::abs(lnormal[0]) > std::abs(lnormal[1])) { - if(std::abs(lnormal[0]) > std::abs(lnormal[2])) + if(std::abs(lnormal[0]) > std::abs(lnormal[2])) { longest = 0; + } } else { - if(std::abs(lnormal[1]) > std::abs(lnormal[2])) + if(std::abs(lnormal[1]) > std::abs(lnormal[2])) { longest = 1; + } } for(int i = 1; i <= lconnect[0]; ++i) { @@ -436,14 +438,15 @@ void EBToPVD::EBGridCoverage(const int myID, const Real* problo, const Real* dx, for(int j = lo.y; j <= hi.y; ++j) { for(int i = lo.x; i <= hi.x; ++i) { - if(flag(i,j,k).isSingleValued()) + if(flag(i,j,k).isSingleValued()) { lc1 = lc1 + 1; + } } } }; ++m_grid; - if(lc1 == 0) return; + if(lc1 == 0) { return; } std::stringstream ss; ss << std::setw(4) << std::setfill('0') << myID; diff --git a/Src/EB/AMReX_EB_FluxRedistribute.cpp b/Src/EB/AMReX_EB_FluxRedistribute.cpp index f0f38aee929..d5d2b22c0c7 100644 --- a/Src/EB/AMReX_EB_FluxRedistribute.cpp +++ b/Src/EB/AMReX_EB_FluxRedistribute.cpp @@ -102,17 +102,19 @@ amrex_flux_redistribute ( #if (AMREX_SPACEDIM == 2) int kk(0); #else - for (int kk = -1; kk <= 1; kk++) + for (int kk = -1; kk <= 1; kk++) { #endif - for (int jj = -1; jj <= 1; jj++) - for (int ii = -1; ii <= 1; ii++) - if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) && - dbox.contains(IntVect(AMREX_D_DECL(i+ii,j+jj,k+kk)))) - { - Real wted_frac = vfrac(i+ii,j+jj,k+kk) * wt(i+ii,j+jj,k+kk) * mask(i+ii,j+jj,k+kk); - vtot += wted_frac; - divnc += wted_frac * divc(i+ii,j+jj,k+kk,n); - } + for (int jj = -1; jj <= 1; jj++) { + for (int ii = -1; ii <= 1; ii++) { + if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) && + dbox.contains(IntVect(AMREX_D_DECL(i+ii,j+jj,k+kk)))) + { + Real wted_frac = vfrac(i+ii,j+jj,k+kk) * wt(i+ii,j+jj,k+kk) * mask(i+ii,j+jj,k+kk); + vtot += wted_frac; + divnc += wted_frac * divc(i+ii,j+jj,k+kk,n); + } + AMREX_D_TERM(},},}) + divnc /= vtot; // We need to multiply by mask to make sure optmp is zero for cells @@ -135,17 +137,19 @@ amrex_flux_redistribute ( #if (AMREX_SPACEDIM == 2) int kk(0); #else - for (int kk = -1; kk <= 1; kk++) + for (int kk = -1; kk <= 1; kk++) { #endif - for (int jj = -1; jj <= 1; jj++) - for (int ii = -1; ii <= 1; ii++) - if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) && - dbox.contains(IntVect(AMREX_D_DECL(i+ii,j+jj,k+kk)))) - { - Real unwted_frac = vfrac(i+ii,j+jj,k+kk) * mask(i+ii,j+jj,k+kk); - vtot += unwted_frac; - divnc += unwted_frac*divc(i+ii,j+jj,k+kk,n); - } + for (int jj = -1; jj <= 1; jj++) { + for (int ii = -1; ii <= 1; ii++) { + if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) && + dbox.contains(IntVect(AMREX_D_DECL(i+ii,j+jj,k+kk)))) + { + Real unwted_frac = vfrac(i+ii,j+jj,k+kk) * mask(i+ii,j+jj,k+kk); + vtot += unwted_frac; + divnc += unwted_frac*divc(i+ii,j+jj,k+kk,n); + } + AMREX_D_TERM(},},}) + divnc /= vtot; // We need to multiply by mask to make sure optmp is zero for cells @@ -172,14 +176,16 @@ amrex_flux_redistribute ( #if (AMREX_SPACEDIM == 2) int kk(0); #else - for (int kk = -1; kk <= 1; kk++) + for (int kk = -1; kk <= 1; kk++) { #endif - for (int jj = -1; jj <= 1; jj++) - for (int ii = -1; ii <= 1; ii++) - if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) ) - { - wtot += vfrac(i+ii,j+jj,k+kk)*wt(i+ii,j+jj,k+kk)* mask(i+ii,j+jj,k+kk); - } + for (int jj = -1; jj <= 1; jj++) { + for (int ii = -1; ii <= 1; ii++) { + if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) ) + { + wtot += vfrac(i+ii,j+jj,k+kk)*wt(i+ii,j+jj,k+kk)* mask(i+ii,j+jj,k+kk); + } + AMREX_D_TERM(},},}) + #ifdef AMREX_USE_FLOAT wtot = Real(1.0)/(wtot + Real(1.e-30)); #else @@ -212,71 +218,72 @@ amrex_flux_redistribute ( #else ( (i >= bx_ilo) && (i <= bx_ihi) && (j >= bx_jlo) && (j <= bx_jhi) && (k >= bx_klo) && (k <= bx_khi) ); #endif - if (inside) as_fine_valid_cell = true; + if (inside) { as_fine_valid_cell = true; } as_fine_ghost_cell = (levmsk(i,j,k) == level_mask_not_covered); // not covered by other grids } #if (AMREX_SPACEDIM == 2) kk = 0; #else - for (int kk = -1; kk <= 1; kk++) + for (int kk = -1; kk <= 1; kk++) { #endif - for (int jj = -1; jj <= 1; jj++) - for (int ii = -1; ii <= 1; ii++) - if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) ) - { - int iii = i + ii; - int jjj = j + jj; - int kkk = k + kk; - - Real drho = delm(i,j,k,n)*wtot*wt(iii,jjj,kkk)* mask(iii,jjj,kkk) ; - Gpu::Atomic::Add(&optmp(iii,jjj,kkk,n), drho); - - valid_dst_cell = ( (iii >= bx_ilo) && (iii <= bx_ihi) && - (jjj >= bx_jlo) && (jjj <= bx_jhi) ); + for (int jj = -1; jj <= 1; jj++) { + for (int ii = -1; ii <= 1; ii++) { + if ( (ii != 0 || jj != 0 || kk != 0) && flag(i,j,k).isConnected(ii,jj,kk) ) + { + int iii = i + ii; + int jjj = j + jj; + int kkk = k + kk; + + Real drho = delm(i,j,k,n)*wtot*wt(iii,jjj,kkk)* mask(iii,jjj,kkk) ; + Gpu::Atomic::Add(&optmp(iii,jjj,kkk,n), drho); + + valid_dst_cell = ( (iii >= bx_ilo) && (iii <= bx_ihi) && + (jjj >= bx_jlo) && (jjj <= bx_jhi) ); #if (AMREX_SPACEDIM == 3) - valid_dst_cell &= ( (kkk >= bx_klo) && (kkk <= bx_khi) ); + valid_dst_cell &= ( (kkk >= bx_klo) && (kkk <= bx_khi) ); #endif - if (as_crse_crse_cell) - { - if ( (rr_flag_crse(iii,jjj,kkk) == amrex_yafluxreg_fine_cell) && - (vfrac(i,j,k) > reredistribution_threshold) ) - { - Gpu::Atomic::Add(&rr_drho_crse(i,j,k,n), - dt*drho*(vfrac(iii,jjj,kkk)/vfrac(i,j,k))); - } - } - - if (as_crse_covered_cell && valid_dst_cell) - { - if ( (rr_flag_crse(iii,jjj,kkk) == amrex_yafluxreg_crse_fine_boundary_cell) && - (vfrac(iii,jjj,kkk) > reredistribution_threshold) ) - { + if (as_crse_crse_cell) + { + if ( (rr_flag_crse(iii,jjj,kkk) == amrex_yafluxreg_fine_cell) && + (vfrac(i,j,k) > reredistribution_threshold) ) + { + Gpu::Atomic::Add(&rr_drho_crse(i,j,k,n), + dt*drho*(vfrac(iii,jjj,kkk)/vfrac(i,j,k))); + } + } + + if (as_crse_covered_cell && valid_dst_cell) + { + if ( (rr_flag_crse(iii,jjj,kkk) == amrex_yafluxreg_crse_fine_boundary_cell) && + (vfrac(iii,jjj,kkk) > reredistribution_threshold) ) + { // recipient is a crse/fine boundary cell - Gpu::Atomic::Add(&rr_drho_crse(iii,jjj,kkk,n), -dt*drho); - } - } - - if (as_fine_valid_cell && !valid_dst_cell) - { - Gpu::Atomic::Add(&dm_as_fine(iii,jjj,kkk,n), dt*drho*vfrac(iii,jjj,kkk)); - } - - if (as_fine_ghost_cell && valid_dst_cell) - { - Gpu::Atomic::Add(&dm_as_fine(i,j,k,n), -dt*drho*vfrac(iii,jjj,kkk)); - } - - } // isConnected + Gpu::Atomic::Add(&rr_drho_crse(iii,jjj,kkk,n), -dt*drho); + } + } + + if (as_fine_valid_cell && !valid_dst_cell) + { + Gpu::Atomic::Add(&dm_as_fine(iii,jjj,kkk,n), dt*drho*vfrac(iii,jjj,kkk)); + } + + if (as_fine_ghost_cell && valid_dst_cell) + { + Gpu::Atomic::Add(&dm_as_fine(i,j,k,n), -dt*drho*vfrac(iii,jjj,kkk)); + } + } // isConnected + AMREX_D_TERM(},},}) } // isSingleValued }); amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if (!flag(i,j,k).isCovered()) + if (!flag(i,j,k).isCovered()) { dqdt(i,j,k,icomp+n) = divc(i,j,k,n) + optmp(i,j,k,n); + } }); diff --git a/Src/EB/AMReX_EB_LeastSquares_2D_K.H b/Src/EB/AMReX_EB_LeastSquares_2D_K.H index fb2d6e0673b..7ebfff7e34d 100644 --- a/Src/EB/AMReX_EB_LeastSquares_2D_K.H +++ b/Src/EB/AMReX_EB_LeastSquares_2D_K.H @@ -38,10 +38,11 @@ void decomp_chol_np6(Array2D& aa) p[ii] = std::sqrt(sum1); } } else { - if (ising == 0) + if (ising == 0) { aa(jj,ii) = sum1 / p[ii]; - else + } else { aa(jj,ii) = 0.0; + } } } } @@ -66,9 +67,11 @@ void cholsol_np6(Array2D& Amatrix, Array1D& b) Array2D AtA; - for (int irow = 0; irow < neq; irow++) - for (int icol = 0; icol < neq; icol++) + for (int irow = 0; irow < neq; irow++) { + for (int icol = 0; icol < neq; icol++) { AtA(irow,icol) = 0.0; + } + } for (int irow = 0; irow < 12; irow++) { @@ -95,50 +98,56 @@ void cholsol_np6(Array2D& Amatrix, Array1D& b) AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (y^2)^T (y^2) } - for (int irow = 0; irow < neq-1; irow++) - for (int icol = irow+1; icol < neq; icol++) + for (int irow = 0; irow < neq-1; irow++) { + for (int icol = irow+1; icol < neq; icol++) { AtA(icol,irow) = AtA(irow,icol); + } + } decomp_chol_np6(AtA); - if (AtA(0,0) > 0.) - b(0) = b(0) / AtA(0,0); - else - b(0) = 0.; + if (AtA(0,0) > 0.) { + b(0) = b(0) / AtA(0,0); + } else { + b(0) = 0.; + } for (int ii = 1; ii < neq; ii++) { - if (AtA(ii,ii) > 0.) - { - for (int jj = 0; jj < ii; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); + if (AtA(ii,ii) > 0.) + { + for (int jj = 0; jj < ii; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = 0.0; - } + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = 0.0; + } } - if (AtA(neq-1,neq-1) > 0.) - b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); - else - b(neq-1) = 0.0; + if (AtA(neq-1,neq-1) > 0.) { + b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); + } else { + b(neq-1) = 0.0; + } for (int ii = neq-2; ii >= 0; ii--) { - if (AtA(ii,ii) > 0.) - { - for (int jj = ii+1; jj < neq; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); - - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = 0.0; - } + if (AtA(ii,ii) > 0.) + { + for (int jj = ii+1; jj < neq; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } + + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = 0.0; + } } } @@ -149,9 +158,11 @@ void cholsol_for_eb(Array2D& Amatrix, Array1D& b) Array2D AtA; - for (int irow = 0; irow < neq; irow++) - for (int icol = 0; icol < neq; icol++) + for (int irow = 0; irow < neq; irow++) { + for (int icol = 0; icol < neq; icol++) { AtA(irow,icol) = 0.0; + } + } for (int irow = 0; irow < 18; irow++) { @@ -178,50 +189,56 @@ void cholsol_for_eb(Array2D& Amatrix, Array1D& b) AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (y^2)^T (y^2) } - for (int irow = 0; irow < neq-1; irow++) - for (int icol = irow+1; icol < neq; icol++) - AtA(icol,irow) = AtA(irow,icol); + for (int irow = 0; irow < neq-1; irow++) { + for (int icol = irow+1; icol < neq; icol++) { + AtA(icol,irow) = AtA(irow,icol); + } + } decomp_chol_np6(AtA); - if (AtA(0,0) > 0.) - b(0) = b(0) / AtA(0,0); - else - b(0) = 0.; + if (AtA(0,0) > 0.) { + b(0) = b(0) / AtA(0,0); + } else { + b(0) = 0.; + } for (int ii = 1; ii < neq; ii++) { - if (AtA(ii,ii) > 0.) - { - for (int jj = 0; jj < ii; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); + if (AtA(ii,ii) > 0.) + { + for (int jj = 0; jj < ii; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = 0.0; - } + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = 0.0; + } } - if (AtA(neq-1,neq-1) > 0.) - b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); - else - b(neq-1) = 0.0; + if (AtA(neq-1,neq-1) > 0.) { + b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); + } else { + b(neq-1) = 0.0; + } for (int ii = neq-2; ii >= 0; ii--) { - if (AtA(ii,ii) > 0.) - { - for (int jj = ii+1; jj < neq; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); - - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = 0.0; - } + if (AtA(ii,ii) > 0.) + { + for (int jj = ii+1; jj < neq; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } + + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = 0.0; + } } } @@ -242,12 +259,14 @@ Real grad_x_of_phi_on_centroids(int i,int j,int k,int n, // Order of column -- first six are cell centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1) // Order of column -- second six are EB centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1) - for (int irow = 0; irow < 12; irow++) - for (int icol = 0; icol < 6; icol++) + for (int irow = 0; irow < 12; irow++) { + for (int icol = 0; icol < 6; icol++) { Amatrix(irow,icol) = 0.0; + } + } // Columns: [e x y x*x x*y y*y] - for (int ii = i-1; ii <= i; ii++) // Normal to face + for (int ii = i-1; ii <= i; ii++) { // Normal to face for (int jj = j-1; jj <= j+1; jj++) // Tangential to face { if (!flag(ii,jj,k).isCovered()) @@ -279,13 +298,14 @@ Real grad_x_of_phi_on_centroids(int i,int j,int k,int n, } } } + } // Make the RHS = A^T v for (int irow = 0; irow < 6; irow++) { rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet - for (int ii = i-1; ii <= i; ii++) // Normal to face + for (int ii = i-1; ii <= i; ii++) { // Normal to face for (int jj = j-1; jj <= j+1; jj++) // Tangential to face { if (!flag(ii,jj,k).isCovered()) @@ -294,10 +314,12 @@ Real grad_x_of_phi_on_centroids(int i,int j,int k,int n, rhs(irow) += Amatrix(a_ind ,irow)* phi(ii,jj,k,n); if (flag(ii,jj,k).isSingleValued() && - is_eb_dirichlet && is_eb_inhomog) + is_eb_dirichlet && is_eb_inhomog) { rhs(irow) += Amatrix(a_ind+6,irow)*phieb(ii,jj,k,n); + } } } + } } cholsol_np6(Amatrix, rhs); @@ -322,12 +344,14 @@ Real grad_y_of_phi_on_centroids(int i,int j,int k,int n, // Order of column -- first six are cell centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1) // Order of column -- second six are EB centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1) - for (int irow = 0; irow < 12; irow++) - for (int icol = 0; icol < 6; icol++) + for (int irow = 0; irow < 12; irow++) { + for (int icol = 0; icol < 6; icol++) { Amatrix(irow,icol) = 0.0; + } + } // Columns: [e x y x*x x*y y*y] - for (int jj = j-1; jj <= j; jj++) // Normal to face + for (int jj = j-1; jj <= j; jj++) { // Normal to face for (int ii = i-1; ii <= i+1; ii++) // Tangential to face { if (!flag(ii,jj,k).isCovered()) @@ -356,23 +380,27 @@ Real grad_y_of_phi_on_centroids(int i,int j,int k,int n, } } } + } // Make the RHS = A^T v for (int irow = 0; irow < 6; irow++) { rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet - for (int jj = j-1; jj <= j; jj++) // Normal to face - for (int ii = i-1; ii <= i+1; ii++) // Tangential to face + for (int jj = j-1; jj <= j; jj++) { // Normal to face + for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face if (!flag(ii,jj,k).isCovered()) { int a_ind = (ii-(i-1)) + 3*(jj-(j-1)); rhs(irow) += Amatrix(a_ind ,irow)* phi(ii,jj,k,n); if (flag(ii,jj,k).isSingleValued() && - is_eb_dirichlet && is_eb_inhomog) + is_eb_dirichlet && is_eb_inhomog) { rhs(irow) += Amatrix(a_ind+6,irow)*phieb(ii,jj,k,n); + } } + } + } } cholsol_np6(Amatrix, rhs); @@ -395,12 +423,14 @@ Real grad_eb_of_phi_on_centroids(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, next 9 are EB centroids - for (int irow = 0; irow < 18; irow++) - for (int icol = 0; icol < 6; icol++) + for (int irow = 0; irow < 18; irow++) { + for (int icol = 0; icol < 6; icol++) { Amatrix(irow,icol) = 0.0; + } + } // Column 0-2: [e x y] - for (int ii = i-1; ii <= i+1; ii++) + for (int ii = i-1; ii <= i+1; ii++) { for (int jj = j-1; jj <= j+1; jj++) { if (!flag(ii,jj,k).isCovered()) @@ -425,6 +455,7 @@ Real grad_eb_of_phi_on_centroids(int i,int j,int k,int n, } } } + } // Columns 3 : [x*x x*y y*y] @@ -440,15 +471,18 @@ Real grad_eb_of_phi_on_centroids(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int ii = i-1; ii <= i+1; ii++) - for (int jj = j-1; jj <= j+1; jj++) + for (int ii = i-1; ii <= i+1; ii++) { + for (int jj = j-1; jj <= j+1; jj++) { if (!flag(ii,jj,k).isCovered()) { int a_ind = (jj-(j-1)) + 3*(ii-(i-1)); rhs(irow) += Amatrix(a_ind,irow) * phi(ii,jj,k,n); - if (flag(ii,jj,k).isSingleValued() && is_eb_inhomog) + if (flag(ii,jj,k).isSingleValued() && is_eb_inhomog) { rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n); + } } + } + } } cholsol_for_eb(Amatrix, rhs); @@ -477,15 +511,17 @@ Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Order of column -- first six are cell centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1) // Order of column -- second six are EB centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1) - for (int irow = 0; irow < 18; irow++) - for (int icol = 0; icol < 6; icol++) + for (int irow = 0; irow < 18; irow++) { + for (int icol = 0; icol < 6; icol++) { Amatrix(irow,icol) = 0.0; + } + } const int im = (i > domhi_x) ? 2 : 1; const int ip = 2 - im; // Columns: [e x y x*x x*y y*y] - for (int ii = i-im; ii <= i+ip; ii++) // Normal to face + for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face for (int jj = j-1; jj <= j+1; jj++) // Tangential to face { @@ -498,8 +534,9 @@ Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n, continue; } - if ( !phi.contains(ii,jj,k) ) + if ( !phi.contains(ii,jj,k) ) { continue; + } if (!flag(ii,jj,k).isCovered()) { @@ -513,17 +550,21 @@ Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n, Real y_off = static_cast(jj-j); if(on_x_face){ - if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) + if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) { continue; - if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) + } + if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) { continue; + } } if(on_y_face){ - if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) + if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) { continue; - if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) + } + if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) { continue; + } } Amatrix(a_ind,0) = 1.0; @@ -551,16 +592,18 @@ Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n, } } } + } for (int irow = 0; irow < 6; irow++) { rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet - for (int ii = i-im; ii <= i+ip; ii++) // Normal to face + for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face for (int jj = j-1; jj <= j+1; jj++) // Tangential to face { - if ( !phi.contains(ii,jj,k) ) + if ( !phi.contains(ii,jj,k) ) { continue; + } if (!flag(ii,jj,k).isCovered()) { @@ -572,10 +615,12 @@ Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n, rhs(irow) += Amatrix(a_ind ,irow)*phi_val; if (flag(ii,jj,k).isSingleValued() && - is_eb_dirichlet && is_eb_inhomog && Amatrix(a_ind+9,irow) != 0.0) + is_eb_dirichlet && is_eb_inhomog && Amatrix(a_ind+9,irow) != 0.0) { rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n); + } } } + } } cholsol_for_eb(Amatrix, rhs); @@ -603,15 +648,17 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Order of column -- first six are cell centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1) // Order of column -- second six are EB centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1) - for (int irow = 0; irow < 18; irow++) - for (int icol = 0; icol < 6; icol++) + for (int irow = 0; irow < 18; irow++) { + for (int icol = 0; icol < 6; icol++) { Amatrix(irow,icol) = 0.0; + } + } const int jm = (j > domhi_y) ? 2 : 1; const int jp = 2 - jm; // Columns: [e x y x*x x*y y*y] - for (int jj = j-jm; jj <= j+jp; jj++) // Normal to face + for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face for (int ii = i-1; ii <= i+1; ii++) // Tangential to face { @@ -620,11 +667,12 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y))) { - continue; + continue; } - if ( !phi.contains(ii,jj,k) ) - continue; + if ( !phi.contains(ii,jj,k) ) { + continue; + } if (!flag(ii,jj,k).isCovered()) { @@ -635,17 +683,21 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, Real y_off = static_cast(jj-j) + 0.5; if(on_x_face){ - if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) + if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) { continue; - if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) + } + if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) { continue; + } } if(on_y_face){ - if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) + if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) { continue; - if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) + } + if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) { continue; + } } Amatrix(a_ind,0) = 1.0; @@ -672,6 +724,7 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, } } } + } // Make the RHS = A^T v @@ -679,10 +732,11 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, { rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet - for (int jj = j-jm; jj <= j+jp; jj++) // Normal to face + for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face - if ( !phi.contains(ii,jj,k) ) + if ( !phi.contains(ii,jj,k) ) { continue; + } if (!flag(ii,jj,k).isCovered()) { int a_ind = (ii-(i-1)) + 3*(jj-(j-jm)); @@ -694,10 +748,11 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, if (flag(ii,jj,k).isSingleValued() && is_eb_dirichlet && is_eb_inhomog && Amatrix(a_ind+9,irow) != Real(0.0)){ - rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n); + rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n); } } } + } } @@ -724,12 +779,14 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, next 9 are EB centroids - for (int irow = 0; irow < 18; irow++) - for (int icol = 0; icol < 6; icol++) + for (int irow = 0; irow < 18; irow++) { + for (int icol = 0; icol < 6; icol++) { Amatrix(irow,icol) = 0.0; + } + } // Column 0-2: [e x y] - for (int ii = i-1; ii <= i+1; ii++) + for (int ii = i-1; ii <= i+1; ii++) { for (int jj = j-1; jj <= j+1; jj++) { @@ -738,7 +795,7 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, // by removing the on_?_face restrictions. if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y))){ - continue; + continue; } @@ -750,17 +807,21 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, Real y_off = static_cast(jj-j); if(on_x_face){ - if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) + if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) { continue; - if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) + } + if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) { continue; + } } if(on_y_face){ - if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) + if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) { continue; - if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) + } + if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) { continue; + } } @@ -781,14 +842,15 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, } } } + } // Columns 3 : [x*x x*y y*y] for (int irow = 0; irow < 18; irow++) { - Amatrix(irow,3) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,5) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,3) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,5) = Amatrix(irow,2) * Amatrix(irow,2); } // Make the RHS = A^T v @@ -796,10 +858,11 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int ii = i-1; ii <= i+1; ii++) + for (int ii = i-1; ii <= i+1; ii++) { for (int jj = j-1; jj <= j+1; jj++) { - if ( !phi.contains(ii,jj,k) ) + if ( !phi.contains(ii,jj,k) ) { continue; + } if (!flag(ii,jj,k).isCovered()) { int a_ind = (jj-(j-1)) + 3*(ii-(i-1)); @@ -810,10 +873,12 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, rhs(irow) += Amatrix(a_ind,irow) * phi_val; - if (flag(ii,jj,k).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+9,irow) != 0.0) + if (flag(ii,jj,k).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+9,irow) != 0.0) { rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n); + } } } + } } cholsol_for_eb(Amatrix, rhs); diff --git a/Src/EB/AMReX_EB_LeastSquares_3D_K.H b/Src/EB/AMReX_EB_LeastSquares_3D_K.H index 83c0f509842..bbdd20bdba8 100644 --- a/Src/EB/AMReX_EB_LeastSquares_3D_K.H +++ b/Src/EB/AMReX_EB_LeastSquares_3D_K.H @@ -38,10 +38,11 @@ void decomp_chol_np10(Array2D& aa) p[ii] = std::sqrt(sum1); } } else { - if (ising == 0) + if (ising == 0) { aa(jj,ii) = sum1 / p[ii]; - else + } else { aa(jj,ii) = Real(0.0); + } } } } @@ -66,9 +67,11 @@ void cholsol_np10(Array2D& Amatrix, Array1D& b) Array2D AtA; - for (int irow = 0; irow < neq; irow++) - for (int icol = 0; icol < neq; icol++) + for (int irow = 0; irow < neq; irow++) { + for (int icol = 0; icol < neq; icol++) { AtA(irow,icol) = Real(0.0); + } + } for (int irow = 0; irow < 36; irow++) { @@ -138,50 +141,56 @@ void cholsol_np10(Array2D& Amatrix, Array1D& b) AtA(9,9) += Amatrix(irow,9)*Amatrix(irow,9); // (z^2)^T z^2 } - for (int irow = 0; irow < neq-1; irow++) - for (int icol = irow+1; icol < neq; icol++) - AtA(icol,irow) = AtA(irow,icol); + for (int irow = 0; irow < neq-1; irow++) { + for (int icol = irow+1; icol < neq; icol++) { + AtA(icol,irow) = AtA(irow,icol); + } + } decomp_chol_np10(AtA); - if (AtA(0,0) > 0.) - b(0) = b(0) / AtA(0,0); - else - b(0) = 0.; + if (AtA(0,0) > 0.) { + b(0) = b(0) / AtA(0,0); + } else { + b(0) = 0.; + } for (int ii = 1; ii < neq; ii++) { - if (AtA(ii,ii) > 0.) - { - for (int jj = 0; jj < ii; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); - - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = Real(0.0); - } + if (AtA(ii,ii) > 0.) + { + for (int jj = 0; jj < ii; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } + + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = Real(0.0); + } } - if (AtA(neq-1,neq-1) > 0.) - b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); - else - b(neq-1) = Real(0.0); + if (AtA(neq-1,neq-1) > 0.) { + b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); + } else { + b(neq-1) = Real(0.0); + } for (int ii = neq-2; ii >= 0; ii--) { - if (AtA(ii,ii) > 0.) - { - for (int jj = ii+1; jj < neq; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); - - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = Real(0.0); - } + if (AtA(ii,ii) > 0.) + { + for (int jj = ii+1; jj < neq; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } + + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = Real(0.0); + } } } @@ -193,122 +202,130 @@ void cholsol_for_eb(Array2D& Amatrix, Array1D& b) Array2D AtA; - for (int irow = 0; irow < neq; irow++) - for (int icol = 0; icol < neq; icol++) + for (int irow = 0; irow < neq; irow++) { + for (int icol = 0; icol < neq; icol++) { AtA(irow,icol) = Real(0.0); + } + } for (int irow = 0; irow < 54; irow++) { - AtA(0,0) += Amatrix(irow,0)*Amatrix(irow,0); // e^T e - AtA(0,1) += Amatrix(irow,0)*Amatrix(irow,1); // e^T x - AtA(0,2) += Amatrix(irow,0)*Amatrix(irow,2); // e^T y - AtA(0,3) += Amatrix(irow,0)*Amatrix(irow,3); // e^T z - AtA(0,4) += Amatrix(irow,0)*Amatrix(irow,4); // e^T x^2 - AtA(0,5) += Amatrix(irow,0)*Amatrix(irow,5); // e^T x*y - AtA(0,6) += Amatrix(irow,0)*Amatrix(irow,6); // e^T y^2 - AtA(0,7) += Amatrix(irow,0)*Amatrix(irow,7); // e^T x*z - AtA(0,8) += Amatrix(irow,0)*Amatrix(irow,8); // e^T y*z - AtA(0,9) += Amatrix(irow,0)*Amatrix(irow,9); // e^T z^2 - - AtA(1,1) += Amatrix(irow,1)*Amatrix(irow,1); // x^T x - AtA(1,2) += Amatrix(irow,1)*Amatrix(irow,2); // x^T y - AtA(1,3) += Amatrix(irow,1)*Amatrix(irow,3); // x^T y - AtA(1,4) += Amatrix(irow,1)*Amatrix(irow,4); // x^T (x^2) - AtA(1,5) += Amatrix(irow,1)*Amatrix(irow,5); // x^T (xy) - AtA(1,6) += Amatrix(irow,1)*Amatrix(irow,6); // x^T (y^2) - AtA(1,7) += Amatrix(irow,1)*Amatrix(irow,7); // x^T x*z - AtA(1,8) += Amatrix(irow,1)*Amatrix(irow,8); // x^T y*z - AtA(1,9) += Amatrix(irow,1)*Amatrix(irow,9); // x^T z^2 - - AtA(2,2) += Amatrix(irow,2)*Amatrix(irow,2); // y^T y - AtA(2,3) += Amatrix(irow,2)*Amatrix(irow,3); // y^T z - AtA(2,4) += Amatrix(irow,2)*Amatrix(irow,4); // y^T (x^2) - AtA(2,5) += Amatrix(irow,2)*Amatrix(irow,5); // y^T (xy) - AtA(2,6) += Amatrix(irow,2)*Amatrix(irow,6); // y^T (y^2) - AtA(2,7) += Amatrix(irow,2)*Amatrix(irow,7); // y^T x*z - AtA(2,8) += Amatrix(irow,2)*Amatrix(irow,8); // y^T y*z - AtA(2,9) += Amatrix(irow,2)*Amatrix(irow,9); // y^T z^2 - - AtA(3,3) += Amatrix(irow,3)*Amatrix(irow,3); // z^T z - AtA(3,4) += Amatrix(irow,3)*Amatrix(irow,4); // z^T (x^2) - AtA(3,5) += Amatrix(irow,3)*Amatrix(irow,5); // z^T (xy) - AtA(3,6) += Amatrix(irow,3)*Amatrix(irow,6); // z^T (y^2) - AtA(3,7) += Amatrix(irow,3)*Amatrix(irow,7); // z^T x*z - AtA(3,8) += Amatrix(irow,3)*Amatrix(irow,8); // z^T y*z - AtA(3,9) += Amatrix(irow,3)*Amatrix(irow,9); // z^T z^2 - - AtA(4,4) += Amatrix(irow,4)*Amatrix(irow,4); // (x^2)^T (x^2) - AtA(4,5) += Amatrix(irow,4)*Amatrix(irow,5); // (x^2)^T (xy) - AtA(4,6) += Amatrix(irow,4)*Amatrix(irow,6); // (x^2)^T (y^2) - AtA(4,7) += Amatrix(irow,4)*Amatrix(irow,7); // (x^2)^T x*z - AtA(4,8) += Amatrix(irow,4)*Amatrix(irow,8); // (x^2)^T y*z - AtA(4,9) += Amatrix(irow,4)*Amatrix(irow,9); // (x^2)^T z^2 - - AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (xy)^T (xy) - AtA(5,6) += Amatrix(irow,5)*Amatrix(irow,6); // (xy)^T (y^2) - AtA(5,7) += Amatrix(irow,5)*Amatrix(irow,7); // (xy)^T x*z - AtA(5,8) += Amatrix(irow,5)*Amatrix(irow,8); // (xy)^T y*z - AtA(5,9) += Amatrix(irow,5)*Amatrix(irow,9); // (xy)^T z^2 - - AtA(6,6) += Amatrix(irow,6)*Amatrix(irow,6); // (y^2)^T (y^2) - AtA(6,7) += Amatrix(irow,6)*Amatrix(irow,7); // (y^2)^T x*z - AtA(6,8) += Amatrix(irow,6)*Amatrix(irow,8); // (y^2)^T y*z - AtA(6,9) += Amatrix(irow,6)*Amatrix(irow,9); // (y^2)^T z^2 - - AtA(7,7) += Amatrix(irow,7)*Amatrix(irow,7); // (xz)^T x*z - AtA(7,8) += Amatrix(irow,7)*Amatrix(irow,8); // (xz)^T y*z - AtA(7,9) += Amatrix(irow,7)*Amatrix(irow,9); // (xz)^T z^2 - - AtA(8,8) += Amatrix(irow,8)*Amatrix(irow,8); // (yz)^T y*z - AtA(8,9) += Amatrix(irow,8)*Amatrix(irow,9); // (yz)^T z^2 - - AtA(9,9) += Amatrix(irow,9)*Amatrix(irow,9); // (z^2)^T z^2 + AtA(0,0) += Amatrix(irow,0)*Amatrix(irow,0); // e^T e + AtA(0,1) += Amatrix(irow,0)*Amatrix(irow,1); // e^T x + AtA(0,2) += Amatrix(irow,0)*Amatrix(irow,2); // e^T y + AtA(0,3) += Amatrix(irow,0)*Amatrix(irow,3); // e^T z + AtA(0,4) += Amatrix(irow,0)*Amatrix(irow,4); // e^T x^2 + AtA(0,5) += Amatrix(irow,0)*Amatrix(irow,5); // e^T x*y + AtA(0,6) += Amatrix(irow,0)*Amatrix(irow,6); // e^T y^2 + AtA(0,7) += Amatrix(irow,0)*Amatrix(irow,7); // e^T x*z + AtA(0,8) += Amatrix(irow,0)*Amatrix(irow,8); // e^T y*z + AtA(0,9) += Amatrix(irow,0)*Amatrix(irow,9); // e^T z^2 + + AtA(1,1) += Amatrix(irow,1)*Amatrix(irow,1); // x^T x + AtA(1,2) += Amatrix(irow,1)*Amatrix(irow,2); // x^T y + AtA(1,3) += Amatrix(irow,1)*Amatrix(irow,3); // x^T y + AtA(1,4) += Amatrix(irow,1)*Amatrix(irow,4); // x^T (x^2) + AtA(1,5) += Amatrix(irow,1)*Amatrix(irow,5); // x^T (xy) + AtA(1,6) += Amatrix(irow,1)*Amatrix(irow,6); // x^T (y^2) + AtA(1,7) += Amatrix(irow,1)*Amatrix(irow,7); // x^T x*z + AtA(1,8) += Amatrix(irow,1)*Amatrix(irow,8); // x^T y*z + AtA(1,9) += Amatrix(irow,1)*Amatrix(irow,9); // x^T z^2 + + AtA(2,2) += Amatrix(irow,2)*Amatrix(irow,2); // y^T y + AtA(2,3) += Amatrix(irow,2)*Amatrix(irow,3); // y^T z + AtA(2,4) += Amatrix(irow,2)*Amatrix(irow,4); // y^T (x^2) + AtA(2,5) += Amatrix(irow,2)*Amatrix(irow,5); // y^T (xy) + AtA(2,6) += Amatrix(irow,2)*Amatrix(irow,6); // y^T (y^2) + AtA(2,7) += Amatrix(irow,2)*Amatrix(irow,7); // y^T x*z + AtA(2,8) += Amatrix(irow,2)*Amatrix(irow,8); // y^T y*z + AtA(2,9) += Amatrix(irow,2)*Amatrix(irow,9); // y^T z^2 + + AtA(3,3) += Amatrix(irow,3)*Amatrix(irow,3); // z^T z + AtA(3,4) += Amatrix(irow,3)*Amatrix(irow,4); // z^T (x^2) + AtA(3,5) += Amatrix(irow,3)*Amatrix(irow,5); // z^T (xy) + AtA(3,6) += Amatrix(irow,3)*Amatrix(irow,6); // z^T (y^2) + AtA(3,7) += Amatrix(irow,3)*Amatrix(irow,7); // z^T x*z + AtA(3,8) += Amatrix(irow,3)*Amatrix(irow,8); // z^T y*z + AtA(3,9) += Amatrix(irow,3)*Amatrix(irow,9); // z^T z^2 + + AtA(4,4) += Amatrix(irow,4)*Amatrix(irow,4); // (x^2)^T (x^2) + AtA(4,5) += Amatrix(irow,4)*Amatrix(irow,5); // (x^2)^T (xy) + AtA(4,6) += Amatrix(irow,4)*Amatrix(irow,6); // (x^2)^T (y^2) + AtA(4,7) += Amatrix(irow,4)*Amatrix(irow,7); // (x^2)^T x*z + AtA(4,8) += Amatrix(irow,4)*Amatrix(irow,8); // (x^2)^T y*z + AtA(4,9) += Amatrix(irow,4)*Amatrix(irow,9); // (x^2)^T z^2 + + AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (xy)^T (xy) + AtA(5,6) += Amatrix(irow,5)*Amatrix(irow,6); // (xy)^T (y^2) + AtA(5,7) += Amatrix(irow,5)*Amatrix(irow,7); // (xy)^T x*z + AtA(5,8) += Amatrix(irow,5)*Amatrix(irow,8); // (xy)^T y*z + AtA(5,9) += Amatrix(irow,5)*Amatrix(irow,9); // (xy)^T z^2 + + AtA(6,6) += Amatrix(irow,6)*Amatrix(irow,6); // (y^2)^T (y^2) + AtA(6,7) += Amatrix(irow,6)*Amatrix(irow,7); // (y^2)^T x*z + AtA(6,8) += Amatrix(irow,6)*Amatrix(irow,8); // (y^2)^T y*z + AtA(6,9) += Amatrix(irow,6)*Amatrix(irow,9); // (y^2)^T z^2 + + AtA(7,7) += Amatrix(irow,7)*Amatrix(irow,7); // (xz)^T x*z + AtA(7,8) += Amatrix(irow,7)*Amatrix(irow,8); // (xz)^T y*z + AtA(7,9) += Amatrix(irow,7)*Amatrix(irow,9); // (xz)^T z^2 + + AtA(8,8) += Amatrix(irow,8)*Amatrix(irow,8); // (yz)^T y*z + AtA(8,9) += Amatrix(irow,8)*Amatrix(irow,9); // (yz)^T z^2 + + AtA(9,9) += Amatrix(irow,9)*Amatrix(irow,9); // (z^2)^T z^2 } - for (int irow = 0; irow < neq-1; irow++) - for (int icol = irow+1; icol < neq; icol++) - AtA(icol,irow) = AtA(irow,icol); + for (int irow = 0; irow < neq-1; irow++) { + for (int icol = irow+1; icol < neq; icol++) { + AtA(icol,irow) = AtA(irow,icol); + } + } decomp_chol_np10(AtA); - if (AtA(0,0) > 0.) - b(0) = b(0) / AtA(0,0); - else - b(0) = 0.; + if (AtA(0,0) > 0.) { + b(0) = b(0) / AtA(0,0); + } else { + b(0) = 0.; + } for (int ii = 1; ii < neq; ii++) { - if (AtA(ii,ii) > 0.) - { - for (int jj = 0; jj < ii; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); - - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = Real(0.0); - } + if (AtA(ii,ii) > 0.) + { + for (int jj = 0; jj < ii; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } + + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = Real(0.0); + } } - if (AtA(neq-1,neq-1) > 0.) - b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); - else - b(neq-1) = Real(0.0); + if (AtA(neq-1,neq-1) > 0.) { + b(neq-1) = b(neq-1) / AtA(neq-1,neq-1); + } else { + b(neq-1) = Real(0.0); + } for (int ii = neq-2; ii >= 0; ii--) { - if (AtA(ii,ii) > 0.) - { - for (int jj = ii+1; jj < neq; jj++) - b(ii) = b(ii) - AtA(ii,jj)*b(jj); - - b(ii) = b(ii) / AtA(ii,ii); - } - else - { - b(ii) = Real(0.0); - } + if (AtA(ii,ii) > 0.) + { + for (int jj = ii+1; jj < neq; jj++) { + b(ii) = b(ii) - AtA(ii,jj)*b(jj); + } + + b(ii) = b(ii) / AtA(ii,ii); + } + else + { + b(ii) = Real(0.0); + } } } @@ -328,45 +345,50 @@ Real grad_x_of_phi_on_centroids(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, second 9 are EB centroids - for (int irow = 0; irow < 36; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 36; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } // Columns 0-3: [e x y z] - for (int ii = i-1; ii <= i; ii++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int jj = j-1; jj <= j+1; jj++) // Tangential to face - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-1)); - - Real x_off = static_cast(ii-i) + Real(0.5); - Real y_off = static_cast(jj-j); - Real z_off = static_cast(kk-k); - - Amatrix(a_ind,0) = Real(1.0); - Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0); - Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_xface; - Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_xface; - - if (!flag(ii,jj,kk).isRegular()) + for (int ii = i-1; ii <= i; ii++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face + if (!flag(ii,jj,kk).isCovered()) { - Amatrix(a_ind+18,0) = Real(1.0); - Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0); - Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_xface; - Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_xface; + int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-1)); + + Real x_off = static_cast(ii-i) + Real(0.5); + Real y_off = static_cast(jj-j); + Real z_off = static_cast(kk-k); + + Amatrix(a_ind,0) = Real(1.0); + Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0); + Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_xface; + Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_xface; + + if (!flag(ii,jj,kk).isRegular()) + { + Amatrix(a_ind+18,0) = Real(1.0); + Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0); + Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_xface; + Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_xface; + } } } + } + } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 36; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -374,18 +396,22 @@ Real grad_x_of_phi_on_centroids(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int ii = i-1; ii <= i; ii++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int jj = j-1; jj <= j+1; jj++) // Tangential to face - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-1)); + for (int ii = i-1; ii <= i; ii++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-1)); - rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); + rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) - rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n); - } + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) { + rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n); + } + } + } + } + } } cholsol_np10(Amatrix, rhs); @@ -409,45 +435,50 @@ Real grad_y_of_phi_on_centroids(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, second 9 are EB centroids - for (int irow = 0; irow < 36; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 36; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } // Columns 0-2: [e x y] - for (int jj = j-1; jj <= j; jj++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) // Tangential to face - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-1)); - - Real x_off = static_cast(ii-i); - Real y_off = static_cast(jj-j) + Real(0.5); - Real z_off = static_cast(kk-k); - - Amatrix(a_ind,0) = Real(1.0); - Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_yface; - Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1); - Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_yface; - - if (!flag(ii,jj,kk).isRegular()) + for (int jj = j-1; jj <= j; jj++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face + if (!flag(ii,jj,kk).isCovered()) { - Amatrix(a_ind+18,0) = Real(1.0); - Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_yface; - Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1); - Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_yface; + int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-1)); + + Real x_off = static_cast(ii-i); + Real y_off = static_cast(jj-j) + Real(0.5); + Real z_off = static_cast(kk-k); + + Amatrix(a_ind,0) = Real(1.0); + Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_yface; + Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1); + Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_yface; + + if (!flag(ii,jj,kk).isRegular()) + { + Amatrix(a_ind+18,0) = Real(1.0); + Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_yface; + Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1); + Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_yface; + } } } + } + } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 36; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -455,18 +486,22 @@ Real grad_y_of_phi_on_centroids(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int jj = j-1; jj <= j; jj++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) // Tangential to face - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-1)); + for (int jj = j-1; jj <= j; jj++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-1)); - rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); + rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) - rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n); + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) { + rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n); + } + } } + } + } } cholsol_np10(Amatrix, rhs); @@ -490,49 +525,52 @@ Real grad_z_of_phi_on_centroids(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, second 9 are EB centroids - for (int irow = 0; irow < 36; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 36; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } // Columns 0-3: [e x y z] for (int kk = k-1; kk <= k; kk++) // Normal to face { - for (int jj = j-1; jj <= j+1; jj++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) // Tangential to face - { - if (!flag(ii,jj,kk).isCovered()) + for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) // Tangential to face { - int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); - - Real x_off = static_cast(ii-i); - Real y_off = static_cast(jj-j); - Real z_off = static_cast(kk-k) + Real(0.5); - - Amatrix(a_ind,0) = Real(1.0); - Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_zface; - Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_zface; - Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2); - - if (!flag(ii,jj,kk).isRegular()) + if (!flag(ii,jj,kk).isCovered()) { - Amatrix(a_ind+18,0) = Real(1.0); - Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_zface; - Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_zface; - Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2); + int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); + + Real x_off = static_cast(ii-i); + Real y_off = static_cast(jj-j); + Real z_off = static_cast(kk-k) + Real(0.5); + + Amatrix(a_ind,0) = Real(1.0); + Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_zface; + Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_zface; + Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2); + + if (!flag(ii,jj,kk).isRegular()) + { + Amatrix(a_ind+18,0) = Real(1.0); + Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_zface; + Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_zface; + Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2); + } } } - } + } } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 36; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -540,18 +578,22 @@ Real grad_z_of_phi_on_centroids(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int kk = k-1; kk <= k; kk++) // Normal to face - for (int jj = j-1; jj <= j+1; jj++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) // Tangential to face - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); + for (int kk = k-1; kk <= k; kk++) { // Normal to face + for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); - rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); + rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) - rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n); + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) { + rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n); + } + } } + } + } } cholsol_np10(Amatrix, rhs); @@ -574,47 +616,51 @@ Real grad_eb_of_phi_on_centroids(int i,int j,int k,int n, // Order of column -- first 27 are cell centroids, second 27 are EB centroids - for (int irow = 0; irow < 54; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 54; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } // Columns 0-3: [e x y z] - for (int kk = k-1; kk <= k+1; kk++) - for (int jj = j-1; jj <= j+1; jj++) - for (int ii = i-1; ii <= i+1; ii++) - { - if (!flag(ii,jj,kk).isCovered()) + for (int kk = k-1; kk <= k+1; kk++) { + for (int jj = j-1; jj <= j+1; jj++) { + for (int ii = i-1; ii <= i+1; ii++) { - int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); - - Real x_off = static_cast(ii-i) - bcent(i,j,k,0); - Real y_off = static_cast(jj-j) - bcent(i,j,k,1); - Real z_off = static_cast(kk-k) - bcent(i,j,k,2); - - Amatrix(a_ind,0) = Real(1.0); - Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0); - Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1); - Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2); - - if (!flag(ii,jj,kk).isRegular()) + if (!flag(ii,jj,kk).isCovered()) { - Amatrix(a_ind+27,0) = Real(1.0); - Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0); - Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1); - Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2); + int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); + + Real x_off = static_cast(ii-i) - bcent(i,j,k,0); + Real y_off = static_cast(jj-j) - bcent(i,j,k,1); + Real z_off = static_cast(kk-k) - bcent(i,j,k,2); + + Amatrix(a_ind,0) = Real(1.0); + Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0); + Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1); + Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2); + + if (!flag(ii,jj,kk).isRegular()) + { + Amatrix(a_ind+27,0) = Real(1.0); + Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0); + Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1); + Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2); + } } } - } + } + } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 54; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -622,18 +668,22 @@ Real grad_eb_of_phi_on_centroids(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int kk = k-1; kk <= k+1; kk++) - for (int jj = j-1; jj <= j+1; jj++) - for (int ii = i-1; ii <= i+1; ii++) - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); + for (int kk = k-1; kk <= k+1; kk++) { + for (int jj = j-1; jj <= j+1; jj++) { + for (int ii = i-1; ii <= i+1; ii++) { + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); - rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); + rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n); - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) - rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) { + rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + } + } } + } + } } cholsol_for_eb(Amatrix, rhs); @@ -663,97 +713,108 @@ Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, second 9 are EB centroids - for (int irow = 0; irow < 54; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 54; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } const int im = (i > domhi_x) ? 2 : 1; const int ip = 2 - im; // Columns 0-3: [e x y z] - for (int ii = i-im; ii <= i+ip; ii++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int jj = j-1; jj <= j+1; jj++) // Tangential to face - { - // Don't include corner cells. Could make this even more strict - // by removing the on_?_face restrictions. - if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { - continue; - } - - if ( !phi.contains(ii,jj,kk) ) - continue; - - if (!flag(ii,jj,kk).isCovered()) + for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int jj = j-1; jj <= j+1; jj++) // Tangential to face { - - int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-im)); - - Real x_off = static_cast(ii-i) + Real(0.5); - Real y_off = static_cast(jj-j); - Real z_off = static_cast(kk-k); - - if(on_x_face){ - if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) - continue; - if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) - continue; + // Don't include corner cells. Could make this even more strict + // by removing the on_?_face restrictions. + if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { + continue; } - if(on_y_face){ - if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) - continue; - if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) - continue; + if ( !phi.contains(ii,jj,kk) ) { + continue; } - if(on_z_face){ - if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) - continue; - if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) - continue; - } + if (!flag(ii,jj,kk).isCovered()) + { - Amatrix(a_ind,0) = Real(1.0); - Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0); - Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_xface; - Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_xface; + int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-im)); - // Add in information about the location of the EB. Exclude - // EBs that are outside the domain. - if (flag(ii,jj,kk).isSingleValued() && - domlo_x <= ii && ii <= domhi_x && - domlo_y <= jj && jj <= domhi_y && - domlo_z <= kk && kk <= domhi_z) - { - Amatrix(a_ind+27,0) = Real(1.0); - Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0); - Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_xface; - Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_xface; + Real x_off = static_cast(ii-i) + Real(0.5); + Real y_off = static_cast(jj-j); + Real z_off = static_cast(kk-k); + + if(on_x_face){ + if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) { + continue; + } + if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) { + continue; + } + } + + if(on_y_face){ + if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) { + continue; + } + if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) { + continue; + } + } + + if(on_z_face){ + if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) { + continue; + } + if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) { + continue; + } + } + + Amatrix(a_ind,0) = Real(1.0); + Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0); + Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_xface; + Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_xface; + + // Add in information about the location of the EB. Exclude + // EBs that are outside the domain. + if (flag(ii,jj,kk).isSingleValued() && + domlo_x <= ii && ii <= domhi_x && + domlo_y <= jj && jj <= domhi_y && + domlo_z <= kk && kk <= domhi_z) + { + Amatrix(a_ind+27,0) = Real(1.0); + Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0); + Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_xface; + Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_xface; + } } } - } + } + } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 54; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -761,23 +822,27 @@ Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int ii = i-im; ii <= i+ip; ii++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face - if ( !phi.contains(ii,jj,kk) ) - continue; + for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face + if ( !phi.contains(ii,jj,kk) ) { + continue; + } - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-im)); - Real phi_val = Amatrix(a_ind,0)*phi(ii,jj,kk,n); + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-im)); + Real phi_val = Amatrix(a_ind,0)*phi(ii,jj,kk,n); - rhs(irow) += Amatrix(a_ind,irow)* phi_val; + rhs(irow) += Amatrix(a_ind,irow)* phi_val; - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) - rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); - } + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) { + rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + } + } + } } + } } cholsol_for_eb(Amatrix, rhs); @@ -805,96 +870,107 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, second 9 are EB centroids - for (int irow = 0; irow < 54; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 54; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } const int jm = (j > domhi_y) ? 2 : 1; const int jp = 2 - jm; // Columns 0-2: [e x y] - for (int jj = j-jm; jj <= j+jp; jj++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) // Tangential to face - { - // Don't include corner cells. Could make this even more strict - // by removing the on_?_face restrictions. - if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { - continue; - } - - if ( !phi.contains(ii,jj,kk) ) - continue; - - if (!flag(ii,jj,kk).isCovered()) + for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) // Tangential to face { - int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-jm)); - - Real x_off = static_cast(ii-i); - Real y_off = static_cast(jj-j) + Real(0.5); - Real z_off = static_cast(kk-k); - - if(on_x_face){ - if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) - continue; - if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) - continue; + // Don't include corner cells. Could make this even more strict + // by removing the on_?_face restrictions. + if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { + continue; } - if(on_y_face){ - if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) - continue; - if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) - continue; + if ( !phi.contains(ii,jj,kk) ) { + continue; } - if(on_z_face){ - if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) - continue; - if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) - continue; - } - - Amatrix(a_ind,0) = Real(1.0); - Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_yface; - Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1); - Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_yface; - - // Add in information about the location of the EB. Exclude - // EBs that are outside the domain. - if (flag(ii,jj,kk).isSingleValued() && - domlo_x <= ii && ii <= domhi_x && - domlo_y <= jj && jj <= domhi_y && - domlo_z <= kk && kk <= domhi_z) + if (!flag(ii,jj,kk).isCovered()) { - Amatrix(a_ind+27,0) = Real(1.0); - Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_yface; - Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1); - Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_yface; + int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-jm)); + + Real x_off = static_cast(ii-i); + Real y_off = static_cast(jj-j) + Real(0.5); + Real z_off = static_cast(kk-k); + + if(on_x_face){ + if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) { + continue; + } + if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) { + continue; + } + } + + if(on_y_face){ + if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) { + continue; + } + if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) { + continue; + } + } + + if(on_z_face){ + if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) { + continue; + } + if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) { + continue; + } + } + + Amatrix(a_ind,0) = Real(1.0); + Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_yface; + Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1); + Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_yface; + + // Add in information about the location of the EB. Exclude + // EBs that are outside the domain. + if (flag(ii,jj,kk).isSingleValued() && + domlo_x <= ii && ii <= domhi_x && + domlo_y <= jj && jj <= domhi_y && + domlo_z <= kk && kk <= domhi_z) + { + Amatrix(a_ind+27,0) = Real(1.0); + Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_yface; + Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1); + Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_yface; + } } } } + } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 54; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -902,23 +978,27 @@ Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int jj = j-jm; jj <= j+jp; jj++) // Normal to face - for (int kk = k-1; kk <= k+1; kk++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face - if ( !phi.contains(ii,jj,kk) ) - continue; + for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face + for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face + if ( !phi.contains(ii,jj,kk) ) { + continue; + } - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-jm)); - Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n); + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-jm)); + Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n); - rhs(irow) += Amatrix(a_ind,irow)* phi_val; + rhs(irow) += Amatrix(a_ind,irow)* phi_val; - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) - rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) { + rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + } + } } - } + } + } } cholsol_for_eb(Amatrix, rhs); @@ -946,9 +1026,11 @@ Real grad_z_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Order of column -- first 9 are cell centroids, second 9 are EB centroids - for (int irow = 0; irow < 54; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 54; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } const int km = (k > domhi_z) ? 2 : 1; const int kp = 2 - km; @@ -956,88 +1038,96 @@ Real grad_z_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Columns 0-3: [e x y z] for (int kk = k-km; kk <= k+kp; kk++) // Normal to face { - for (int jj = j-1; jj <= j+1; jj++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) // Tangential to face - { - // Don't include corner cells. Could make this even more strict - // by removing the on_?_face restrictions. - if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { - continue; - } - - if (!phi.contains(ii,jj,kk)) - continue; - - if (!flag(ii,jj,kk).isCovered()) + for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) // Tangential to face { - int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-km)); - - Real x_off = static_cast(ii-i); - Real y_off = static_cast(jj-j); - Real z_off = static_cast(kk-k) + Real(0.5); - - if(on_x_face){ - if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) - continue; - if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) - continue; - } - - if(on_y_face){ - if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) - continue; - if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) - continue; + // Don't include corner cells. Could make this even more strict + // by removing the on_?_face restrictions. + if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { + continue; } - if(on_z_face){ - if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) - continue; - if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) - continue; + if (!phi.contains(ii,jj,kk)) { + continue; } - Amatrix(a_ind,0) = Real(1.0); - Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_zface; - Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_zface ; - Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2); - - // Add in information about the location of the EB. Exclude - // EBs that are outside the domain. - if (flag(ii,jj,kk).isSingleValued() && - domlo_x <= ii && ii <= domhi_x && - domlo_y <= jj && jj <= domhi_y && - domlo_z <= kk && kk <= domhi_z) + if (!flag(ii,jj,kk).isCovered()) { - Amatrix(a_ind+27,0) = Real(1.0); - Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_zface; - Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_zface; - Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2); + int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-km)); + + Real x_off = static_cast(ii-i); + Real y_off = static_cast(jj-j); + Real z_off = static_cast(kk-k) + Real(0.5); + + if(on_x_face){ + if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) { + continue; + } + if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) { + continue; + } + } + + if(on_y_face){ + if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) { + continue; + } + if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) { + continue; + } + } + + if(on_z_face){ + if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) { + continue; + } + if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) { + continue; + } + } + + Amatrix(a_ind,0) = Real(1.0); + Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_zface; + Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_zface ; + Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2); + + // Add in information about the location of the EB. Exclude + // EBs that are outside the domain. + if (flag(ii,jj,kk).isSingleValued() && + domlo_x <= ii && ii <= domhi_x && + domlo_y <= jj && jj <= domhi_y && + domlo_z <= kk && kk <= domhi_z) + { + Amatrix(a_ind+27,0) = Real(1.0); + Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_zface; + Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_zface; + Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2); + } } } - } + } } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 54; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -1045,23 +1135,27 @@ Real grad_z_of_phi_on_centroids_extdir(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int kk = k-km; kk <= k+kp; kk++) // Normal to face - for (int jj = j-1; jj <= j+1; jj++) // Tangential to face - for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face - if ( !phi.contains(ii,jj,kk) ) - continue; + for (int kk = k-km; kk <= k+kp; kk++) { // Normal to face + for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face + for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face + if ( !phi.contains(ii,jj,kk) ) { + continue; + } - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-km)); - Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n); + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-km)); + Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n); - rhs(irow) += Amatrix(a_ind,irow)* phi_val; + rhs(irow) += Amatrix(a_ind,irow)* phi_val; - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) - rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) { + rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + } + } } - } + } + } } cholsol_for_eb(Amatrix, rhs); @@ -1088,37 +1182,40 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, // Order of column -- first 27 are cell centroids, second 27 are EB centroids - for (int irow = 0; irow < 54; irow++) - for (int icol = 0; icol < 10; icol++) + for (int irow = 0; irow < 54; irow++) { + for (int icol = 0; icol < 10; icol++) { Amatrix(irow,icol) = Real(0.0); + } + } // Columns 0-3: [e x y z] - for (int kk = k-1; kk <= k+1; kk++) - for (int jj = j-1; jj <= j+1; jj++) - for (int ii = i-1; ii <= i+1; ii++) - { - // This is likely overkill for EB grads. - // Don't include corner cells. Could make this even more strict - // by removing the on_?_face restrictions. - if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || - ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || - ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || - ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { - continue; - } + for (int kk = k-1; kk <= k+1; kk++) { + for (int jj = j-1; jj <= j+1; jj++) { + for (int ii = i-1; ii <= i+1; ii++) + { + // This is likely overkill for EB grads. + // Don't include corner cells. Could make this even more strict + // by removing the on_?_face restrictions. + if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) || + ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) || + ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) || + ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) { + continue; + } - if ( !phi.contains(ii,jj,kk) ) - continue; + if ( !phi.contains(ii,jj,kk) ) { + continue; + } - if (!flag(ii,jj,kk).isCovered()) + if (!flag(ii,jj,kk).isCovered()) { int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); @@ -1127,24 +1224,30 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, Real z_off = static_cast(kk-k) - bcent(i,j,k,2); if(on_x_face){ - if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) + if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) { continue; - if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) + } + if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) { continue; + } } if(on_y_face){ - if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) + if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) { continue; - if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) + } + if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) { continue; + } } if(on_z_face){ - if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) + if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) { continue; - if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) + } + if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) { continue; + } } Amatrix(a_ind,0) = Real(1.0); @@ -1163,17 +1266,19 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2); } } - } + } + } + } // Columns 4-9 : [x*x x*y y*y x*z y*z z*z] for (int irow = 0; irow < 54; irow++) { - Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); - Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); - Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); - Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); - Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); - Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); + Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1); + Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2); + Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2); + Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3); + Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3); + Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3); } // Make the RHS = A^T v @@ -1181,24 +1286,28 @@ Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n, { rhs(irow) = 0.; - for (int kk = k-1; kk <= k+1; kk++) - for (int jj = j-1; jj <= j+1; jj++) - for (int ii = i-1; ii <= i+1; ii++) { - if ( !phi.contains(ii,jj,kk) ) - continue; + for (int kk = k-1; kk <= k+1; kk++) { + for (int jj = j-1; jj <= j+1; jj++) { + for (int ii = i-1; ii <= i+1; ii++) { + if ( !phi.contains(ii,jj,kk) ) { + continue; + } - if (!flag(ii,jj,kk).isCovered()) - { - int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); + if (!flag(ii,jj,kk).isCovered()) + { + int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1)); - Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n); + Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n); - rhs(irow) += Amatrix(a_ind,irow)* phi_val; + rhs(irow) += Amatrix(a_ind,irow)* phi_val; - if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) - rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) { + rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n); + } + } } - } + } + } } cholsol_for_eb(Amatrix, rhs); diff --git a/Src/EB/AMReX_EB_Redistribution.cpp b/Src/EB/AMReX_EB_Redistribution.cpp index 28fd8b70776..54fcea3889d 100644 --- a/Src/EB/AMReX_EB_Redistribution.cpp +++ b/Src/EB/AMReX_EB_Redistribution.cpp @@ -31,12 +31,14 @@ namespace amrex { const Real* dx = geom.CellSize(); #if (AMREX_SPACEDIM == 2) - if (! amrex::almostEqual(dx[0], dx[1])) + if (! amrex::almostEqual(dx[0], dx[1])) { amrex::Abort("apply_eb_redistribution(): grid spacing must be uniform"); + } #elif (AMREX_SPACEDIM == 3) if( ! amrex::almostEqual(dx[0],dx[1]) || - ! amrex::almostEqual(dx[1],dx[2]) ) + ! amrex::almostEqual(dx[1],dx[2]) ) { amrex::Abort("apply_eb_redistribution(): grid spacing must be uniform"); + } #endif // diff --git a/Src/EB/AMReX_EB_RedistributionApply.cpp b/Src/EB/AMReX_EB_RedistributionApply.cpp index d7c26a69378..d1900a435c9 100644 --- a/Src/EB/AMReX_EB_RedistributionApply.cpp +++ b/Src/EB/AMReX_EB_RedistributionApply.cpp @@ -91,29 +91,30 @@ void ApplyRedistribution ( Box const& bx, int ncomp, Array4 cent_hat_const = cent_hat_fab.const_array(); Box domain_per_grown = lev_geom.Domain(); - AMREX_D_TERM(if (lev_geom.isPeriodic(0)) domain_per_grown.grow(0,1);, - if (lev_geom.isPeriodic(1)) domain_per_grown.grow(1,1);, - if (lev_geom.isPeriodic(2)) domain_per_grown.grow(2,1);); + AMREX_D_TERM(if (lev_geom.isPeriodic(0)) { domain_per_grown.grow(0,1); }, + if (lev_geom.isPeriodic(1)) { domain_per_grown.grow(1,1); }, + if (lev_geom.isPeriodic(2)) { domain_per_grown.grow(2,1); }) // At any external Dirichlet domain boundaries we need to set dUdt_in to 0 // in the cells just outside the domain because those values will be used // in the slope computation in state redistribution. We assume here that // the ext_dir values of U_in itself have already been set. - if (!domain_per_grown.contains(bxg1)) + if (!domain_per_grown.contains(bxg1)) { amrex::ParallelFor(bxg1,ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - if (!domain_per_grown.contains(IntVect(AMREX_D_DECL(i,j,k)))) - dUdt_in(i,j,k,n) = 0.; - }); + { + if (!domain_per_grown.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + dUdt_in(i,j,k,n) = 0.; + } + }); + } amrex::ParallelFor(Box(scratch), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const Real scale = (srd_update_scale) ? srd_update_scale(i,j,k) : Real(1.0); - scratch(i,j,k,n) = U_in(i,j,k,n) + dt * dUdt_in(i,j,k,n) / scale; - } - ); + { + const Real scale = (srd_update_scale) ? srd_update_scale(i,j,k) : Real(1.0); + scratch(i,j,k,n) = U_in(i,j,k,n) + dt * dUdt_in(i,j,k,n) / scale; + }); MakeITracker(bx, AMREX_D_DECL(apx, apy, apz), vfrac, itr, lev_geom, target_volfrac); @@ -127,40 +128,38 @@ void ApplyRedistribution ( Box const& bx, int ncomp, amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + // Only update the values which actually changed -- this makes + // the results insensitive to tiling -- otherwise cells that aren't + // changed but are in a tile on which StateRedistribute gets called + // will have precision-level changes due to adding/subtracting U_in + // and multiplying/dividing by dt. Here we test on whether (i,j,k) + // has at least one neighbor and/or whether (i,j,k) is in the + // neighborhood of another cell -- if either of those is true the + // value may have changed + + if (itr(i,j,k,0) > 0 || nrs(i,j,k) > 1.) { - // Only update the values which actually changed -- this makes - // the results insensitive to tiling -- otherwise cells that aren't - // changed but are in a tile on which StateRedistribute gets called - // will have precision-level changes due to adding/subtracting U_in - // and multiplying/dividing by dt. Here we test on whether (i,j,k) - // has at least one neighbor and/or whether (i,j,k) is in the - // neighborhood of another cell -- if either of those is true the - // value may have changed - - if (itr(i,j,k,0) > 0 || nrs(i,j,k) > 1.) - { - const Real scale = (srd_update_scale) ? srd_update_scale(i,j,k) : Real(1.0); + const Real scale = (srd_update_scale) ? srd_update_scale(i,j,k) : Real(1.0); - dUdt_out(i,j,k,n) = scale * (dUdt_out(i,j,k,n) - U_in(i,j,k,n)) / dt; + dUdt_out(i,j,k,n) = scale * (dUdt_out(i,j,k,n) - U_in(i,j,k,n)) / dt; - } - else - { - dUdt_out(i,j,k,n) = dUdt_in(i,j,k,n); - } } - ); + else + { + dUdt_out(i,j,k,n) = dUdt_in(i,j,k,n); + } + }); } else if (redistribution_type == "NoRedist") { amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - dUdt_out(i,j,k,n) = dUdt_in(i,j,k,n); - } - ); + { + dUdt_out(i,j,k,n) = dUdt_in(i,j,k,n); + }); } else { - amrex::Error("Not a legit redist_type"); + amrex::Error("Not a legit redist_type"); } } @@ -238,7 +237,7 @@ ApplyMLRedistribution ( Box const& bx, int ncomp, ); } else { - amrex::Error("Not a legit redist_type in ApplyML"); + amrex::Error("Not a legit redist_type in ApplyML"); } } diff --git a/Src/EB/AMReX_EB_Slopes_2D_K.H b/Src/EB/AMReX_EB_Slopes_2D_K.H index 7bc346c04a7..06f1e6cba6a 100644 --- a/Src/EB/AMReX_EB_Slopes_2D_K.H +++ b/Src/EB/AMReX_EB_Slopes_2D_K.H @@ -440,95 +440,93 @@ amrex_calc_slopes_extdir_eb (int i, int j, int k, int n, } else { - amrex::Real A[dim_a][AMREX_SPACEDIM]; + amrex::Real A[dim_a][AMREX_SPACEDIM]; - int lc=0; - int kk = 0; - { - for(int jj(-1); jj<=1; jj++) - for(int ii(-1); ii<=1; ii++) - { + int lc=0; + int kk = 0; + + for(int jj(-1); jj<=1; jj++) { + for(int ii(-1); ii<=1; ii++) + { if( flag(i,j,k).isConnected(ii,jj,kk) && ! (ii==0 && jj==0 && kk==0)) + { + bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1); + bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1); + + bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1); + bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1); + + bool klo_test = false; + bool khi_test = false; + + // These are the default values if no physical boundary + A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0); + A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1); + // Do corrections for entire x-face + if (ilo_test) + { + if (!jlo_test && !jhi_test && !klo_test && !khi_test) + { + A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0); + } + A[lc][0] = amrex::Real(-0.5) ; + } else if (ihi_test) { + + if (!jlo_test && !jhi_test && !klo_test && !khi_test) + { + A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0); + } + A[lc][0] = amrex::Real(0.5) ; + } + + // Do corrections for entire y-face + if (jlo_test) { - bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1); - bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1); - - bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1); - bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1); - - bool klo_test = false; - bool khi_test = false; - - // These are the default values if no physical boundary - A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0); - A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1); - // Do corrections for entire x-face - if (ilo_test) - { - if (!jlo_test && !jhi_test && !klo_test && !khi_test) - { - A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0); - } - A[lc][0] = amrex::Real(-0.5) ; - } else if (ihi_test) { - - if (!jlo_test && !jhi_test && !klo_test && !khi_test) - { - A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0); - } - A[lc][0] = amrex::Real(0.5) ; - } - - // Do corrections for entire y-face - if (jlo_test) - { - if (!ilo_test && !ihi_test && !klo_test && !khi_test) - { - A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0); - } - A[lc][1] = amrex::Real(-0.5) ; - - } else if (jhi_test) { - - if (!ilo_test && !ihi_test && !klo_test && !khi_test) - { - A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0); - } - A[lc][1] = amrex::Real(0.5) ; + if (!ilo_test && !ihi_test && !klo_test && !khi_test) + { + A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0); } + A[lc][1] = amrex::Real(-0.5) ; - A[lc][0] -= ccent(i,j,k,0); - A[lc][1] -= ccent(i,j,k,1); + } else if (jhi_test) { - } else { - A[lc][0] = amrex::Real(0.0); - A[lc][1] = amrex::Real(0.0); + if (!ilo_test && !ihi_test && !klo_test && !khi_test) + { + A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0); + } + A[lc][1] = amrex::Real(0.5) ; } - lc++; - } // i,j - } // k - const auto& slopes = amrex_calc_slopes_eb_given_A (i,j,k,n,A,state,flag); - xslope = slopes[0]; - yslope = slopes[1]; + A[lc][0] -= ccent(i,j,k,0); + A[lc][1] -= ccent(i,j,k,1); - // This will over-write the values of xslope and yslope if appropriate - amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,state,vfrac, - edlo_x,edlo_y,edhi_x,edhi_y, - domlo_x,domlo_y,domhi_x,domhi_y,max_order); + } else { + A[lc][0] = amrex::Real(0.0); + A[lc][1] = amrex::Real(0.0); + } + lc++; + }} // i,j - } // end of needs_bndry_stencil + const auto& slopes = amrex_calc_slopes_eb_given_A (i,j,k,n,A,state,flag); + xslope = slopes[0]; + yslope = slopes[1]; - // Zero out slopes outside of an extdir (or hoextrap) boundary - // TODO: is this the right thing to do at a HOEXTRAP boundary?? - if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) || - (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) ) - { - xslope = 0.; - yslope = 0.; + // This will over-write the values of xslope and yslope if appropriate + amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,state,vfrac, + edlo_x,edlo_y,edhi_x,edhi_y, + domlo_x,domlo_y,domhi_x,domhi_y,max_order); + + // Zero out slopes outside of an extdir (or hoextrap) boundary + // TODO: is this the right thing to do at a HOEXTRAP boundary?? + if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) || + (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) ) + { + xslope = 0.; + yslope = 0.; + } + return {xslope,yslope}; } - return {xslope,yslope}; } // amrex_calc_slopes_extdir_eb_grown calculates the slope in each coordinate direction using a @@ -574,94 +572,92 @@ amrex_calc_slopes_extdir_eb_grown (int i, int j, int k, int n, int nx, int ny, } else { - amrex::Real A[dim_a][AMREX_SPACEDIM]; + amrex::Real A[dim_a][AMREX_SPACEDIM]; - int lc=0; - int kk = 0; - { - for(int jj(-ny); jj<=ny; jj++) - for(int ii(-nx); ii<=nx; ii++) - { + int lc=0; + int kk = 0; + + for(int jj(-ny); jj<=ny; jj++) { + for(int ii(-nx); ii<=nx; ii++) + { if ( !flag(i+ii,j+jj,k).isCovered() && !(ii==0 && jj==0 && kk==0) ) - { - bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1); - bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1); + { + bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1); + bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1); - bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1); - bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1); + bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1); + bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1); - bool klo_test = false; - bool khi_test = false; + bool klo_test = false; + bool khi_test = false; - // These are the default values if no physical boundary - A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0); - A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1); - // Do corrections for entire x-face - if (ilo_test) + // These are the default values if no physical boundary + A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0); + A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1); + // Do corrections for entire x-face + if (ilo_test) + { + if (!jlo_test && !jhi_test && !klo_test && !khi_test) { - if (!jlo_test && !jhi_test && !klo_test && !khi_test) - { - A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0); - } - A[lc][0] = amrex::Real(-0.5) ; - } else if (ihi_test) { - - if (!jlo_test && !jhi_test && !klo_test && !khi_test) - { - A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0); - } - A[lc][0] = amrex::Real(0.5) ; + A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0); } + A[lc][0] = amrex::Real(-0.5) ; + } else if (ihi_test) { - // Do corrections for entire y-face - if (jlo_test) + if (!jlo_test && !jhi_test && !klo_test && !khi_test) { - if (!ilo_test && !ihi_test && !klo_test && !khi_test) - { - A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0); - } - A[lc][1] = amrex::Real(-0.5) ; - - } else if (jhi_test) { - - if (!ilo_test && !ihi_test && !klo_test && !khi_test) - { - A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0); - } - A[lc][1] = amrex::Real(0.5) ; - } - - A[lc][0] -= ccent(i,j,k,0); - A[lc][1] -= ccent(i,j,k,1); - - } else { - A[lc][0] = amrex::Real(0.0); - A[lc][1] = amrex::Real(0.0); + A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0); + } + A[lc][0] = amrex::Real(0.5) ; } - lc++; - } // i,j - } // k - const auto& slopes = amrex_calc_slopes_eb_given_A_grown (i,j,k,n,nx,ny,A,state,flag); - xslope = slopes[0]; - yslope = slopes[1]; + // Do corrections for entire y-face + if (jlo_test) + { + if (!ilo_test && !ihi_test && !klo_test && !khi_test) + { + A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0); + } + A[lc][1] = amrex::Real(-0.5) ; - // This will over-write the values of xslope and yslope if appropriate - amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,state,vfrac, - edlo_x,edlo_y,edhi_x,edhi_y, - domlo_x,domlo_y,domhi_x,domhi_y,max_order); + } else if (jhi_test) { - } // end of needs_bndry_stencil + if (!ilo_test && !ihi_test && !klo_test && !khi_test) + { + A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0); + } + A[lc][1] = amrex::Real(0.5) ; + } - // Zero out slopes outside of an extdir (or hoextrap) boundary - // TODO: is this the right thing to do at a HOEXTRAP boundary?? - if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) || - (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) ) - { - xslope = 0.; - yslope = 0.; + A[lc][0] -= ccent(i,j,k,0); + A[lc][1] -= ccent(i,j,k,1); + + } else { + A[lc][0] = amrex::Real(0.0); + A[lc][1] = amrex::Real(0.0); + } + lc++; + }} // i,j + + const auto& slopes = amrex_calc_slopes_eb_given_A_grown (i,j,k,n,nx,ny,A,state,flag); + xslope = slopes[0]; + yslope = slopes[1]; + + // This will over-write the values of xslope and yslope if appropriate + amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,state,vfrac, + edlo_x,edlo_y,edhi_x,edhi_y, + domlo_x,domlo_y,domhi_x,domhi_y,max_order); + + // Zero out slopes outside of an extdir (or hoextrap) boundary + // TODO: is this the right thing to do at a HOEXTRAP boundary?? + if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) || + (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) ) + { + xslope = 0.; + yslope = 0.; + } + return {xslope,yslope}; } - return {xslope,yslope}; } AMREX_GPU_DEVICE AMREX_FORCE_INLINE @@ -705,8 +701,8 @@ amrex_calc_alpha_limiter(int i, int j, int k, int n, for(int ii(-1); ii<=1; ii++){ if( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0)) { - if ((ii==-1 || ii==1) && jj==0) cuts_x++; - if ((jj==-1 || jj==1) && ii==0) cuts_y++; + if ((ii==-1 || ii==1) && jj==0) { cuts_x++; } + if ((jj==-1 || jj==1) && ii==0) { cuts_y++; } } } } @@ -716,101 +712,101 @@ amrex_calc_alpha_limiter(int i, int j, int k, int n, //Reconstruct values at the face centroids and compute the limiter if(flag(i,j,k).isConnected(0,1,0)) { - amrex::Real xf = fcy(i,j+1,k,0); // local (x,z) of centroid of y-face we are extrapolating to + amrex::Real xf = fcy(i,j+1,k,0); // local (x,z) of centroid of y-face we are extrapolating to - amrex::Real delta_x = xf - xc; - amrex::Real delta_y = amrex::Real(0.5) - yc; + amrex::Real delta_x = xf - xc; + amrex::Real delta_y = amrex::Real(0.5) - yc; - amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + delta_y * slopes[1]; + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + delta_y * slopes[1]; - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); - // Compute max and min values in a 3x2 stencil - for(int jj(0); jj<=1; jj++){ - for(int ii(-1); ii<=1; ii++){ - if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); + // Compute max and min values in a 3x2 stencil + for(int jj(0); jj<=1; jj++){ + for(int ii(-1); ii<=1; ii++){ + if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } } - } - } + } + } - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if (flag(i,j,k).isConnected(0,-1,0)){ - amrex::Real xf = fcy(i,j,k,0); // local (x,z) of centroid of y-face we are extrapolating to + amrex::Real xf = fcy(i,j,k,0); // local (x,z) of centroid of y-face we are extrapolating to - amrex::Real delta_x = xf - xc; - amrex::Real delta_y = amrex::Real(0.5) + yc; + amrex::Real delta_x = xf - xc; + amrex::Real delta_y = amrex::Real(0.5) + yc; - amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] - delta_y * slopes[1]; + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] - delta_y * slopes[1]; - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); - // Compute max and min values in a 3x2 stencil - for(int jj(-1); jj<=0; jj++){ - for(int ii(-1); ii<=1; ii++){ - if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); + // Compute max and min values in a 3x2 stencil + for(int jj(-1); jj<=0; jj++){ + for(int ii(-1); ii<=1; ii++){ + if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } } - } - } + } + } - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if(flag(i,j,k).isConnected(1,0,0)) { - amrex::Real yf = fcx(i+1,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to + amrex::Real yf = fcx(i+1,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to - amrex::Real delta_x = amrex::Real(0.5) - xc; - amrex::Real delta_y = yf - yc; + amrex::Real delta_x = amrex::Real(0.5) - xc; + amrex::Real delta_y = yf - yc; - amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + delta_y * slopes[1]; + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + delta_y * slopes[1]; - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); - for(int jj(-1); jj<=1; jj++){ - for(int ii(0); ii<=1; ii++){ - if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); + for(int jj(-1); jj<=1; jj++){ + for(int ii(0); ii<=1; ii++){ + if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } } - } - } + } + } - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if(flag(i,j,k).isConnected(-1,0,0)) { - amrex::Real yf = fcx(i,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to + amrex::Real yf = fcx(i,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to - amrex::Real delta_x = amrex::Real(0.5) + xc; - amrex::Real delta_y = yf - yc; + amrex::Real delta_x = amrex::Real(0.5) + xc; + amrex::Real delta_y = yf - yc; - amrex::Real q_hat = state(i,j,k,n) - delta_x * slopes[0] + delta_y * slopes[1]; + amrex::Real q_hat = state(i,j,k,n) - delta_x * slopes[0] + delta_y * slopes[1]; - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); - for(int jj(-1); jj<=1; jj++){ - for(int ii(-1); ii<=0; ii++){ - if( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0)) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); + for(int jj(-1); jj<=1; jj++){ + for(int ii(-1); ii<=0; ii++){ + if( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0)) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } } - } - } - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + } + } + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } amrex::Real xalpha = alpha; amrex::Real yalpha = alpha; //Zeroing out the slopes in the direction where a covered face exists. - if (cuts_x<2) xalpha = 0; - if (cuts_y<2) yalpha = 0; + if (cuts_x<2) { xalpha = 0; } + if (cuts_y<2) { yalpha = 0; } return {xalpha,yalpha}; } @@ -838,11 +834,13 @@ amrex_lim_slopes_eb (int i, int j, int k, int n, // Setting limiter to 1 for stencils that just consists of non-EB cells because // amrex_calc_slopes_eb routine will call the slope routine for non-EB stencils that has already a limiter - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) - alpha_lim[0] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) { + alpha_lim[0] = 1.0; + } - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) - alpha_lim[1] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) { + alpha_lim[1] = 1.0; + } return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1]}; } @@ -873,11 +871,13 @@ amrex_lim_slopes_extdir_eb (int i, int j, int k, int n, // Setting limiter to 1 for stencils that just consists of non-EB cells because // amrex_calc_slopes_extdir_eb routine will call the slope routine for non-EB stencils that has already a limiter - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) { alpha_lim[0] = 1.0; + } - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) - alpha_lim[1] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) { + alpha_lim[1] = 1.0; + } return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1]}; } diff --git a/Src/EB/AMReX_EB_Slopes_3D_K.H b/Src/EB/AMReX_EB_Slopes_3D_K.H index a7e91c7b905..3de0e92287f 100644 --- a/Src/EB/AMReX_EB_Slopes_3D_K.H +++ b/Src/EB/AMReX_EB_Slopes_3D_K.H @@ -27,30 +27,26 @@ amrex_calc_slopes_eb_given_A (int i, int j, int k, int n, amrex::Real du[dim_a]; int ll=0; - for(int kk(-1); kk<=1; kk++) - { - for(int jj(-1); jj<=1; jj++){ - for(int ii(-1); ii<=1; ii++){ - - if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) - { - du[ll] = state(i+ii,j+jj,k+kk,n) - state(i,j,k,n); - } else { - du[ll] = amrex::Real(0.0); - } - ll++; - } + for(int kk(-1); kk<=1; kk++) { + for(int jj(-1); jj<=1; jj++){ + for(int ii(-1); ii<=1; ii++){ + if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) + { + du[ll] = state(i+ii,j+jj,k+kk,n) - state(i,j,k,n); + } else { + du[ll] = amrex::Real(0.0); } - } + ll++; + }}} amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM]; amrex::Real Atb[AMREX_SPACEDIM]; for(int jj(0); jj q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); - } - } - } - } - - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + amrex::Real xf = fcy(i,j+1,k,0); // local (x,z) of centroid of y-face we are extrapolating to + amrex::Real zf = fcy(i,j+1,k,1); // local (x,z) of centroid of y-face we are extrapolating to + + amrex::Real delta_x = xf - xc; + amrex::Real delta_y = amrex::Real(0.5) - yc; + amrex::Real delta_z = zf - zc; + + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + + delta_y * slopes[1] + + delta_z * slopes[2]; + + // Compute max and min values in a 3x2x3 stencil + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); + for(int kk(-1); kk<=1; kk++) { + for(int jj(0); jj<=1; jj++){ + for(int ii(-1); ii<=1; ii++){ + if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } + } + } + } + } + + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if (flag(i,j,k).isConnected(0,-1,0)){ - amrex::Real xf = fcy(i,j,k,0); // local (x,z) of centroid of y-face we are extrapolating to - amrex::Real zf = fcy(i,j,k,1); // local (x,z) of centroid of y-face we are extrapolating to - - amrex::Real delta_x = xf - xc; - amrex::Real delta_y = amrex::Real(0.5) + yc; - amrex::Real delta_z = zf - zc; - - amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] - - delta_y * slopes[1] - + delta_z * slopes[2]; - - // Compute max and min values in a 3x2x3 stencil - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); - for(int kk(-1); kk<=1; kk++) { - for(int jj(-1); jj<=0; jj++){ - for(int ii(-1); ii<=1; ii++){ - if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); - } - } - } - } - - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + amrex::Real xf = fcy(i,j,k,0); // local (x,z) of centroid of y-face we are extrapolating to + amrex::Real zf = fcy(i,j,k,1); // local (x,z) of centroid of y-face we are extrapolating to + + amrex::Real delta_x = xf - xc; + amrex::Real delta_y = amrex::Real(0.5) + yc; + amrex::Real delta_z = zf - zc; + + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + - delta_y * slopes[1] + + delta_z * slopes[2]; + + // Compute max and min values in a 3x2x3 stencil + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); + for(int kk(-1); kk<=1; kk++) { + for(int jj(-1); jj<=0; jj++){ + for(int ii(-1); ii<=1; ii++){ + if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } + } + } + } + } + + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if (flag(i,j,k).isConnected(1,0,0)) { - amrex::Real yf = fcx(i+1,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to - amrex::Real zf = fcx(i+1,j,k,1); // local (y,z) of centroid of x-face we are extrapolating to - - amrex::Real delta_x = amrex::Real(0.5) - xc; - amrex::Real delta_y = yf - yc; - amrex::Real delta_z = zf - zc; - - amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] - + delta_y * slopes[1] - + delta_z * slopes[2]; - - // Compute max and min values in a 2x3x3 stencil - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); - for(int kk(-1); kk<=1; kk++) { - for(int jj(-1); jj<=1; jj++){ - for(int ii(0); ii<=1; ii++){ - if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); - } - } - } - } - - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + amrex::Real yf = fcx(i+1,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to + amrex::Real zf = fcx(i+1,j,k,1); // local (y,z) of centroid of x-face we are extrapolating to + + amrex::Real delta_x = amrex::Real(0.5) - xc; + amrex::Real delta_y = yf - yc; + amrex::Real delta_z = zf - zc; + + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + + delta_y * slopes[1] + + delta_z * slopes[2]; + + // Compute max and min values in a 2x3x3 stencil + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); + for(int kk(-1); kk<=1; kk++) { + for(int jj(-1); jj<=1; jj++){ + for(int ii(0); ii<=1; ii++){ + if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } + } + } + } + } + + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if (flag(i,j,k).isConnected(-1,0,0)) { - amrex::Real yf = fcx(i,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to - amrex::Real zf = fcx(i,j,k,1); // local (y,z) of centroid of x-face we are extrapolating to - - amrex::Real delta_x = amrex::Real(0.5) + xc; - amrex::Real delta_y = yf - yc; - amrex::Real delta_z = zf - zc; - - amrex::Real q_hat = state(i,j,k,n) - delta_x * slopes[0] - + delta_y * slopes[1] - + delta_z * slopes[2]; - - // Compute max and min values in a 3x2x3 stencil - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); - for(int kk(-1); kk<=1; kk++) { - for(int jj(-1); jj<=1; jj++){ - for(int ii(-1); ii<=0; ii++){ - if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); - } - } - } - } - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + amrex::Real yf = fcx(i,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to + amrex::Real zf = fcx(i,j,k,1); // local (y,z) of centroid of x-face we are extrapolating to + + amrex::Real delta_x = amrex::Real(0.5) + xc; + amrex::Real delta_y = yf - yc; + amrex::Real delta_z = zf - zc; + + amrex::Real q_hat = state(i,j,k,n) - delta_x * slopes[0] + + delta_y * slopes[1] + + delta_z * slopes[2]; + + // Compute max and min values in a 3x2x3 stencil + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); + for(int kk(-1); kk<=1; kk++) { + for(int jj(-1); jj<=1; jj++){ + for(int ii(-1); ii<=0; ii++){ + if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } + } + } + } + } + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if (flag(i,j,k).isConnected(0,0,1)) { - amrex::Real xf = fcz(i,j,k+1,0); // local (x,y) of centroid of z-face we are extrapolating to - amrex::Real yf = fcz(i,j,k+1,1); // local (x,y) of centroid of z-face we are extrapolating to - - amrex::Real delta_x = xf - xc; - amrex::Real delta_y = yf - yc; - amrex::Real delta_z = amrex::Real(0.5) - zc; - - amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] - + delta_y * slopes[1] - + delta_z * slopes[2]; - - // Compute max and min values in a 3x3x2 stencil - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); - for(int kk(0); kk<=1; kk++) { - for(int jj(-1); jj<=1; jj++){ - for(int ii(-1); ii<=1; ii++){ - if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); - } - } - } - } - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + amrex::Real xf = fcz(i,j,k+1,0); // local (x,y) of centroid of z-face we are extrapolating to + amrex::Real yf = fcz(i,j,k+1,1); // local (x,y) of centroid of z-face we are extrapolating to + + amrex::Real delta_x = xf - xc; + amrex::Real delta_y = yf - yc; + amrex::Real delta_z = amrex::Real(0.5) - zc; + + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + + delta_y * slopes[1] + + delta_z * slopes[2]; + + // Compute max and min values in a 3x3x2 stencil + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); + for(int kk(0); kk<=1; kk++) { + for(int jj(-1); jj<=1; jj++){ + for(int ii(-1); ii<=1; ii++){ + if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } + } + } + } + } + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } if (flag(i,j,k).isConnected(0,0,-1)) { - amrex::Real xf = fcz(i,j,k,0); // local (x,y) of centroid of z-face we are extrapolating to - amrex::Real yf = fcz(i,j,k,1); // local (x,y) of centroid of z-face we are extrapolating to - - amrex::Real delta_x = xf - xc; - amrex::Real delta_y = yf - yc; - amrex::Real delta_z = amrex::Real(0.5) + zc; - - amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] - + delta_y * slopes[1] - - delta_z * slopes[2]; - - // Compute max and min values in a 3x2x3 stencil - amrex::Real q_min = state(i,j,k,n); - amrex::Real q_max = state(i,j,k,n); - for(int kk(-1); kk<=0; kk++) { - for(int jj(-1); jj<=1; jj++){ - for(int ii(-1); ii<=1; ii++){ - if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { - if (state(i+ii,j+jj,k+kk,n) > q_max) q_max = state(i+ii,j+jj,k+kk,n); - if (state(i+ii,j+jj,k+kk,n) < q_min) q_min = state(i+ii,j+jj,k+kk,n); - } - } - } - } - alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); + amrex::Real xf = fcz(i,j,k,0); // local (x,y) of centroid of z-face we are extrapolating to + amrex::Real yf = fcz(i,j,k,1); // local (x,y) of centroid of z-face we are extrapolating to + + amrex::Real delta_x = xf - xc; + amrex::Real delta_y = yf - yc; + amrex::Real delta_z = amrex::Real(0.5) + zc; + + amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + + delta_y * slopes[1] + - delta_z * slopes[2]; + + // Compute max and min values in a 3x2x3 stencil + amrex::Real q_min = state(i,j,k,n); + amrex::Real q_max = state(i,j,k,n); + for(int kk(-1); kk<=0; kk++) { + for(int jj(-1); jj<=1; jj++){ + for(int ii(-1); ii<=1; ii++){ + if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) { + if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); } + if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); } + } + } + } + } + alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha); } amrex::Real xalpha = alpha; @@ -1037,9 +1029,9 @@ amrex_calc_alpha_limiter(int i, int j, int k, int n, amrex::Real zalpha = alpha; //Zeroing out the slopes in the direction where a covered face exists. - if (cuts_x<2) xalpha = 0; - if (cuts_y<2) yalpha = 0; - if (cuts_z<2) zalpha = 0; + if (cuts_x<2) { xalpha = 0; } + if (cuts_y<2) { yalpha = 0; } + if (cuts_z<2) { zalpha = 0; } return {xalpha,yalpha,zalpha}; } @@ -1068,14 +1060,17 @@ amrex_lim_slopes_eb (int i, int j, int k, int n, // Setting limiter to 1 for stencils that just consists of non-EB cells because // amrex_calc_slopes_eb routine will call the slope routine for non-EB stencils that has already a limiter - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) - alpha_lim[0] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) { + alpha_lim[0] = 1.0; + } - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) - alpha_lim[1] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) { + alpha_lim[1] = 1.0; + } - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) - alpha_lim[2] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) { + alpha_lim[2] = 1.0; + } return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1],alpha_lim[2]*slopes[2]}; } @@ -1109,14 +1104,17 @@ amrex_lim_slopes_extdir_eb (int i, int j, int k, int n, // Setting limiter to 1 for stencils that just consists of non-EB cells because // amrex_calc_slopes_extdir_eb routine will call the slope routine for non-EB stencils that has already a limiter - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) { alpha_lim[0] = 1.0; + } - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) - alpha_lim[1] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) { + alpha_lim[1] = 1.0; + } - if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) - alpha_lim[2] = 1.0; + if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) { + alpha_lim[2] = 1.0; + } return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1],alpha_lim[2]*slopes[2]}; } diff --git a/Src/EB/AMReX_EB_StateRedistItracker.cpp b/Src/EB/AMReX_EB_StateRedistItracker.cpp index 50c2360161b..538a917c72a 100644 --- a/Src/EB/AMReX_EB_StateRedistItracker.cpp +++ b/Src/EB/AMReX_EB_StateRedistItracker.cpp @@ -53,8 +53,8 @@ MakeITracker ( Box const& bx, }); Box domain_per_grown = domain; - if (is_periodic_x) domain_per_grown.grow(0,4); - if (is_periodic_y) domain_per_grown.grow(1,4); + if (is_periodic_x) { domain_per_grown.grow(0,4); } + if (is_periodic_y) { domain_per_grown.grow(1,4); } Box const& bxg4 = amrex::grow(bx,4); Box bx_per_g4= domain_per_grown & bxg4; @@ -78,16 +78,18 @@ MakeITracker ( Box const& bx, // As a first pass, choose just based on the normal if (std::abs(nx) > std::abs(ny)) { - if (nx > 0) + if (nx > 0) { itracker(i,j,k,1) = 5; - else + } else { itracker(i,j,k,1) = 4; + } } else { - if (ny > 0) + if (ny > 0) { itracker(i,j,k,1) = 7; - else + } else { itracker(i,j,k,1) = 2; + } } bool xdir_mns_ok = (is_periodic_x || (i > domain.smallEnd(0))); @@ -115,8 +117,9 @@ MakeITracker ( Box const& bx, int joff = jmap[itracker(i,j,k,1)]; // Sanity check - if (vfrac(i+ioff,j+joff,k) == 0.) + if (vfrac(i+ioff,j+joff,k) == 0.) { amrex::Abort(" Trying to merge with covered cell"); + } Real sum_vol = vfrac(i,j,k) + vfrac(i+ioff,j+joff,k); @@ -184,14 +187,15 @@ MakeITracker ( Box const& bx, ioff = imap[itracker(i,j,k,1)] + imap[itracker(i,j,k,2)]; joff = jmap[itracker(i,j,k,1)] + jmap[itracker(i,j,k,2)]; - if (ioff > 0 && joff > 0) + if (ioff > 0 && joff > 0) { itracker(i,j,k,3) = 8; - else if (ioff < 0 && joff > 0) + } else if (ioff < 0 && joff > 0) { itracker(i,j,k,3) = 6; - else if (ioff > 0 && joff < 0) + } else if (ioff > 0 && joff < 0) { itracker(i,j,k,3) = 3; - else + } else { itracker(i,j,k,3) = 1; + } // (i,j) merges with at least three cells now itracker(i,j,k,0) += 1; @@ -270,10 +274,10 @@ MakeITracker ( Box const& bx, }); Box domain_per_grown = domain; - if (is_periodic_x) domain_per_grown.grow(0,4); - if (is_periodic_y) domain_per_grown.grow(1,4); + if (is_periodic_x) { domain_per_grown.grow(0,4); } + if (is_periodic_y) { domain_per_grown.grow(1,4); } #if (AMREX_SPACEDIM == 3) - if (is_periodic_z) domain_per_grown.grow(2,4); + if (is_periodic_z) { domain_per_grown.grow(2,4); } #endif Box const& bxg4 = amrex::grow(bx,4); @@ -312,54 +316,59 @@ MakeITracker ( Box const& bx, if ( (std::abs(nx) > std::abs(ny)) && (std::abs(nx) > std::abs(nz)) ) { - if (nx > 0) + if (nx > 0) { itracker(i,j,k,1) = 5; - else + } else { itracker(i,j,k,1) = 4; + } // y-component of normal is greatest } else if ( (std::abs(ny) >= std::abs(nx)) && (std::abs(ny) > std::abs(nz)) ) { - if (ny > 0) + if (ny > 0) { itracker(i,j,k,1) = 7; - else - + } else { itracker(i,j,k,1) = 2; + } // z-component of normal is greatest } else { - if (nz > 0) + if (nz > 0) { itracker(i,j,k,1) = 22; - else + } else { itracker(i,j,k,1) = 13; + } } // Override above logic if trying to reach outside a domain boundary (and non-periodic) if ( (!xdir_mns_ok && (itracker(i,j,k,1) == 4)) || (!xdir_pls_ok && (itracker(i,j,k,1) == 5)) ) { - if ( (std::abs(ny) > std::abs(nz)) ) + if ( (std::abs(ny) > std::abs(nz)) ) { itracker(i,j,k,1) = (ny > 0) ? 7 : 2; - else + } else { itracker(i,j,k,1) = (nz > 0) ? 22 : 13; + } } if ( (!ydir_mns_ok && (itracker(i,j,k,1) == 2)) || (!ydir_pls_ok && (itracker(i,j,k,1) == 7)) ) { - if ( (std::abs(nx) > std::abs(nz)) ) + if ( (std::abs(nx) > std::abs(nz)) ) { itracker(i,j,k,1) = (nx > 0) ? 5 : 4; - else + } else { itracker(i,j,k,1) = (nz > 0) ? 22 : 13; + } } if ( (!zdir_mns_ok && (itracker(i,j,k,1) == 13)) || (!zdir_pls_ok && (itracker(i,j,k,1) == 22)) ) { - if ( (std::abs(nx) > std::abs(ny)) ) + if ( (std::abs(nx) > std::abs(ny)) ) { itracker(i,j,k,1) = (nx > 0) ? 5 : 4; - else + } else { itracker(i,j,k,1) = (ny > 0) ? 7 : 2; + } } // (i,j,k) merges with at least one cell now @@ -466,36 +475,39 @@ MakeITracker ( Box const& bx, // Both nbors are in the koff=0 plane if (koff == 0) { - if (ioff > 0 && joff > 0) + if (ioff > 0 && joff > 0) { itracker(i,j,k,3) = 8; - else if (ioff < 0 && joff > 0) + } else if (ioff < 0 && joff > 0) { itracker(i,j,k,3) = 6; - else if (ioff > 0 && joff < 0) + } else if (ioff > 0 && joff < 0) { itracker(i,j,k,3) = 3; - else + } else { itracker(i,j,k,3) = 1; + } // Both nbors are in the joff=0 plane } else if (joff == 0) { - if (ioff > 0 && koff > 0) + if (ioff > 0 && koff > 0) { itracker(i,j,k,3) = 23; - else if (ioff < 0 && koff > 0) + } else if (ioff < 0 && koff > 0) { itracker(i,j,k,3) = 21; - else if (ioff > 0 && koff < 0) + } else if (ioff > 0 && koff < 0) { itracker(i,j,k,3) = 14; - else + } else { itracker(i,j,k,3) = 12; + } // Both nbors are in the ioff=0 plane } else { - if (joff > 0 && koff > 0) + if (joff > 0 && koff > 0) { itracker(i,j,k,3) = 25; - else if (joff < 0 && koff > 0) + } else if (joff < 0 && koff > 0) { itracker(i,j,k,3) = 19; - else if (joff > 0 && koff < 0) + } else if (joff > 0 && koff < 0) { itracker(i,j,k,3) = 16; - else + } else { itracker(i,j,k,3) = 10; + } } // (i,j,k) merges with at least three cells now @@ -539,14 +551,16 @@ MakeITracker ( Box const& bx, { itracker(i,j,k,4) = 22; - if (ioff > 0) + if (ioff > 0) { itracker(i,j,k,5) = 23; - else + } else { itracker(i,j,k,5) = 21; - if (joff > 0) + } + if (joff > 0) { itracker(i,j,k,6) = 25; - else + } else { itracker(i,j,k,6) = 19; + } if (ioff > 0 && joff > 0) { itracker(i,j,k,7) = 26; @@ -561,14 +575,16 @@ MakeITracker ( Box const& bx, itracker(i,j,k,4) = 13; - if (ioff > 0) + if (ioff > 0) { itracker(i,j,k,5) = 14; - else + } else { itracker(i,j,k,5) = 12; - if (joff > 0) + } + if (joff > 0) { itracker(i,j,k,6) = 16; - else + } else { itracker(i,j,k,6) = 10; + } if (ioff > 0 && joff > 0) { itracker(i,j,k,7) = 17; @@ -585,14 +601,16 @@ MakeITracker ( Box const& bx, { itracker(i,j,k,4) = 7; - if (ioff > 0) + if (ioff > 0) { itracker(i,j,k,5) = 8; - else + } else { itracker(i,j,k,5) = 6; - if (koff > 0) + } + if (koff > 0) { itracker(i,j,k,6) = 25; - else + } else { itracker(i,j,k,6) = 16; + } if (ioff > 0 && koff > 0) { itracker(i,j,k,7) = 26; @@ -608,14 +626,16 @@ MakeITracker ( Box const& bx, itracker(i,j,k,4) = 2; - if (ioff > 0) + if (ioff > 0) { itracker(i,j,k,5) = 3; - else + } else { itracker(i,j,k,5) = 1; - if (koff > 0) + } + if (koff > 0) { itracker(i,j,k,6) = 19; - else + } else { itracker(i,j,k,6) = 10; + } if (ioff > 0 && koff > 0) { itracker(i,j,k,7) = 20; @@ -633,14 +653,16 @@ MakeITracker ( Box const& bx, { itracker(i,j,k,4) = 5; - if (joff > 0) + if (joff > 0) { itracker(i,j,k,5) = 8; - else + } else { itracker(i,j,k,5) = 3; - if (koff > 0) + } + if (koff > 0) { itracker(i,j,k,6) = 23; - else + } else { itracker(i,j,k,6) = 14; + } if (joff > 0 && koff > 0) { itracker(i,j,k,7) = 26; @@ -655,14 +677,16 @@ MakeITracker ( Box const& bx, itracker(i,j,k,4) = 4; - if (joff > 0) + if (joff > 0) { itracker(i,j,k,5) = 6; - else + } else { itracker(i,j,k,5) = 1; - if (koff > 0) + } + if (koff > 0) { itracker(i,j,k,6) = 21; - else + } else { itracker(i,j,k,6) = 12; + } if (joff > 0 && koff > 0) { itracker(i,j,k,7) = 24; diff --git a/Src/EB/AMReX_EB_StateRedistSlopeLimiter_K.H b/Src/EB/AMReX_EB_StateRedistSlopeLimiter_K.H index 7ce6c5cdd0a..90ce4a07016 100644 --- a/Src/EB/AMReX_EB_StateRedistSlopeLimiter_K.H +++ b/Src/EB/AMReX_EB_StateRedistSlopeLimiter_K.H @@ -74,10 +74,10 @@ amrex_calc_centroid_limiter(int i, int j, int k, int n, { Real new_lim = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n)); - if (amrex::Math::abs(delta_x) > epsilon) xalpha = amrex::min(xalpha,new_lim); - if (amrex::Math::abs(delta_y) > epsilon) yalpha = amrex::min(yalpha,new_lim); + if (amrex::Math::abs(delta_x) > epsilon) { xalpha = amrex::min(xalpha,new_lim); } + if (amrex::Math::abs(delta_y) > epsilon) { yalpha = amrex::min(yalpha,new_lim); } #if (AMREX_SPACEDIM == 3) - if (amrex::Math::abs(delta_z) > epsilon) zalpha = amrex::min(zalpha,new_lim); + if (amrex::Math::abs(delta_z) > epsilon) { zalpha = amrex::min(zalpha,new_lim); } #endif } } diff --git a/Src/EB/AMReX_EB_StateRedistUtils.cpp b/Src/EB/AMReX_EB_StateRedistUtils.cpp index 80539761516..6c819beb860 100644 --- a/Src/EB/AMReX_EB_StateRedistUtils.cpp +++ b/Src/EB/AMReX_EB_StateRedistUtils.cpp @@ -62,10 +62,10 @@ MakeStateRedistUtils ( Box const& bx, const Box domain = lev_geom.Domain(); Box domain_per_grown = domain; - if (is_periodic_x) domain_per_grown.grow(0,2); - if (is_periodic_y) domain_per_grown.grow(1,2); + if (is_periodic_x) { domain_per_grown.grow(0,2); } + if (is_periodic_y) { domain_per_grown.grow(1,2); } #if (AMREX_SPACEDIM == 3) - if (is_periodic_z) domain_per_grown.grow(2,2); + if (is_periodic_z) { domain_per_grown.grow(2,2); } #endif amrex::ParallelFor(bxg3, @@ -114,8 +114,9 @@ MakeStateRedistUtils ( Box const& bx, vol_of_nbors += vfrac(r,s,t); } - if (itracker(i,j,k,0) > 0) + if (itracker(i,j,k,0) > 0) { alpha(i,j,k,1) = (target_vol - vfrac(i,j,k)) / vol_of_nbors; + } } else { nbhd_vol(i,j,k) = 0.; diff --git a/Src/EB/AMReX_EB_StateRedistribute.cpp b/Src/EB/AMReX_EB_StateRedistribute.cpp index af6c255bf63..023c8fe0716 100644 --- a/Src/EB/AMReX_EB_StateRedistribute.cpp +++ b/Src/EB/AMReX_EB_StateRedistribute.cpp @@ -78,10 +78,10 @@ StateRedistribute ( Box const& bx, int ncomp, Box const& bxg3 = amrex::grow(bx,3); Box domain_per_grown = domain; - if (is_periodic_x) domain_per_grown.grow(0,2); - if (is_periodic_y) domain_per_grown.grow(1,2); + if (is_periodic_x) { domain_per_grown.grow(0,2); } + if (is_periodic_y) { domain_per_grown.grow(1,2); } #if (AMREX_SPACEDIM == 3) - if (is_periodic_z) domain_per_grown.grow(2,2); + if (is_periodic_z) { domain_per_grown.grow(2,2); } #endif // Solution at the centroid of my nbhd @@ -96,15 +96,17 @@ StateRedistribute ( Box const& bx, int ncomp, amrex::ParallelFor(bxg3, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - for (int n = 0; n < ncomp; n++) + for (int n = 0; n < ncomp; n++) { soln_hat(i,j,k,n) = U_in(i,j,k,n); + } if (vfrac(i,j,k) > 0.0 && bxg2.contains(IntVect(AMREX_D_DECL(i,j,k))) && domain_per_grown.contains(IntVect(AMREX_D_DECL(i,j,k)))) { // Start with U_in(i,j,k) itself - for (int n = 0; n < ncomp; n++) + for (int n = 0; n < ncomp; n++) { soln_hat(i,j,k,n) = U_in(i,j,k,n) * alpha(i,j,k,0) * vfrac(i,j,k); + } // This loops over the neighbors of (i,j,k), and doesn't include (i,j,k) itself for (int i_nbor = 1; i_nbor <= itracker(i,j,k,0); i_nbor++) @@ -115,12 +117,14 @@ StateRedistribute ( Box const& bx, int ncomp, if (domain_per_grown.contains(IntVect(AMREX_D_DECL(r,s,t)))) { - for (int n = 0; n < ncomp; n++) + for (int n = 0; n < ncomp; n++) { soln_hat(i,j,k,n) += U_in(r,s,t,n) * alpha(i,j,k,1) * vfrac(r,s,t) / nrs(r,s,t); + } } } - for (int n = 0; n < ncomp; n++) + for (int n = 0; n < ncomp; n++) { soln_hat(i,j,k,n) /= nbhd_vol(i,j,k); + } } }); @@ -135,8 +139,9 @@ StateRedistribute ( Box const& bx, int ncomp, { if (bx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { - for (int n = 0; n < ncomp; n++) + for (int n = 0; n < ncomp; n++) { amrex::Gpu::Atomic::Add(&U_out(i,j,k,n),alpha(i,j,k,0)*nrs(i,j,k)*soln_hat(i,j,k,n)); + } } } else { @@ -170,11 +175,10 @@ StateRedistribute ( Box const& bx, int ncomp, #if (AMREX_SPACEDIM == 2) int kk = 0; #elif (AMREX_SPACEDIM == 3) - for(int kk(-1); kk<=1; kk++) + for(int kk(-1); kk<=1; kk++) { #endif - { - for(int jj(-1); jj<=1; jj++) - for(int ii(-1); ii<=1; ii++) + for(int jj(-1); jj<=1; jj++) { + for(int ii(-1); ii<=1; ii++) { if (flag(i,j,k).isConnected(ii,jj,kk)) { int r = i+ii; int s = j+jj; int t = k+kk; @@ -188,13 +192,13 @@ StateRedistribute ( Box const& bx, int ncomp, z_min = amrex::min(z_min, cent_hat(r,s,t,2)+static_cast(kk)); #endif } - } + AMREX_D_TERM(},},}) // If we need to grow the stencil, we let it be -nx:nx in the x-direction, // for example. Note that nx,ny,nz are either 1 or 2 - if ( (x_max-x_min) < slope_stencil_min_width ) nx = 2; - if ( (y_max-y_min) < slope_stencil_min_width ) ny = 2; + if ( (x_max-x_min) < slope_stencil_min_width ) { nx = 2; } + if ( (y_max-y_min) < slope_stencil_min_width ) { ny = 2; } #if (AMREX_SPACEDIM == 3) - if ( (z_max-z_min) < slope_stencil_min_width ) nz = 2; + if ( (z_max-z_min) < slope_stencil_min_width ) { nz = 2; } #endif amrex::GpuArray slopes_eb; diff --git a/Src/EB/AMReX_EB_chkpt_file.cpp b/Src/EB/AMReX_EB_chkpt_file.cpp index 3521a47d638..1acb6b4764d 100644 --- a/Src/EB/AMReX_EB_chkpt_file.cpp +++ b/Src/EB/AMReX_EB_chkpt_file.cpp @@ -40,8 +40,9 @@ ChkptFile::writeHeader (const BoxArray& cut_ba, const BoxArray& covered_ba, std::ofstream::trunc | std::ofstream::binary); - if ( ! HeaderFile.good() ) + if ( ! HeaderFile.good() ) { FileOpenFailed(HeaderFileName); + } HeaderFile.precision(17); @@ -51,17 +52,20 @@ ChkptFile::writeHeader (const BoxArray& cut_ba, const BoxArray& covered_ba, HeaderFile << nlevels << "\n"; // Geometry - for (int i = 0; i < AMREX_SPACEDIM; ++i) + for (int i = 0; i < AMREX_SPACEDIM; ++i) { HeaderFile << geom.ProbLo(i) << ' '; + } HeaderFile << '\n'; - for (int i = 0; i < AMREX_SPACEDIM; ++i) + for (int i = 0; i < AMREX_SPACEDIM; ++i) { HeaderFile << geom.ProbHi(i) << ' '; + } HeaderFile << '\n'; // ngrow - for (int i = 0; i < AMREX_SPACEDIM; ++i) + for (int i = 0; i < AMREX_SPACEDIM; ++i) { HeaderFile << ngrow[i] << ' '; + } HeaderFile << '\n'; // extend domain face @@ -117,7 +121,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, std::string File(m_restart_file + "/Header"); - if (amrex::Verbose()) amrex::Print() << "file=" << File << std::endl; + if (amrex::Verbose()) { amrex::Print() << "file=" << File << std::endl; } VisMF::IO_Buffer io_buffer(VisMF::GetIOBufferSize()); @@ -191,12 +195,12 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, AMREX_ALWAYS_ASSERT_WITH_MESSAGE(max_grid_size == mgs_chkptfile, "EB2::ChkptFile cannot read from different max_grid_size"); - if (amrex::Verbose()) amrex::Print() << "Loading cut_grids\n"; + if (amrex::Verbose()) { amrex::Print() << "Loading cut_grids\n"; } cut_grids.readFrom(is); gotoNextLine(is); if (is.peek() != EOF) { - if (amrex::Verbose()) amrex::Print() << "Loading covered_grids\n"; + if (amrex::Verbose()) { amrex::Print() << "Loading covered_grids\n"; } covered_grids.readFrom(is); gotoNextLine(is); } @@ -205,7 +209,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // volfrac { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_volfrac_name << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_volfrac_name << std::endl; } volfrac.define(cut_grids, dmap, 1, ng_gfab); @@ -215,7 +219,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // centroid { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_centroid_name << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_centroid_name << std::endl; } centroid.define(cut_grids, dmap, AMREX_SPACEDIM, ng_gfab); @@ -225,7 +229,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // bndryarea { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_bndryarea_name << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_bndryarea_name << std::endl; } bndryarea.define(cut_grids, dmap, 1, ng_gfab); @@ -235,7 +239,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // bndrycent { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_bndrycent_name << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_bndrycent_name << std::endl; } bndrycent.define(cut_grids, dmap, AMREX_SPACEDIM, ng_gfab); @@ -245,7 +249,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // bndrynorm { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_bndrynorm_name << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_bndrynorm_name << std::endl; } bndrynorm.define(cut_grids, dmap, AMREX_SPACEDIM, ng_gfab); @@ -256,7 +260,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { // areafrac { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_areafrac_name[idim] << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_areafrac_name[idim] << std::endl; } areafrac[idim].define(convert(cut_grids, IntVect::TheDimensionVector(idim)), dmap, 1, ng_gfab); @@ -266,7 +270,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // facecent { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_facecent_name[idim] << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_facecent_name[idim] << std::endl; } facecent[idim].define(convert(cut_grids, IntVect::TheDimensionVector(idim)), dmap, AMREX_SPACEDIM-1, ng_gfab); @@ -276,7 +280,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // edgecent { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_edgecent_name[idim] << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_edgecent_name[idim] << std::endl; } IntVect edge_type{1}; edge_type[idim] = 0; edgecent[idim].define(convert(cut_grids, edge_type), dmap, 1, ng_gfab); @@ -288,7 +292,7 @@ ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, // levelset { - if (amrex::Verbose()) amrex::Print() << " Loading " << m_levelset_name << std::endl; + if (amrex::Verbose()) { amrex::Print() << " Loading " << m_levelset_name << std::endl; } levelset.define(convert(cut_grids,IntVect::TheNodeVector()), dmap, 1, ng_gfab); diff --git a/Src/EB/AMReX_EB_utils.cpp b/Src/EB/AMReX_EB_utils.cpp index ff8c16f5d76..948a3d5db20 100644 --- a/Src/EB/AMReX_EB_utils.cpp +++ b/Src/EB/AMReX_EB_utils.cpp @@ -139,19 +139,19 @@ facets_nearest_pt (IntVect const& ind_pt, IntVect const& ind_loop, RealVect cons if ( std::abs(edge_v[0]) > eps ) { cx_lo = -( edge_p0[0] - static_cast( ind_loop[0] ) * dx[0] ) / edge_v[0]; cx_hi = -( edge_p0[0] - static_cast( ind_loop[0] + 1 ) * dx[0] ) / edge_v[0]; - if ( edge_v[0] < 0._rt ) amrex::Swap(cx_lo, cx_hi); + if ( edge_v[0] < 0._rt ) { amrex::Swap(cx_lo, cx_hi); } } // if ( std::abs(edge_v[1]) > eps ) { cy_lo = -( edge_p0[1] - static_cast( ind_loop[1] ) * dx[1] ) / edge_v[1]; cy_hi = -( edge_p0[1] - static_cast( ind_loop[1] + 1 ) * dx[1] ) / edge_v[1]; - if ( edge_v[1] < 0._rt ) amrex::Swap(cy_lo, cy_hi); + if ( edge_v[1] < 0._rt ) { amrex::Swap(cy_lo, cy_hi); } } // if ( std::abs(edge_v[2]) > eps ) { cz_lo = -( edge_p0[2] - static_cast( ind_loop[2] ) * dx[2] ) / edge_v[2]; cz_hi = -( edge_p0[2] - static_cast( ind_loop[2] + 1 ) * dx[2] ) / edge_v[2]; - if ( edge_v[2] < 0._rt ) amrex::Swap(cz_lo, cz_hi); + if ( edge_v[2] < 0._rt ) { amrex::Swap(cz_lo, cz_hi); } } // Real lambda_min = amrex::max(cx_lo, cy_lo, cz_lo); diff --git a/Src/EB/AMReX_WriteEBSurface.cpp b/Src/EB/AMReX_WriteEBSurface.cpp index 52b50d6c182..75a0421c94d 100644 --- a/Src/EB/AMReX_WriteEBSurface.cpp +++ b/Src/EB/AMReX_WriteEBSurface.cpp @@ -28,7 +28,7 @@ void WriteEBSurface (const BoxArray & ba, const DistributionMapping & dmap, cons const Box & bx = mfi.validbox(); if (my_flag.getType(bx) == FabType::covered || - my_flag.getType(bx) == FabType::regular) continue; + my_flag.getType(bx) == FabType::regular) { continue; } std::array areafrac; const MultiCutFab * bndrycent; @@ -62,7 +62,7 @@ void WriteEBSurface (const BoxArray & ba, const DistributionMapping & dmap, cons const Box & bx = mfi.validbox(); if (my_flag.getType(bx) == FabType::covered || - my_flag.getType(bx) == FabType::regular) continue; + my_flag.getType(bx) == FabType::regular) { continue; } eb_to_pvd.EBGridCoverage(cpu, problo, dx, bx, my_flag.const_array()); } diff --git a/Src/EB/AMReX_algoim.cpp b/Src/EB/AMReX_algoim.cpp index 08a4e2f5d54..254e15dab0f 100644 --- a/Src/EB/AMReX_algoim.cpp +++ b/Src/EB/AMReX_algoim.cpp @@ -30,7 +30,7 @@ compute_integrals (MultiFab& intgmf, IntVect nghost) const auto& flags = my_factory.getMultiEBCellFlagFab(); MFItInfo mfi_info; - if (Gpu::notInLaunchRegion()) mfi_info.EnableTiling().SetDynamic(true); + if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); } #ifdef AMREX_USE_OMP #pragma omp parallel if(Gpu::notInLaunchRegion()) @@ -73,7 +73,7 @@ compute_integrals (MultiFab& intgmf, IntVect nghost) if (ebflag.isRegular()) { set_regular(i,j,k,intg); } else if (ebflag.isCovered()) { - for (int n = 0; n < numIntgs; ++n) intg(i,j,k,n) = 0.0; + for (int n = 0; n < numIntgs; ++n) { intg(i,j,k,n) = 0.0; } } else { EBPlane phi(bc(i,j,k,0),bc(i,j,k,1),bc(i,j,k,2), bn(i,j,k,0),bn(i,j,k,1),bn(i,j,k,2)); @@ -125,15 +125,15 @@ compute_integrals (MultiFab& intgmf, IntVect nghost) { const auto lo = amrex::lbound(bx); const auto hi = amrex::ubound(bx); - for (int k = lo.z; k <= hi.z; ++k) - for (int j = lo.y; j <= hi.y; ++j) + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { const auto ebflag = fg(i,j,k); if (ebflag.isRegular()) { set_regular(i,j,k,intg); } else if (ebflag.isCovered()) { - for (int n = 0; n < numIntgs; ++n) intg(i,j,k,n) = 0.0; + for (int n = 0; n < numIntgs; ++n) { intg(i,j,k,n) = 0.0; } } else { EBPlane phi(bc(i,j,k,0),bc(i,j,k,1),bc(i,j,k,2), bn(i,j,k,0),bn(i,j,k,1),bn(i,j,k,2)); @@ -179,7 +179,7 @@ compute_integrals (MultiFab& intgmf, IntVect nghost) intg(i,j,k,i_S_xyz ) = q.eval([](Real x, Real y, Real z) noexcept { return x*y*z; }); } - } + }}} } } } @@ -213,7 +213,7 @@ compute_surface_integrals (MultiFab& sintgmf, IntVect nghost) const auto& barea = my_factory.getBndryArea(); MFItInfo mfi_info; - if (Gpu::notInLaunchRegion()) mfi_info.EnableTiling().SetDynamic(true); + if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); } #ifdef AMREX_USE_OMP #pragma omp parallel if(Gpu::notInLaunchRegion()) @@ -260,12 +260,12 @@ compute_surface_integrals (MultiFab& sintgmf, IntVect nghost) if (ebflag.isRegular()) { set_regular_surface(i,j,k,sintg); } else if (ebflag.isCovered()) { - for (int n = 0; n < numSurfIntgs; ++n) sintg(i,j,k,n) = 0.0; + for (int n = 0; n < numSurfIntgs; ++n) { sintg(i,j,k,n) = 0.0; } } else { constexpr Real almostone = Real(1.) - Real(100.)*std::numeric_limits::epsilon(); if (vf(i,j,k) >= almostone) { - for(int n = 0; n < numSurfIntgs; ++n) sintg(i,j,k,n) = 0.0; + for(int n = 0; n < numSurfIntgs; ++n) { sintg(i,j,k,n) = 0.0; } Real apxm = apx(i ,j ,k ); Real apxp = apx(i+1,j ,k ); @@ -317,20 +317,20 @@ compute_surface_integrals (MultiFab& sintgmf, IntVect nghost) { const auto lo = amrex::lbound(bx); const auto hi = amrex::ubound(bx); - for (int k = lo.z; k <= hi.z; ++k) - for (int j = lo.y; j <= hi.y; ++j) + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { const auto ebflag = fg(i,j,k); if (ebflag.isRegular()) { set_regular_surface(i,j,k,sintg); } else if (ebflag.isCovered()) { - for (int n = 0; n < numSurfIntgs; ++n) sintg(i,j,k,n) = 0.0; + for (int n = 0; n < numSurfIntgs; ++n) { sintg(i,j,k,n) = 0.0; } } else { constexpr Real almostone = Real(1.) - Real(100.)*std::numeric_limits::epsilon(); if (vf(i,j,k) >= almostone) { - for(int n = 0; n < numSurfIntgs; ++n) sintg(i,j,k,n) = 0.0; + for(int n = 0; n < numSurfIntgs; ++n) { sintg(i,j,k,n) = 0.0; } Real apxm = apx(i ,j ,k ); Real apxp = apx(i+1,j ,k ); @@ -376,7 +376,7 @@ compute_surface_integrals (MultiFab& sintgmf, IntVect nghost) { return x*y*z; }); } } - } + }}} } } } diff --git a/Src/EB/AMReX_algoim_K.H b/Src/EB/AMReX_algoim_K.H index c8a47614a08..09b7f1bd715 100644 --- a/Src/EB/AMReX_algoim_K.H +++ b/Src/EB/AMReX_algoim_K.H @@ -442,7 +442,7 @@ struct ImplicitIntegral // Loop over segments of divided interval for (int i = 0; i < nroots - 1; ++i) { - if (roots[i+1] - roots[i] < tol) continue; + if (roots[i+1] - roots[i] < tol) { continue; } // Evaluate sign of phi within segment and check for consistency with psi bool okay = true; @@ -457,7 +457,7 @@ struct ImplicitIntegral bool new_ok = (phi(x) > 0.0) ? (psi[j].sign() >= 0) : (psi[j].sign() <= 0); okay = okay && new_ok; } - if (!okay) continue; + if (!okay) { continue; } for (int j = 0; j < p; ++j) { @@ -521,8 +521,9 @@ struct ImplicitIntegral // integral domain is the entire hyperrectangle. if (psiCount == 0) { - if (!S) + if (!S) { tensorProductIntegral(); + } return; }