Skip to content

Commit

Permalink
modify the test
Browse files Browse the repository at this point in the history
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Jan 19, 2024
1 parent 1ff845b commit 69bd06d
Show file tree
Hide file tree
Showing 4 changed files with 81 additions and 35 deletions.
2 changes: 1 addition & 1 deletion Tests/LinearSolvers/CurlCurl/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
foreach(D IN LISTS AMReX_SPACEDIM)
if (NOT (D EQUAL 3))
if (D EQUAL 1)
return()
endif ()

Expand Down
32 changes: 24 additions & 8 deletions Tests/LinearSolvers/CurlCurl/MyTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,11 +24,11 @@ MyTest::solve ()
MLCurlCurl mlcc({geom}, {grids}, {dmap}, info);

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

mlcc.setScalars(alpha, beta);

Expand All @@ -40,6 +40,20 @@ MyTest::solve ()
mf.setVal(Real(0));
}
mlmg.solve({&solution}, {&rhs}, Real(1.0e-10), Real(0));

amrex::Print() << " Number of cells: " << n_cell << std::endl;
auto dvol = AMREX_D_TERM(geom.CellSize(0), *geom.CellSize(1), *geom.CellSize(2));
Array<std::string,3> names{"Ex", "Ey", "Ez"};
for (int idim = 0; idim < 3; ++idim) {
MultiFab::Subtract(solution[idim], exact[idim], 0, 0, 1, 0);
auto e0 = solution[idim].norminf();
auto e1 = solution[idim].norm1(0,geom.periodicity());
e1 *= dvol;
auto e2 = solution[idim].norm2(0,geom.periodicity());
e2 *= std::sqrt(dvol);
amrex::Print() << " " << names[idim] << " errors (max, L1, L2): "
<< e0 << " " << e1 << " " << e2 << std::endl;
}
}

void
Expand All @@ -64,7 +78,7 @@ void
MyTest::initData ()
{
RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(1.,1.,1.)});
Array<int,AMREX_SPACEDIM> is_periodic{AMREX_D_DECL(0,1,0)};
Array<int,AMREX_SPACEDIM> is_periodic{AMREX_D_DECL(0,0,1)};
Geometry::Setup(&rb, 0, is_periodic.data());
Box domain(IntVect(0), IntVect(n_cell-1));
geom.define(domain);
Expand All @@ -76,9 +90,11 @@ MyTest::initData ()
grids.maxSize(max_grid_size);
dmap.define(grids);

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
for (int idim = 0; idim < 3; ++idim) {
IntVect itype(1);
itype[idim] = 0;
if (idim < AMREX_SPACEDIM) {
itype[idim] = 0;
}
BoxArray const& ba = amrex::convert(grids, itype);
solution[idim].define(ba,dmap,1,1);
exact [idim].define(ba,dmap,1,1);
Expand All @@ -87,7 +103,7 @@ MyTest::initData ()

initProb();

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
for (int idim = 0; idim < 3; ++idim) {
exact[idim].LocalCopy(solution[idim], 0, 0, 1, IntVect(1));
}
}
14 changes: 6 additions & 8 deletions Tests/LinearSolvers/CurlCurl/initProb.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,12 @@ MyTest::initProb ()
for (MFIter mfi(rhs[0], TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& gbx = mfi.tilebox(IntVect(1),IntVect(1));
GpuArray<Array4<Real>,AMREX_SPACEDIM> rhsfab
{AMREX_D_DECL(rhs[0].array(mfi),
rhs[1].array(mfi),
rhs[2].array(mfi))};
GpuArray<Array4<Real>,AMREX_SPACEDIM> solfab
{AMREX_D_DECL(solution[0].array(mfi),
solution[1].array(mfi),
solution[2].array(mfi))};
GpuArray<Array4<Real>,3> rhsfab{rhs[0].array(mfi),
rhs[1].array(mfi),
rhs[2].array(mfi)};
GpuArray<Array4<Real>,3> solfab{solution[0].array(mfi),
solution[1].array(mfi),
solution[2].array(mfi)};
amrex::ParallelFor(gbx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Expand Down
68 changes: 50 additions & 18 deletions Tests/LinearSolvers/CurlCurl/initProb_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,70 +5,102 @@

AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void actual_init_prob (int i, int j, int k,
amrex::GpuArray<amrex::Array4<amrex::Real>,AMREX_SPACEDIM> const& rhs,
amrex::GpuArray<amrex::Array4<amrex::Real>,AMREX_SPACEDIM> const& sol,
amrex::GpuArray<amrex::Array4<amrex::Real>,3> const& rhs,
amrex::GpuArray<amrex::Array4<amrex::Real>,3> const& sol,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dx,
amrex::Real alpha, amrex::Real beta)
{
using namespace amrex;

constexpr amrex::Real tpi = Real(2.)*amrex::Math::pi<Real>();
constexpr amrex::Real fpi = Real(4.)*amrex::Math::pi<Real>();
constexpr Real pi = amrex::Math::pi<Real>();

Real xnd = problo[0] + Real(i)*dx[0];
Real ynd = problo[1] + Real(j)*dx[1];
Real znd = problo[2] + Real(k)*dx[2];
Real xcc = xnd + Real(0.5)*dx[0];
Real ycc = ynd + Real(0.5)*dx[1];
#if (AMREX_SPACEDIM == 3)
Real znd = problo[2] + Real(k)*dx[2];
Real zcc = znd + Real(0.5)*dx[2];
#endif

if (sol[0].contains(i,j,k)) {
Real x = xcc;
Real y = ynd;
Real Ex = std::sin(pi*x) * std::sin(Real(2.5)*pi*y);
#if (AMREX_SPACEDIM == 3)
Real z = znd;
sol[0](i,j,k) = std::cos(tpi*x) * std::sin(fpi*y) * std::sin(fpi*z);
}
Ex *= std::sin(Real(2.0)*pi*z + Real(1./3.)*pi);
#endif
sol[0](i,j,k) = Ex;
}

if (sol[1].contains(i,j,k)) {
Real x = xnd;
Real y = ycc;
Real Ey = std::cos(Real(2.5)*pi*x) * std::sin(Real(3.)*pi*y);
#if (AMREX_SPACEDIM == 3)
Real z = znd;
sol[1](i,j,k) = std::sin(fpi*x) * std::cos(tpi*y) * std::sin(fpi*z);
Ey *= std::sin(Real(4.)*pi*z + Real(0.25)*pi);
#endif
sol[1](i,j,k) = Ey;
}

if (sol[2].contains(i,j,k)) {
Real x = xnd;
Real y = ynd;
Real Ez = std::cos(Real(3.5)*pi*x) * std::sin(Real(3.5)*pi*y);
#if (AMREX_SPACEDIM == 3)
Real z = zcc;
sol[2](i,j,k) = std::sin(fpi*x) * std::sin(fpi*y) * std::cos(tpi*z);
Ez *= std::sin(Real(4.)*pi*z + Real(1./6.)*pi);
#endif
sol[2](i,j,k) = Ez;
}

if (rhs[0].contains(i,j,k)) {
Real x = xcc;
Real y = ynd;
#if (AMREX_SPACEDIM == 2)
Real cce = Real(-7.5)*pi*pi*std::sin(Real(2.5)*pi*x)*std::cos(Real(3.)*pi*y)
+ Real(6.25)*pi*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y);
#else
Real z = znd;
rhs[0](i,j,k) = (beta + alpha*Real(4.0)*tpi*fpi) * sol[0](i,j,k)
- alpha*tpi*fpi * std::cos(fpi*x) * (std::sin(tpi*y) * std::sin(fpi*z) +
std::sin(fpi*y) * std::sin(tpi*z));
Real cce = Real(-7.5)*pi*pi*std::sin(Real(2.5)*pi*x)*std::cos(Real(3.)*pi*y)*std::sin(Real(4.)*pi*z+Real(0.25)*pi)
+ Real(6.25)*pi*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y)*std::sin(Real(2.)*pi*z+Real(1./3.)*pi)
- Real(14.)*pi*pi*std::sin(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y)*std::cos(Real(4.)*pi*z+Real(1./6.)*pi)
+ Real(4.)*pi*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y)*std::sin(Real(2.)*pi*z+Real(1./3.)*pi);
#endif
rhs[0](i,j,k) = alpha*cce + beta*sol[0](i,j,k);
}

if (rhs[1].contains(i,j,k)) {
Real x = xnd;
Real y = ycc;
#if (AMREX_SPACEDIM == 2)
Real cce = Real(6.25)*pi*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y)
+ Real(2.5)*pi*pi*std::cos(pi*x)*std::cos(Real(2.5)*pi*y);
#else
Real z = znd;
rhs[1](i,j,k) = (beta + alpha*Real(4.0)*tpi*fpi) * sol[1](i,j,k)
- alpha*tpi*fpi * std::cos(fpi*y) * (std::sin(tpi*x) * std::sin(fpi*z) +
std::sin(fpi*x) * std::sin(tpi*z));
Real cce = Real(6.25)*pi*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y)*std::sin(Real(4.)*pi*z+Real(0.25)*pi)
+ Real(2.5)*pi*pi*std::cos(pi*x)*std::cos(Real(2.5)*pi*y)*std::sin(Real(2.)*pi*z+Real(1./3.)*pi)
+ Real(14.)*pi*pi*std::cos(Real(3.5)*pi*x)*std::cos(Real(3.5)*pi*y)*std::cos(Real(4.)*pi*z+Real(1./6.)*pi)
+ Real(16.)*pi*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y)*std::sin(Real(4.)*pi*z+Real(0.25)*pi);
#endif
rhs[1](i,j,k) = alpha*cce + beta*sol[1](i,j,k);
}

if (rhs[2].contains(i,j,k)) {
Real x = xnd;
Real y = ynd;
#if (AMREX_SPACEDIM == 2)
Real cce = Real(24.5)*pi*pi*std::cos(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y);
#else
Real z = zcc;
rhs[2](i,j,k) = (beta + alpha*Real(4.0)*tpi*fpi) * sol[2](i,j,k)
- alpha*tpi*fpi * std::cos(fpi*z) * (std::sin(tpi*x) * std::sin(fpi*y) +
std::sin(fpi*x) * std::sin(tpi*y));
Real cce = Real(24.5)*pi*pi*std::cos(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y)*std::sin(Real(4.)*pi*z+Real(1./6.)*pi)
+ Real(2.)*pi*pi*std::cos(pi*x)*std::sin(Real(2.5)*pi*y)*std::cos(Real(2.)*pi*z+Real(1./3.)*pi)
+ Real(12.)*pi*pi*std::cos(Real(2.5)*pi*x)*std::cos(Real(3.)*pi*y)*std::cos(Real(4.)*pi*z+Real(0.25)*pi);
#endif
rhs[2](i,j,k) = alpha*cce + beta*sol[2](i,j,k);
}
}

Expand Down

0 comments on commit 69bd06d

Please sign in to comment.