Skip to content

Commit

Permalink
Symmetry boundary
Browse files Browse the repository at this point in the history
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Jan 17, 2024
1 parent 149ab53 commit 71cb2b1
Show file tree
Hide file tree
Showing 7 changed files with 122 additions and 142 deletions.
3 changes: 2 additions & 1 deletion Src/Boundary/AMReX_LO_BCTYPES.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#define AMREX_LO_INFLOW 106
#define AMREX_LO_INHOMOG_NEUMANN 107
#define AMREX_LO_ROBIN 108
#define AMREX_LO_SYMMETRY 109
#define AMREX_LO_PERIODIC 200
#define AMREX_LO_BOGUS 1729

Expand All @@ -38,6 +39,7 @@ namespace amrex {
inflow = AMREX_LO_INFLOW,
inhomogNeumann = AMREX_LO_INHOMOG_NEUMANN,
Robin = AMREX_LO_ROBIN,
symmetry = AMREX_LO_SYMMETRY,
Periodic = AMREX_LO_PERIODIC,
bogus = AMREX_LO_BOGUS
};
Expand All @@ -48,4 +50,3 @@ namespace amrex {
#endif

#endif

5 changes: 5 additions & 0 deletions Src/Boundary/AMReX_LO_BCTYPES.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,11 @@ std::ostream& operator<< (std::ostream& os, const LinOpBCType& t)
os << "Robin";
break;
}
case LinOpBCType::symmetry:
{
os << "symmetry";
break;
}
case LinOpBCType::Periodic:
{
os << "Periodic";
Expand Down
16 changes: 4 additions & 12 deletions Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,25 +116,17 @@ public:

// public for cuda

void smooth_x (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const;
void smooth_y (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const;
void smooth_z (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const;
void smooth (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const;

void compresid (int amrlev, int mglev, MF& resid, MF const& b) const;

void applyPhysBC_x (int amrlev, int mglev, MultiFab& mf) const;
void applyPhysBC_y (int amrlev, int mglev, MultiFab& mf) const;
void applyPhysBC_z (int amrlev, int mglev, MultiFab& mf) const;
void applyPhysBC (int amrlev, int mglev, MultiFab& mf) const;

private:

void applyBC (int amrlev, int mglev, MF& in) const;
void applyBC_x (int amrlev, int mglev, MultiFab& mf) const;
void applyBC_y (int amrlev, int mglev, MultiFab& mf) const;
void applyBC_z (int amrlev, int mglev, MultiFab& mf) const;
void applyBC (int amrlev, int mglev, MultiFab& mf) const;

[[nodiscard]] iMultiFab const& getDotMask (int amrlev, int mglev, int idim) const;

Expand Down
211 changes: 85 additions & 126 deletions Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,102 +153,32 @@ void MLCurlCurl::smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
applyBC(amrlev, mglev, sol);
}

smooth_x(amrlev, mglev, sol, rhs[0], 0);
applyBC_x(amrlev, mglev, sol[0]);
smooth(amrlev, mglev, sol, rhs[0], 0); // Ex red
applyBC(amrlev, mglev, sol[0]);

smooth_y(amrlev, mglev, sol, rhs[1], 0);
applyBC_y(amrlev, mglev, sol[1]);
smooth(amrlev, mglev, sol, rhs[1], 0); // Ey red
applyBC(amrlev, mglev, sol[1]);

smooth_z(amrlev, mglev, sol, rhs[2], 0);
applyBC_z(amrlev, mglev, sol[2]);
smooth(amrlev, mglev, sol, rhs[2], 0); // Ez red
applyBC(amrlev, mglev, sol[2]);

smooth_x(amrlev, mglev, sol, rhs[0], 1);
applyBC_x(amrlev, mglev, sol[0]);
smooth(amrlev, mglev, sol, rhs[0], 1); // Ex black
applyBC(amrlev, mglev, sol[0]);

smooth_y(amrlev, mglev, sol, rhs[1], 1);
applyBC_y(amrlev, mglev, sol[1]);
smooth(amrlev, mglev, sol, rhs[1], 1); // Ey black
applyBC(amrlev, mglev, sol[1]);

smooth_z(amrlev, mglev, sol, rhs[2], 1);
smooth(amrlev, mglev, sol, rhs[2], 1); // Ez black

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
amrex::OverrideSync(sol[idim], getDotMask(amrlev,mglev,idim),
this->m_geom[amrlev][mglev].periodicity());
}
}

void MLCurlCurl::smooth_x (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const
void MLCurlCurl::smooth (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const
{
amrex::ignore_unused(amrlev,mglev);

auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray();
auto const a = m_alpha;
auto const b = m_beta;

int const dirichlet_ylo = getDirichlet(amrlev, mglev, 1, 0);
int const dirichlet_yhi = getDirichlet(amrlev, mglev, 1, 1);
int const dirichlet_zlo = getDirichlet(amrlev, mglev, 2, 0);
int const dirichlet_zhi = getDirichlet(amrlev, mglev, 2, 1);

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(sol[0],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& bx = mfi.tilebox();
auto const& ex = sol[0].array(mfi);
auto const& ey = sol[1].const_array(mfi);
auto const& ez = sol[2].const_array(mfi);
auto const& rh = rhs.const_array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (j != dirichlet_ylo && j != dirichlet_yhi &&
k != dirichlet_zlo && k != dirichlet_zhi) {
mlcurlcurl_gsrb_x(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack);
}
});
}
}

void MLCurlCurl::smooth_y (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const
{
amrex::ignore_unused(amrlev,mglev);

auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray();
auto const a = m_alpha;
auto const b = m_beta;

int const dirichlet_xlo = getDirichlet(amrlev, mglev, 0, 0);
int const dirichlet_xhi = getDirichlet(amrlev, mglev, 0, 1);
int const dirichlet_zlo = getDirichlet(amrlev, mglev, 2, 0);
int const dirichlet_zhi = getDirichlet(amrlev, mglev, 2, 1);

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(sol[1],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& bx = mfi.tilebox();
auto const& ex = sol[0].const_array(mfi);
auto const& ey = sol[1]. array(mfi);
auto const& ez = sol[2].const_array(mfi);
auto const& rh = rhs.const_array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (i != dirichlet_xlo && i != dirichlet_xhi &&
k != dirichlet_zlo && k != dirichlet_zhi) {
mlcurlcurl_gsrb_y(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack);
}
});
}
}

void MLCurlCurl::smooth_z (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int redblack) const
{
amrex::ignore_unused(amrlev,mglev);

auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray();
auto const a = m_alpha;
auto const b = m_beta;
Expand All @@ -257,24 +187,50 @@ void MLCurlCurl::smooth_z (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
int const dirichlet_xhi = getDirichlet(amrlev, mglev, 0, 1);
int const dirichlet_ylo = getDirichlet(amrlev, mglev, 1, 0);
int const dirichlet_yhi = getDirichlet(amrlev, mglev, 1, 1);
int const dirichlet_zlo = getDirichlet(amrlev, mglev, 2, 0);
int const dirichlet_zhi = getDirichlet(amrlev, mglev, 2, 1);

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(sol[2],TilingIfNotGPU()); mfi.isValid(); ++mfi)
for (MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& bx = mfi.tilebox();
auto const& ex = sol[0].const_array(mfi);
auto const& ey = sol[1].const_array(mfi);
auto const& ez = sol[2]. array(mfi);
auto const& rh = rhs.const_array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (i != dirichlet_xlo && i != dirichlet_xhi &&
j != dirichlet_ylo && j != dirichlet_yhi) {
mlcurlcurl_gsrb_z(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack);
}
});
if (rhs.ixType() == sol[0].ixType()) {
auto const& ex = sol[0].array(mfi);
auto const& ey = sol[1].const_array(mfi);
auto const& ez = sol[2].const_array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (j != dirichlet_ylo && j != dirichlet_yhi &&
k != dirichlet_zlo && k != dirichlet_zhi) {
mlcurlcurl_gsrb_x(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack);
}
});
} else if (rhs.ixType() == sol[1].ixType()) {
auto const& ex = sol[0].const_array(mfi);
auto const& ey = sol[1].array(mfi);
auto const& ez = sol[2].const_array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (i != dirichlet_xlo && i != dirichlet_xhi &&
k != dirichlet_zlo && k != dirichlet_zhi) {
mlcurlcurl_gsrb_y(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack);
}
});
} else {
auto const& ex = sol[0].const_array(mfi);
auto const& ey = sol[1].const_array(mfi);
auto const& ez = sol[2].array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (i != dirichlet_xlo && i != dirichlet_xhi &&
j != dirichlet_ylo && j != dirichlet_yhi) {
mlcurlcurl_gsrb_z(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack);
}
});
}
}
}

Expand Down Expand Up @@ -483,42 +439,45 @@ void MLCurlCurl::applyBC (int amrlev, int mglev, MF& in) const
{
Vector<MultiFab*> mfs{AMREX_D_DECL(in.data(),&(in[1]),&(in[2]))};
FillBoundary(mfs, this->m_geom[amrlev][mglev].periodicity());
applyPhysBC_x(amrlev, mglev, in[0]);
applyPhysBC_y(amrlev, mglev, in[1]);
applyPhysBC_z(amrlev, mglev, in[2]);
applyPhysBC(amrlev, mglev, in[0]);
applyPhysBC(amrlev, mglev, in[1]);
applyPhysBC(amrlev, mglev, in[2]);
}

void MLCurlCurl::applyBC_x (int amrlev, int mglev, MultiFab& mf) const
void MLCurlCurl::applyBC (int amrlev, int mglev, MultiFab& mf) const
{
mf.FillBoundary(this->m_geom[amrlev][mglev].periodicity());
applyPhysBC_x(amrlev, mglev, mf);
}

void MLCurlCurl::applyBC_y (int amrlev, int mglev, MultiFab& mf) const
{
mf.FillBoundary(this->m_geom[amrlev][mglev].periodicity());
applyPhysBC_y(amrlev, mglev, mf);
}

void MLCurlCurl::applyBC_z (int amrlev, int mglev, MultiFab& mf) const
{
mf.FillBoundary(this->m_geom[amrlev][mglev].periodicity());
applyPhysBC_z(amrlev, mglev, mf);
}

void MLCurlCurl::applyPhysBC_x (int amrlev, int mglev, MultiFab& mf) const
{

}

void MLCurlCurl::applyPhysBC_y (int amrlev, int mglev, MultiFab& mf) const
{

}

void MLCurlCurl::applyPhysBC_z (int amrlev, int mglev, MultiFab& mf) const
{

applyPhysBC(amrlev, mglev, mf);
}

void MLCurlCurl::applyPhysBC (int amrlev, int mglev, MultiFab& mf) const
{
auto const idxtype = mf.ixType();
Box const domain = amrex::convert(this->m_geom[amrlev][mglev].Domain(), idxtype);
for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
auto const& vbx = mfi.validbox();
auto const& a = mf.array(mfi);
for (OrientationIter oit; oit; ++oit) {
Orientation const face = oit();
int const idim = face.coordDir();
bool is_symmetric = face.isLow()
? m_lobc[0][idim] == LinOpBCType::symmetry
: m_hibc[0][idim] == LinOpBCType::symmetry;
if (domain[face] == vbx[face] && is_symmetric) {
Box b = vbx;
int shift = face.isLow() ? -1 : 1;
b.setRange(idim, domain[face] + shift, 1);
#ifdef AMREX_USE_GPU
static_assert(false, "MLCurlCurl: todo");
#else
amrex::LoopOnCpu(b, [&] (int i, int j, int k)
{
mlcurlcurl_bc_symmetry(i, j, k, face, idxtype, a);
});
#endif
}
}
}
}

iMultiFab const& MLCurlCurl::getDotMask (int amrlev, int mglev, int idim) const
Expand Down
23 changes: 23 additions & 0 deletions Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,6 +257,29 @@ void mlcurlcurl_interpadd (int dir, int i, int j, int k,
}
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void mlcurlcurl_bc_symmetry (int i, int j, int k, Orientation face, IndexType it,
Array4<Real> const& a)
{
int const idir = face.coordDir();
int offset = face.isLow() ? 1 : -1;
Real sign;
if (it.cellCentered()) {
sign = Real(-1.0);
} else {
sign = Real(1.0);
offset *= 2;
}

if (idir == 0) {
a(i,j,k) = sign * a(i+offset,j,k);
} else if (idir == 1) {
a(i,j,k) = sign * a(i,j+offset,k);
} else {
a(i,j,k) = sign * a(i,j,k+offset);
}
}

}

#endif
2 changes: 1 addition & 1 deletion Tests/LinearSolvers/CurlCurl/MyTest.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ private:
amrex::Array<amrex::MultiFab,AMREX_SPACEDIM> exact;
amrex::Array<amrex::MultiFab,AMREX_SPACEDIM> rhs;

amrex::Real alpha_over_dx2 = 1.e-1;
amrex::Real alpha_over_dx2 = 1.0;
amrex::Real alpha;
amrex::Real beta = 1.0;
};
Expand Down
4 changes: 2 additions & 2 deletions Tests/LinearSolvers/CurlCurl/MyTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,12 @@ MyTest::solve ()

MLCurlCurl mlcc({geom}, {grids}, {dmap}, info);

mlcc.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet,
mlcc.setDomainBC({AMREX_D_DECL(LinOpBCType::symmetry,
LinOpBCType::Periodic,
LinOpBCType::Dirichlet)},
{AMREX_D_DECL(LinOpBCType::Dirichlet,
LinOpBCType::Periodic,
LinOpBCType::Dirichlet)});
LinOpBCType::symmetry)});

mlcc.setScalars(alpha, beta);

Expand Down

0 comments on commit 71cb2b1

Please sign in to comment.