1 /** @file gsPreconditioner_test.cpp
2
3 @brief Tests for various preconditioners for Poisson problems.
4
5 This Source Code Form is subject to the terms of the Mozilla Public
6 License, v. 2.0. If a copy of the MPL was not distributed with this
7 file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
9 Author(s): J. Sogn, S. Takacs
10 */
11
12 #include "gismo_unittest.h"
13 #include <gsSolver/gsSolverUtils.h>
14
15 using namespace gismo;
16
runPreconditionerTest(index_t testcase)17 void runPreconditionerTest( index_t testcase )
18 {
19 // Define Geometry
20 gsMultiPatch<> mp( *gsNurbsCreator<>::NurbsQuarterAnnulus() );
21
22 // Create mulibasis
23 gsMultiBasis<> mb(mp);
24
25 // Refine multibais
26 const index_t numRefine = 4;
27 for (int i = 0; i < numRefine; ++i)
28 mb.uniformRefine();
29
30 // Define Boundary conditions
31 gsConstantFunction<> one(1,mp.geoDim());
32 gsBoundaryConditions<> bc;
33 bc.addCondition( boundary::west, condition_type::neumann, &one );
34 bc.addCondition( boundary::east, condition_type::neumann, &one );
35 bc.addCondition( boundary::south, condition_type::neumann, &one );
36 bc.addCondition( boundary::north, condition_type::dirichlet, &one );
37
38 // Initilize Assembler and assemble
39 gsOptionList opt = gsAssembler<>::defaultOptions();
40 gsPoissonAssembler<> assembler(
41 mp,
42 mb,
43 bc,
44 one,
45 (dirichlet::strategy) opt.getInt("DirichletStrategy"),
46 (iFace::strategy) opt.getInt("InterfaceStrategy")
47 );
48
49 assembler.assemble();
50 gsSparseMatrix<> mat = assembler.matrix();
51 gsMatrix<> rhs = assembler.rhs();
52 gsMatrix<> sol;
53 sol.setRandom(rhs.rows(), rhs.cols());
54
55 if (testcase==0)
56 {
57 gsConjugateGradient<> solver(mat, makeJacobiOp(mat));
58 solver.setTolerance( 1.e-8 );
59 solver.setMaxIterations( 110 );
60 solver.solve(rhs,sol);
61 CHECK ( solver.error() <= solver.tolerance() );
62 }
63 else if (testcase==1)
64 {
65 gsConjugateGradient<> solver(mat, makeSymmetricGaussSeidelOp(mat));
66 solver.setTolerance( 1.e-8 );
67 solver.setMaxIterations( 35 );
68 solver.solve(rhs,sol);
69 CHECK ( solver.error() <= solver.tolerance() );
70 }
71 else if (testcase==2)
72 {
73 gsConjugateGradient<> solver(mat, gsPatchPreconditionersCreator<>::fastDiagonalizationOp(mb[0],bc));
74 solver.setTolerance( 1.e-8 );
75 solver.setMaxIterations( 35 );
76 solver.solve(rhs,sol);
77 CHECK ( solver.error() <= solver.tolerance() );
78 }
79 else if (testcase==3)
80 {
81 const real_t h = mb[0].getMinCellLength();
82 gsGenericAssembler<> gAssembler(mp,mb,opt,&bc);
83 mat += (1/(h*h)) * gAssembler.assembleMass();
84 gsConjugateGradient<> solver(mat, gsPatchPreconditionersCreator<>::subspaceCorrectedMassSmootherOp(mb[0],bc));
85 solver.setTolerance( 1.e-8 );
86 solver.setMaxIterations( 50 );
87 solver.solve(rhs,sol);
88 CHECK ( solver.error() <= solver.tolerance() );
89 }
90 }
91
92
SUITE(gsPreconditioner_test)93 SUITE(gsPreconditioner_test)
94 {
95
96 TEST(gsJacobiPreconditioner_test)
97 {
98 runPreconditionerTest(0);
99 }
100 TEST(gsSymmetricGaussSeidelPreconditioner_test)
101 {
102 runPreconditionerTest(1);
103 }
104 TEST(gsFastDiagonalizationPreconditioner_test)
105 {
106 runPreconditionerTest(2);
107 }
108 TEST(gsSubspaceCorrectedMassPreconditioner_test)
109 {
110 runPreconditionerTest(3);
111 }
112
113 TEST(gsPatchPreconditioner_stiff_test)
114 {
115 // Define Geometry
116 gsMultiPatch<> mp( *gsNurbsCreator<>::BSplineSquare() );
117
118 // Create mulibasis
119 gsMultiBasis<> mb(mp);
120
121 // Refine multibasis
122 mb.uniformRefine();
123 dynamic_cast< gsTensorBSplineBasis<2>& >(mb[0]).component(0).uniformRefine();
124
125 // Set degree
126 mb[0].setDegreePreservingMultiplicity(3);
127
128 // Define Boundary conditions
129 gsConstantFunction<> one(1,mp.geoDim());
130 gsBoundaryConditions<> bc;
131 bc.addCondition( boundary::west, condition_type::neumann, &one );
132 bc.addCondition( boundary::east, condition_type::neumann, &one );
133 bc.addCondition( boundary::south, condition_type::neumann, &one );
134 bc.addCondition( boundary::north, condition_type::dirichlet, &one );
135
136 // Initilize Assembler and assemble
137 gsOptionList opt = gsAssembler<>::defaultOptions();
138 gsPoissonAssembler<> assembler(
139 mp,
140 mb,
141 bc,
142 one,
143 (dirichlet::strategy) opt.getInt("DirichletStrategy"),
144 (iFace::strategy) opt.getInt("InterfaceStrategy")
145 );
146 assembler.assemble();
147 const gsSparseMatrix<>& stiff0 = assembler.matrix();
148
149 // Get stiffness matrix from gsPatchPreconditionersCreator
150 gsSparseMatrix<> stiff1 = gsPatchPreconditionersCreator<>::stiffnessMatrix(mb[0],bc,opt);
151 CHECK ( (stiff1-stiff0 ).norm() < 1/(real_t)(10000) );
152
153 // Get stiffness matrix operator from gsPatchPreconditionersCreator
154 gsMatrix<> stiff2;
155 gsPatchPreconditionersCreator<>::stiffnessMatrixOp(mb[0],bc,opt)->toMatrix(stiff2);
156 CHECK ( (stiff2-stiff0 ).norm() < 1/(real_t)(10000) );
157
158 // Get inverse of stiffness matrix from gsPatchPreconditionersCreator
159 gsLinearOperator<>::Ptr stiffInv = gsPatchPreconditionersCreator<>::fastDiagonalizationOp(mb[0],bc,opt);
160 gsMatrix<> result;
161 stiffInv->apply(stiff0,result);
162 for (index_t i=0; i<result.rows(); ++i) result(i,i)-=1;
163 CHECK ( result.norm() < 1/(real_t)(10000) );
164 }
165
166 TEST(gsPatchPreconditioner_mass_test)
167 {
168 // Define Geometry
169 gsMultiPatch<> mp( *gsNurbsCreator<>::BSplineSquare() );
170
171 // Create mulibasis
172 gsMultiBasis<> mb(mp);
173
174 // Refine multibasis
175 mb.uniformRefine();
176 dynamic_cast< gsTensorBSplineBasis<2>& >(mb[0]).component(0).uniformRefine();
177
178 // Set degree
179 mb[0].setDegreePreservingMultiplicity(3);
180
181 // Define Boundary conditions
182 gsConstantFunction<> one(1,mp.geoDim());
183 gsBoundaryConditions<> bc;
184 bc.addCondition( boundary::west, condition_type::neumann, &one );
185 bc.addCondition( boundary::east, condition_type::neumann, &one );
186 bc.addCondition( boundary::south, condition_type::neumann, &one );
187 bc.addCondition( boundary::north, condition_type::dirichlet, &one );
188
189 // Initilize Assembler and assemble
190 gsOptionList opt = gsAssembler<>::defaultOptions();
191 gsGenericAssembler<> assembler(
192 mp,
193 mb,
194 opt,
195 &bc
196 );
197 gsSparseMatrix<> mass0 = assembler.assembleMass();
198
199
200 // Get mass matrix from gsPatchPreconditionersCreator
201 gsSparseMatrix<> mass1 = gsPatchPreconditionersCreator<>::massMatrix(mb[0],bc,opt);
202 CHECK ( ( mass0-mass1 ).norm() < 1/(real_t)(10000) );
203
204 // Get mass matrix operator from gsPatchPreconditionersCreator
205 gsMatrix<> mass2;
206 gsPatchPreconditionersCreator<>::massMatrixOp(mb[0],bc,opt)->toMatrix(mass2);
207 CHECK ( ( mass0-mass2 ).norm() < 1/(real_t)(10000) );
208
209 // Get inverse of mass matrix from gsPatchPreconditionersCreator
210 gsLinearOperator<>::Ptr massInv = gsPatchPreconditionersCreator<>::massMatrixInvOp(mb[0],bc,opt);
211 gsMatrix<> result;
212 massInv->apply(mass0,result);
213 for (index_t i=0; i<result.rows(); ++i) result(i,i)-=1;
214 CHECK ( result.norm() < 1/(real_t)(10000) );
215 }
216
217 TEST(gsAdditiveOp_test)
218 {
219 gsSparseMatrix<real_t,RowMajor> t1(3,2);
220 t1(0,0)=1;
221 t1(1,1)=1;
222 gsSparseMatrix<real_t,RowMajor> t2(3,1);
223 t2(2,0)=1;
224
225 std::vector< gsSparseMatrix<real_t,RowMajor> > t;
226 t.reserve(2);
227 t.push_back(t1);
228 t.push_back(t2);
229
230 gsMatrix<> o1(2,2);
231 o1 << 1,0, 0,2;
232
233 gsMatrix<> o2(1,1);
234 o2 << 3;
235
236 std::vector< gsLinearOperator<>::Ptr > o;
237 o.reserve(2);
238 o.push_back(makeMatrixOp(o1));
239 o.push_back(makeMatrixOp(o2));
240
241 gsMatrix<> in(3,1);
242 in << 9,8,7;
243 gsMatrix<> out(3,1);
244 out << 9,16,21;
245
246 {
247 gsSumOp<> s;
248 for (size_t i=0; i<t.size(); ++i)
249 s.addOperator(
250 gsProductOp<>::make(
251 makeMatrixOp(t[i].transpose()),
252 o[i],
253 makeMatrixOp(t[i])
254 )
255 );
256 gsMatrix<> res;
257 s.apply( in, res );
258 CHECK ( (res-out).norm() < 1/(real_t)(10000) );
259 }
260
261 {
262 gsAdditiveOp<> a(t,o);
263 gsMatrix<> res;
264 a.apply( in, res );
265 CHECK ( (res-out).norm() < 1/(real_t)(10000) );
266 }
267 }
268
269
270 }
271