Skip to content

Commit

Permalink
Modify Tests/LinearSolvers/ABecLaplacian_C (#3888)
Browse files Browse the repository at this point in the history 
Add Dirichlet boundaries to the mix for testing.
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed May 3, 2024
1 parent ee11254 commit c361681
Show file tree
Hide file tree
Showing 3 changed files with 59 additions and 86 deletions.
65 changes: 12 additions & 53 deletions Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,18 +154,16 @@ MyTest::solveABecLaplacian ()

mlabec.setMaxOrder(linop_maxorder);

// This is a 3d problem with homogeneous Neumann BC
mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann,
mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet,
LinOpBCType::Neumann,
LinOpBCType::Neumann)},
{AMREX_D_DECL(LinOpBCType::Neumann,
LinOpBCType::Neumann,
LinOpBCType::Dirichlet,
LinOpBCType::Neumann)});

for (int ilev = 0; ilev < nlevels; ++ilev)
{
// for problem with pure homogeneous Neumann BC, we could pass a nullptr
mlabec.setLevelBC(ilev, nullptr);
mlabec.setLevelBC(ilev, &solution[ilev]);
}

mlabec.setScalars(ascalar, bscalar);
Expand Down Expand Up @@ -213,20 +211,18 @@ MyTest::solveABecLaplacian ()

mlabec.setMaxOrder(linop_maxorder);

// This is a 3d problem with homogeneous Neumann BC
mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann,
mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet,
LinOpBCType::Neumann,
LinOpBCType::Neumann)},
{AMREX_D_DECL(LinOpBCType::Neumann,
LinOpBCType::Neumann,
LinOpBCType::Dirichlet,
LinOpBCType::Neumann)});

if (ilev > 0) {
mlabec.setCoarseFineBC(&solution[ilev-1], ref_ratio);
}

// for problem with pure homogeneous Neumann BC, we could pass a nullptr
mlabec.setLevelBC(0, nullptr);
mlabec.setLevelBC(0, &solution[ilev]);

mlabec.setScalars(ascalar, bscalar);

Expand Down Expand Up @@ -263,18 +259,6 @@ MyTest::solveABecLaplacian ()
mlmg.solve({&solution[ilev]}, {&rhs[ilev]}, tol_rel, tol_abs);
}
}

// Since this problem has Neumann BC, solution + constant is also a
// solution. So we are going to shift the solution by a constant
// for comparison with the "exact solution".
// The statement above is incorrect because we have the a term, albeit small.
const Real npts = grids[0].d_numPts();
const Real avg1 = exact_solution[0].sum();
const Real avg2 = solution[0].sum();
const Real offset = (avg1-avg2)/npts;
for (int ilev = 0; ilev < nlevels; ++ilev) {
solution[ilev].plus(offset, 0, 1, 0);
}
}

void
Expand Down Expand Up @@ -406,18 +390,6 @@ MyTest::solveABecLaplacianInhomNeumann ()
mlmg.solve({&solution[ilev]}, {&rhs[ilev]}, tol_rel, tol_abs);
}
}

// Since this problem has Neumann BC, solution + constant is also a
// solution. So we are going to shift the solution by a constant
// for comparison with the "exact solution".
// The statement above is incorrect because we have the a term, albeit small.
const Real npts = grids[0].d_numPts();
const Real avg1 = exact_solution[0].sum();
const Real avg2 = solution[0].sum();
const Real offset = (avg1-avg2)/npts;
for (int ilev = 0; ilev < nlevels; ++ilev) {
solution[ilev].plus(offset, 0, 1, 0);
}
}

void
Expand Down Expand Up @@ -481,31 +453,30 @@ MyTest::solveABecLaplacianGMRES ()
const auto tol_rel = Real(1.e-10);
const auto tol_abs = Real(0.0);

AMREX_ALWAYS_ASSERT_WITH_MESSAGE(composite_solve == false,
"solveABecLaplacianGMRES does not support composite solve");

const auto nlevels = static_cast<int>(geom.size());

AMREX_ALWAYS_ASSERT_WITH_MESSAGE(composite_solve == false || nlevels == 1,
"solveABecLaplacianGMRES does not support composite solve");

for (int ilev = 0; ilev < nlevels; ++ilev)
{
MLABecLaplacian mlabec({geom[ilev]}, {grids[ilev]}, {dmap[ilev]}, info);

mlabec.setMaxOrder(linop_maxorder);

// This is a 3d problem with homogeneous Neumann BC
mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann,
mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet,
LinOpBCType::Neumann,
LinOpBCType::Neumann)},
{AMREX_D_DECL(LinOpBCType::Neumann,
LinOpBCType::Neumann,
LinOpBCType::Dirichlet,
LinOpBCType::Neumann)});

if (ilev > 0) {
mlabec.setCoarseFineBC(&solution[ilev-1], ref_ratio);
}

// for problem with pure homogeneous Neumann BC, we could pass a nullptr
mlabec.setLevelBC(0, nullptr);
mlabec.setLevelBC(0, &solution[ilev]);

mlabec.setScalars(ascalar, bscalar);

Expand Down Expand Up @@ -536,18 +507,6 @@ MyTest::solveABecLaplacianGMRES ()
<< " " << res.norm1(0) << " " << res.norm2(0) << '\n';
}
}

// Since this problem has Neumann BC, solution + constant is also a
// solution. So we are going to shift the solution by a constant
// for comparison with the "exact solution".
// The statement above is incorrect because we have the a term, albeit small.
const Real npts = grids[0].d_numPts();
const Real avg1 = exact_solution[0].sum();
const Real avg2 = solution[0].sum();
const Real offset = (avg1-avg2)/npts;
for (int ilev = 0; ilev < nlevels; ++ilev) {
solution[ilev].plus(offset, 0, 1, 0);
}
}

void
Expand Down
15 changes: 4 additions & 11 deletions Tests/LinearSolvers/ABecLaplacian_C/initProb.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,29 +45,22 @@ MyTest::initProbABecLaplacian ()
#endif
for (MFIter mfi(rhs[ilev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Box& vbx = mfi.validbox();
const Box& gbx = mfi.growntilebox(1);

auto rhsfab = rhs[ilev].array(mfi);
auto solfab = solution[ilev].array(mfi);
auto exactfab = exact_solution[ilev].array(mfi);
auto acoeffab = acoef[ilev].array(mfi);
auto bcoeffab = bcoef[ilev].array(mfi);

amrex::ParallelFor(gbx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
actual_init_bcoef(i,j,k,bcoeffab,prob_lo,prob_hi,dx);
});

amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
actual_init_abeclap(i,j,k,rhsfab,exactfab,acoeffab,bcoeffab,
a,b,prob_lo,prob_hi,dx);
actual_init_abeclap(i,j,k,rhsfab,solfab,exactfab,acoeffab,bcoeffab,
a,b,prob_lo,prob_hi,dx,vbx);
});
}

solution[ilev].setVal(0.0);
}
}

Expand Down
65 changes: 43 additions & 22 deletions Tests/LinearSolvers/ABecLaplacian_C/initProb_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#define INIT_PROB_K_H_

#include <AMReX_FArrayBox.H>
#include <algorithm>

AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void actual_init_poisson (int i, int j, int k,
Expand Down Expand Up @@ -70,14 +71,16 @@ void actual_init_bcoef (int i, int j, int k,

AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void actual_init_abeclap (int i, int j, int k,
amrex::Array4<amrex::Real > const& rhs,
amrex::Array4<amrex::Real > const& exact,
amrex::Array4<amrex::Real > const& alpha,
amrex::Array4<amrex::Real const> const& beta,
amrex::Array4<amrex::Real> const& rhs,
amrex::Array4<amrex::Real> const& sol,
amrex::Array4<amrex::Real> const& exact,
amrex::Array4<amrex::Real> const& alpha,
amrex::Array4<amrex::Real> const& beta,
amrex::Real a, amrex::Real b,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& prob_lo,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& prob_hi,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dx)
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dx,
amrex::Box const& vbx)
{
constexpr amrex::Real w = 0.05;
constexpr amrex::Real sigma = 10.;
Expand Down Expand Up @@ -105,25 +108,43 @@ void actual_init_abeclap (int i, int j, int k,
#endif

amrex::Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc));
amrex::Real tmp = std::cosh(theta*(r-0.25));
amrex::Real dbdrfac = (sigma-1.)/2./(tmp*tmp) * theta/r;
dbdrfac *= b;

alpha(i,j,k) = 1.;

exact(i,j,k) = std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ .25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z);
beta(i,j,k) = (sigma-1.)/2.*std::tanh(theta*(r-0.25)) + (sigma+1.)/2.;

rhs(i,j,k) = beta(i,j,k)*b*fac*(std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z))
+ dbdrfac*((x-xc)*(tpi*std::sin(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ pi*std::sin(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z))
+ (y-yc)*(tpi*std::cos(tpi*x) * std::sin(tpi*y) * std::cos(tpi*z)
+ pi*std::cos(fpi*x) * std::sin(fpi*y) * std::cos(fpi*z))
+ (z-zc)*(tpi*std::cos(tpi*x) * std::cos(tpi*y) * std::sin(tpi*z)
+ pi*std::cos(fpi*x) * std::cos(fpi*y) * std::sin(fpi*z)))
+ a * (std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ 0.25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z));
if (vbx.contains(i,j,k)) {
amrex::Real tmp = std::cosh(theta*(r-0.25));
amrex::Real dbdrfac = (sigma-1.)/2./(tmp*tmp) * theta/r;
dbdrfac *= b;

alpha(i,j,k) = 1.;

sol(i,j,k) = 0.; // Use zero as initial guess for the linear solver

exact(i,j,k) = std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ .25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z);

rhs(i,j,k) = beta(i,j,k)*b*fac*(std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z))
+ dbdrfac*((x-xc)*(tpi*std::sin(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ pi*std::sin(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z))
+ (y-yc)*(tpi*std::cos(tpi*x) * std::sin(tpi*y) * std::cos(tpi*z)
+ pi*std::cos(fpi*x) * std::sin(fpi*y) * std::cos(fpi*z))
+ (z-zc)*(tpi*std::cos(tpi*x) * std::cos(tpi*y) * std::sin(tpi*z)
+ pi*std::cos(fpi*x) * std::cos(fpi*y) * std::sin(fpi*z)))
+ a * (std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)
+ 0.25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z));
} else {
amrex::Real xb = std::clamp(x, prob_lo[0], prob_hi[0]);
amrex::Real yb = std::clamp(y, prob_lo[1], prob_hi[1]);
#if (AMREX_SPACEDIM == 2)
amrex::Real zb = 0.0;
#else
amrex::Real zb = std::clamp(z, prob_lo[2], prob_hi[2]);
#endif
// provide Dirichlet BC
sol(i,j,k) = std::cos(tpi*xb) * std::cos(tpi*yb) * std::cos(tpi*zb)
+ .25 * std::cos(fpi*xb) * std::cos(fpi*yb) * std::cos(fpi*zb);
}
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE
Expand Down

0 comments on commit c361681

Please sign in to comment.