1/*
2 * petsc_solver_unit_tests.hpp
3 *
4 * Created on: Jun 19, 2017
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_UNIT_TESTS_CPP_
9#define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_UNIT_TESTS_CPP_
10
11#ifdef HAVE_PETSC
12
13#define BOOST_TEST_DYN_LINK
14#include <boost/test/unit_test.hpp>
15
16#include "petsc_solver_report_unit_tests.hpp"
17
18#include "Grid/grid_dist_id.hpp"
19#include "Matrix/SparseMatrix.hpp"
20#include "Vector/Vector.hpp"
21#include "FiniteDifference/Laplacian.hpp"
22#include "FiniteDifference/FDScheme.hpp"
23#include "Solvers/petsc_solver.hpp"
24
25
26BOOST_AUTO_TEST_SUITE( mg_solvers_test )
27
28
29BOOST_AUTO_TEST_CASE( laplacian_3D_int_zero_mg )
30{
31 constexpr unsigned int phi = 0;
32 typedef Field<phi,poisson_nn_helm> phi_f;
33
34 Vcluster & v_cl = create_vcluster();
35 if (v_cl.getProcessingUnits() != 3)
36 return;
37
38 Ghost<3,long int> g(2);
39 Ghost<3,long int> stencil_max(2);
40 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
41
42 periodicity<3> p({PERIODIC,PERIODIC,PERIODIC});
43
44 typedef Lap<phi_f,poisson_nn_helm,CENTRAL_SYM> poisson;
45 grid_dist_id<3,float,aggregate<float>> psi({64,64,64},domain,g,p);
46 grid_dist_id<3,float,aggregate<float>> psi2(psi.getDecomposition(),{64,64,64},g);
47
48 // Fill B
49
50 size_t center_x = psi.size(0) / 2;
51 size_t center_y = psi.size(1) / 2;
52 size_t center_z = psi.size(2) / 2;
53 auto it = psi.getDomainIterator();
54
55 while (it.isNext())
56 {
57 auto key = it.get();
58 auto gkey = it.getGKey(key);
59
60 float sx = (float)(gkey.get(0))-center_x;
61 float sy = (float)(gkey.get(1))-center_y;
62 float sz = (float)(gkey.get(2))-center_z;
63
64 float gs = 100.0*exp(-((sx*sx)+(sy*sy)+(sz*sz))/100.0);
65
66 psi.get<0>(key) = sin(2*M_PI*sx/psi.size(0))*sin(2*M_PI*sy/psi.size(1))*sin(2*M_PI*sz/psi.size(2))*gs;
67 psi2.get<0>(key) = sin(2*M_PI*sx/psi.size(0))*sin(2*M_PI*sy/psi.size(1))*sin(2*M_PI*sz/psi.size(2))*gs;
68
69 ++it;
70 }
71
72 FDScheme<poisson_nn_helm> fd(stencil_max, domain, psi);
73
74 fd.template impose_dit<0>(poisson(),psi,psi.getDomainIterator());
75
76 petsc_solver<double> solver;
77 solver.setSolver(KSPBCGS);
78 solver.setAbsTol(0.01);
79 solver.setMaxIter(500);
80
81 ////////
82
83 solver.setPreconditioner(PCHYPRE_BOOMERAMG);
84 solver.setPreconditionerAMG_nl(6);
85 solver.setPreconditionerAMG_maxit(3);
86 solver.setPreconditionerAMG_relax("SOR/Jacobi");
87 solver.setPreconditionerAMG_cycleType("V",6,6);
88 solver.setPreconditionerAMG_coarsen("HMIS");
89 solver.setPreconditionerAMG_coarsenNodalType(0);
90
91 timer tm_solve;
92 tm_solve.start();
93 auto x_ = solver.solve(fd.getA(),fd.getB());
94 tm_solve.stop();
95
96 fd.template copy<phi>(x_,psi);
97 psi.write("AMG_psi");
98
99 #ifdef HAVE_OSX
100 bool check = compare("AMG_psi_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk","test/AMG_psi_" + std::to_string(v_cl.getProcessUnitID()) + "_test_osx.vtk");
101 #else
102 #if __GNUC__ == 6
103 bool check = compare("AMG_psi_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk","test/AMG_psi_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC6.vtk");
104 #endif
105 #if __GNUC__ == 4
106 bool check = compare("AMG_psi_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk","test/AMG_psi_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC4.vtk");
107 #endif
108 #endif
109
110 BOOST_REQUIRE_EQUAL(check,true);
111
112
113 // Resolve
114 timer tm_solve2;
115 tm_solve2.start();
116 auto x2_ = solver.solve(fd.getB());
117 tm_solve2.stop();
118
119 fd.template copy<phi>(x_,psi2);
120 psi2.write("AMG_psi2");
121
122 #ifdef HAVE_OSX
123 check = compare("AMG_psi2_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk","test/AMG_psi2_" + std::to_string(v_cl.getProcessUnitID()) + "_test_osx.vtk");
124 #else
125 #if __GNUC__ == 6
126 check = compare("AMG_psi2_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk","test/AMG_psi2_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC6.vtk");
127 #endif
128 #if __GNUC__ == 4
129 check = compare("AMG_psi2_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk","test/AMG_psi2_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC4.vtk");
130 #endif
131 #endif
132 BOOST_REQUIRE_EQUAL(check,true);
133}
134
135BOOST_AUTO_TEST_SUITE_END()
136
137#endif
138
139#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_UNIT_TESTS_CPP_ */
140