From 69bd06dddcf1c7df4ed3befc575ccfb3df571d9d Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Thu, 18 Jan 2024 12:08:32 -0800 Subject: [PATCH] modify the test --- Tests/LinearSolvers/CurlCurl/CMakeLists.txt | 2 +- Tests/LinearSolvers/CurlCurl/MyTest.cpp | 32 +++++++--- Tests/LinearSolvers/CurlCurl/initProb.cpp | 14 ++--- Tests/LinearSolvers/CurlCurl/initProb_K.H | 68 +++++++++++++++------ 4 files changed, 81 insertions(+), 35 deletions(-) diff --git a/Tests/LinearSolvers/CurlCurl/CMakeLists.txt b/Tests/LinearSolvers/CurlCurl/CMakeLists.txt index afc06542acd..9dacdeb2fea 100644 --- a/Tests/LinearSolvers/CurlCurl/CMakeLists.txt +++ b/Tests/LinearSolvers/CurlCurl/CMakeLists.txt @@ -1,5 +1,5 @@ foreach(D IN LISTS AMReX_SPACEDIM) - if (NOT (D EQUAL 3)) + if (D EQUAL 1) return() endif () diff --git a/Tests/LinearSolvers/CurlCurl/MyTest.cpp b/Tests/LinearSolvers/CurlCurl/MyTest.cpp index 882a6834595..0aa29d79e1b 100644 --- a/Tests/LinearSolvers/CurlCurl/MyTest.cpp +++ b/Tests/LinearSolvers/CurlCurl/MyTest.cpp @@ -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); @@ -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 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 @@ -64,7 +78,7 @@ void MyTest::initData () { RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(1.,1.,1.)}); - Array is_periodic{AMREX_D_DECL(0,1,0)}; + Array 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); @@ -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); @@ -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)); } } diff --git a/Tests/LinearSolvers/CurlCurl/initProb.cpp b/Tests/LinearSolvers/CurlCurl/initProb.cpp index 0d8d98d36cd..3acf779bd05 100644 --- a/Tests/LinearSolvers/CurlCurl/initProb.cpp +++ b/Tests/LinearSolvers/CurlCurl/initProb.cpp @@ -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,AMREX_SPACEDIM> rhsfab - {AMREX_D_DECL(rhs[0].array(mfi), - rhs[1].array(mfi), - rhs[2].array(mfi))}; - GpuArray,AMREX_SPACEDIM> solfab - {AMREX_D_DECL(solution[0].array(mfi), - solution[1].array(mfi), - solution[2].array(mfi))}; + GpuArray,3> rhsfab{rhs[0].array(mfi), + rhs[1].array(mfi), + rhs[2].array(mfi)}; + GpuArray,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 { diff --git a/Tests/LinearSolvers/CurlCurl/initProb_K.H b/Tests/LinearSolvers/CurlCurl/initProb_K.H index 0401b5c9b7a..66b7d2916cd 100644 --- a/Tests/LinearSolvers/CurlCurl/initProb_K.H +++ b/Tests/LinearSolvers/CurlCurl/initProb_K.H @@ -5,70 +5,102 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void actual_init_prob (int i, int j, int k, - amrex::GpuArray,AMREX_SPACEDIM> const& rhs, - amrex::GpuArray,AMREX_SPACEDIM> const& sol, + amrex::GpuArray,3> const& rhs, + amrex::GpuArray,3> const& sol, amrex::GpuArray const& problo, amrex::GpuArray const& dx, amrex::Real alpha, amrex::Real beta) { using namespace amrex; - constexpr amrex::Real tpi = Real(2.)*amrex::Math::pi(); - constexpr amrex::Real fpi = Real(4.)*amrex::Math::pi(); + constexpr Real pi = amrex::Math::pi(); 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); } }