Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
LinearSolver_def.hpp
1#ifndef LinearSolver_DEF_hpp
2#define LinearSolver_DEF_hpp
3#include "LinearSolver_decl.hpp"
12namespace FEDD {
13
14
15template<class SC,class LO,class GO,class NO>
16LinearSolver<SC,LO,GO,NO>::LinearSolver(){}
17
18
19template<class SC,class LO,class GO,class NO>
20LinearSolver<SC,LO,GO,NO>::~LinearSolver(){}
21
22template<class SC,class LO,class GO,class NO>
23int LinearSolver<SC,LO,GO,NO>::solve(Problem_Type* problem, BlockMultiVectorPtr_Type rhs, std::string type ){
24
25 int its=0;
26 if (!type.compare("Monolithic") || !type.compare("MonolithicConstPrec"))
27 its = solveMonolithic( problem, rhs, type );
28 else if (!type.compare("Teko")){
29#ifdef FEDD_HAVE_TEKO
30// its = solveTeko( problem, rhs );
31 its = solveBlock( problem, rhs, "Teko" );
32#else
33 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Teko not found! Build Trilinos with Teko.");
34#endif
35 }
36 else if (type=="Diagonal" || type=="Triangular" || type=="PCD" || type=="LSC"){
37 its = solveBlock( problem, rhs, type );
38 }
39 else
40 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown solver type.");
41
42 return its;
43}
44
45template<class SC,class LO,class GO,class NO>
46int LinearSolver<SC,LO,GO,NO>::solve(TimeProblem_Type* problem, BlockMultiVectorPtr_Type rhs, std::string type ){
47
48 int its=0;
49
50 if (!type.compare("Monolithic"))
51 its = solveMonolithic( problem, rhs );
52 else if (!type.compare("Teko")){
53#ifdef FEDD_HAVE_TEKO
54// its = solveTeko( problem, rhs );
55 its = solveBlock( problem, rhs, "Teko" );
56#else
57 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Teko not found! Build Trilinos with Teko.");
58#endif
59 }
60 else if(!type.compare("FaCSI") || type == "FaCSI-Teko" )
61 its = solveBlock( problem, rhs, type );
62 else if (type=="Diagonal" || type=="Triangular" || type=="PCD" || type=="LSC")
63 its = solveBlock( problem, rhs, type );
64 else
65 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown solver type.");
66
67 return its;
68}
69
70template<class SC,class LO,class GO,class NO>
71int LinearSolver<SC,LO,GO,NO>::solveMonolithic(Problem_Type* problem, BlockMultiVectorPtr_Type rhs, std::string type ){
72
73 bool verbose(problem->getVerbose());
74 int its=0;
75 if (problem->getParameterList()->get("Zero Initial Guess",true)) {
76 problem->getSolution()->putScalar(0.);
77 }
78 Teuchos::RCP<Thyra::MultiVectorBase<SC> >thyraX = problem->getSolution()->getThyraMultiVector();
79
80 Teuchos::RCP<const Thyra::MultiVectorBase<SC> > thyraB;
81 if (rhs.is_null())
82 thyraB = problem->getRhs()->getThyraMultiVectorConst();
83 else
84 thyraB = rhs->getThyraMultiVectorConst();
85
86 ParameterListPtr_Type pListThyraSolver = sublist( problem->getParameterList(), "ThyraSolver" );
87
88 pListThyraSolver->setParameters( problem->getParameterList()->sublist("ThyraPreconditioner") );
89
90 problem->getLinearSolverBuilder()->setParameterList(pListThyraSolver);
91 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > lowsFactory = problem->getLinearSolverBuilder()->createLinearSolveStrategy("");
92
93 bool iterativeSolve = !pListThyraSolver->get("Linear Solver Type", "Belos").compare("Belos");
94 if (iterativeSolve && (type != "MonolithicConstPrec" || problem->getPreconditioner()->getThyraPrec().is_null()))
95 problem->setupPreconditioner("Monolithic");
96
97 if (!pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive").compare("Multiplicative")) {
98 pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
99
100 Teuchos::RCP<const Thyra::LinearOpBase<SC> > thyra_linOp = problem->getPreconditioner()->getThyraPrec()->getUnspecifiedPrecOp();
101 Thyra::apply( *thyra_linOp, Thyra::NOTRANS, *thyraB, thyraX.ptr() );
102 pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",false);
103 }
104
105 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
106
107 lowsFactory->setOStream(out);
108 lowsFactory->setVerbLevel(Teuchos::VERB_HIGH);
109
110 Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > solver = lowsFactory->createOp();
111// Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > solver = linearOpWithSolve(*lowsFactory, problem->getSystem()->getThyraLinOp());
112 ThyraLinOpConstPtr_Type thyraMatrix = problem->getSystem()->getThyraLinOp();
113 if ( iterativeSolve ) {
114 ThyraPrecPtr_Type thyraPrec = problem->getPreconditioner()->getThyraPrec();
115 Thyra::initializePreconditionedOp<SC>(*lowsFactory, thyraMatrix, thyraPrec.getConst(), solver.ptr());
116 }
117 else{
118 Thyra::initializeOp<SC>(*lowsFactory, thyraMatrix, solver.ptr());
119 }
120
121 {
122 Thyra::SolveStatus<SC> status = Thyra::solve<SC>(*solver, Thyra::NOTRANS, *thyraB, thyraX.ptr());
123 if (verbose)
124 std::cout << status << std::endl;
125 if ( !pListThyraSolver->get("Linear Solver Type","Belos").compare("Belos") )
126 its = status.extraParameters->get("Belos/Iteration Count",0);
127 else
128 its = 0;
129
130 problem->getSolution()->fromThyraMultiVector(thyraX);
131 }
132
133 return its;
134}
135
136template<class SC,class LO,class GO,class NO>
137int LinearSolver<SC,LO,GO,NO>::solveMonolithic(TimeProblem_Type* timeProblem, BlockMultiVectorPtr_Type rhs){
138
139
140 bool verbose(timeProblem->getVerbose());
141 int its=0;
142 ProblemPtr_Type problem = timeProblem->getUnderlyingProblem();
143
144 if (problem->getParameterList()->get("Zero Initial Guess",true)) {
145 problem->getSolution()->putScalar(0.);
146 }
147
148 ParameterListPtr_Type pListThyraSolver = sublist( problem->getParameterList(), "ThyraSolver" );
149
150 Teuchos::RCP<Thyra::MultiVectorBase<SC> >thyraX = problem->getSolution()->getThyraMultiVector();
151 Teuchos::RCP<const Thyra::MultiVectorBase<SC> > thyraB;
152
153 if (rhs.is_null())
154 thyraB = problem->getRhs()->getThyraMultiVectorConst();
155 else
156 thyraB = rhs->getThyraMultiVectorConst();
157
158 pListThyraSolver->setParameters( problem->getParameterList()->sublist("ThyraPreconditioner") );
159
160 problem->getLinearSolverBuilder()->setParameterList(pListThyraSolver);
161 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > lowsFactory = problem->getLinearSolverBuilder()->createLinearSolveStrategy("");
162
163
164 problem->setupPreconditioner( "Monolithic" );
165
166 if (!pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive").compare("Multiplicative")) {
167 pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
168
169 Teuchos::RCP<const Thyra::LinearOpBase<SC> > thyra_linOp = problem->getPreconditioner()->getThyraPrec()->getUnspecifiedPrecOp();
170 Thyra::apply( *thyra_linOp, Thyra::NOTRANS, *thyraB, thyraX.ptr() );
171 pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",false);
172 }
173
174 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
175
176 lowsFactory->setOStream(out);
177 lowsFactory->setVerbLevel(Teuchos::VERB_HIGH);
178
179 Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > solver = lowsFactory->createOp();
180
181 ThyraLinOpConstPtr_Type thyraMatrix = timeProblem->getSystemCombined()->getThyraLinOp();
182
183 if ( !pListThyraSolver->get("Linear Solver Type","Belos").compare("Belos") ) {
184 ThyraPrecPtr_Type thyraPrec = problem->getPreconditioner()->getThyraPrec();
185 Thyra::initializePreconditionedOp<SC>(*lowsFactory, thyraMatrix, thyraPrec.getConst(), solver.ptr());
186 }
187 else{
188 Thyra::initializeOp<SC>(*lowsFactory, thyraMatrix, solver.ptr());
189 }
190 {
191 Thyra::SolveStatus<SC> status = Thyra::solve<SC>(*solver, Thyra::NOTRANS, *thyraB, thyraX.ptr());
192 if (verbose)
193 std::cout << status << std::endl;
194 problem->getSolution()->fromThyraMultiVector(thyraX);
195 // problem->getSolution()->writeMM("solution");
196 if ( !pListThyraSolver->get("Linear Solver Type","Belos").compare("Belos") ){
197 its = status.extraParameters->get("Belos/Iteration Count",0);
198 double achievedTol = status.extraParameters->get("Belos/Achieved Tolerance",-1.);
199 }
200 else
201 its = 0;
202 }
203 return its;
204}
205
206template<class SC,class LO,class GO,class NO>
207int LinearSolver<SC,LO,GO,NO>::solveBlock(Problem_Type* problem, BlockMultiVectorPtr_Type rhs, std::string type ){
208 typedef Thyra::DefaultZeroLinearOp<SC> ZeroOp_Type;
209 typedef Teuchos::RCP<ZeroOp_Type> ZeroOpPtr_Type;
210
211 int rank = problem->getComm()->getRank();
212 bool verbose(problem->getVerbose());
213 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
214
215 int its=0;
216
217 if (problem->getParameterList()->get("Zero Initial Guess",true)) {
218 problem->getSolution()->putScalar(0.);
219 }
220
221
222 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > thyraX = problem->getSolution()->getProdThyraMultiVector();
223
224 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > thyraRHS;
225 if ( rhs.is_null() )
226 thyraRHS = problem->getRhs()->getProdThyraMultiVector();
227 else
228 thyraRHS = rhs->getProdThyraMultiVector();
229
230 ParameterListPtr_Type pListThyraSolver = sublist( problem->getParameterList(), "ThyraSolver" );
231
232 problem->getLinearSolverBuilder()->setParameterList(pListThyraSolver);
233 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > lowsFactory = problem->getLinearSolverBuilder()->createLinearSolveStrategy("");
234
235 problem->setupPreconditioner( type );
236
237 lowsFactory->setOStream(out);
238 lowsFactory->setVerbLevel(Teuchos::VERB_HIGH);
239
240 Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > solver = lowsFactory->createOp();
241
242 ThyraPrecPtr_Type thyraPrec = problem->getPreconditioner()->getThyraPrec();
243
244// BlockMatrixPtr_Type system = problem->getSystem();
245// ThyraLinOpBlockConstPtr_Type thyraMatrixBlocks = system->getThyraLinBlockOp();
246// ThyraLinOpBlockPtr_Type thyraMatrixBlocksNonConst = Teuchos::rcp_const_cast<ThyraLinOpBlock_Type>( thyraMatrixBlocks );
247//
248// for (int i=0; i<system->size(); i++) {
249// for (int j=0; j<system->size(); j++) {
250// if ( !thyraMatrixBlocks->blockExists(i,j) ) {
251// ZeroOpPtr_Type zeroOp =
252// Teuchos::rcp( new ZeroOp_Type( thyraRHS->getMultiVectorBlock(i)->range(),
253// thyraX->getMultiVectorBlock(j)->range() ) ); //(range , domain)
254// thyraMatrixBlocksNonConst->getNonconstBlock(i,j) = zeroOp;
255// }
256// }
257// }
258
259
260
261// BlockMatrixPtr_Type system = problem->getSystem();
262// for (int i=0; i<system->size(); i++) {
263// for (int j=0; j<system->size(); j++) {
264// if (!system->blockExists(1,1)){
265// MatrixPtr_Type dummy = Teuchos::rcp( new Matrix_Type( system->getBlock(i,j)->getMap(), 0 ) );
266// dummy->fillComplete( problem->getSolution()->getBlock( j )->getMap(), problem->getRhs()->getBlock( i )->getMap() ); //(domain, range)
267// system->addBlock( dummy, i, j );
268// }
269// }
270// }
271
272 ThyraLinOpConstPtr_Type thyraMatrix = problem->getSystem()->getThyraLinBlockOp();
273 Thyra::initializePreconditionedOp<SC>(*lowsFactory, thyraMatrix, thyraPrec.getConst(), solver.ptr());
274 {
275
276 Thyra::SolveStatus<SC> status = Thyra::solve<SC>(*solver, Thyra::NOTRANS, *thyraRHS, thyraX.ptr());
277 if (verbose)
278 std::cout << status << std::endl;
279
280 its = status.extraParameters->get("Belos/Iteration Count",0);
281
282 problem->getSolution()->fromThyraProdMultiVector( thyraX );
283 }
284
285 return its;
286}
287
288template<class SC,class LO,class GO,class NO>
289int LinearSolver<SC,LO,GO,NO>::solveBlock(TimeProblem_Type* timeProblem, BlockMultiVectorPtr_Type rhs, std::string type ){
290 typedef Thyra::DefaultZeroLinearOp<SC> ZeroOp_Type;
291 typedef Teuchos::RCP<ZeroOp_Type> ZeroOpPtr_Type;
292
293 int rank = timeProblem->getComm()->getRank();
294 bool verbose(timeProblem->getVerbose());
295 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
296
297 int its=0;
298
299 ProblemPtr_Type problem = timeProblem->getUnderlyingProblem();
300 if (problem->getParameterList()->get("Zero Initial Guess",true)) {
301 problem->getSolution()->putScalar(0.);
302 }
303
304 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > thyraX = problem->getSolution()->getProdThyraMultiVector();
305
306 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > thyraRHS;
307 if ( rhs.is_null() )
308 thyraRHS = problem->getRhs()->getProdThyraMultiVector();
309 else
310 thyraRHS = rhs->getProdThyraMultiVector();
311
312 ParameterListPtr_Type pListThyraSolver = sublist( problem->getParameterList(), "ThyraSolver" );
313
314
315 problem->getLinearSolverBuilder()->setParameterList(pListThyraSolver);
316 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > lowsFactory = problem->getLinearSolverBuilder()->createLinearSolveStrategy("");
317
318 problem->setupPreconditioner( type );
319
320 lowsFactory->setOStream(out);
321 lowsFactory->setVerbLevel(Teuchos::VERB_HIGH);
322
323 Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > solver = lowsFactory->createOp();
324
325 ThyraPrecPtr_Type thyraPrec = problem->getPreconditioner()->getThyraPrec();
326
327// BlockMatrixPtr_Type system = timeProblem->getSystemCombined();
328// ThyraLinOpBlockConstPtr_Type thyraMatrixBlocks = system->getThyraLinBlockOp();
329// ThyraLinOpBlockPtr_Type thyraMatrixBlocksNonConst = Teuchos::rcp_const_cast<ThyraLinOpBlock_Type>( thyraMatrixBlocks );
330//
331// for (int i=0; i<system->size(); i++) {
332// for (int j=0; j<system->size(); j++) {
333// if ( !thyraMatrixBlocks->blockExists(i,j) ) {
334// ZeroOpPtr_Type zeroOp =
335// Teuchos::rcp( new ZeroOp_Type( thyraRHS->getMultiVectorBlock(i)->range(),
336// thyraX->getMultiVectorBlock(j)->range() ) ); //(range , domain)
337// thyraMatrixBlocksNonConst->getNonconstBlock(i,j) = zeroOp;
338// }
339// }
340// }
341
342 BlockMatrixPtr_Type system = timeProblem->getSystemCombined();
343// for (int i=0; i<system->size(); i++) {
344// for (int j=0; j<system->size(); j++) {
345// if (!system->blockExists(1,1)){
346// MatrixPtr_Type dummy;
347// dummy.reset( new Matrix_Type( system->getBlock(i,j)->getMap(), 0 ) );
348// dummy->fillComplete( problem->getSolution()->getBlock( j )->getMap(), problem->getRhs()->getBlock( i )->getMap() ); //(domain, range)
349// system->addBlock( dummy, i, j );
350// }
351// }
352// }
353
354 ThyraLinOpConstPtr_Type thyraMatrix = timeProblem->getSystemCombined()->getThyraLinBlockOp();
355// ThyraLinOpBlockConstPtr_Type thyraMatrixBlock = timeProblem->getSystemCombined()->getThyraLinBlockOp();
356 Thyra::initializePreconditionedOp<SC>(*lowsFactory, thyraMatrix, thyraPrec.getConst(), solver.ptr());
357 {
358
359 Thyra::SolveStatus<SC> status = Thyra::solve<SC>(*solver, Thyra::NOTRANS, *thyraRHS, thyraX.ptr());
360 if (verbose)
361 std::cout << status << std::endl;
362
363 its = status.extraParameters->get("Belos/Iteration Count",0);
364
365 problem->getSolution()->fromThyraProdMultiVector( thyraX );
366 }
367
368 return its;
369}
370
371// #ifdef FEDD_HAVE_TEKO
372// template<class SC,class LO,class GO,class NO>
373// int LinearSolver<SC,LO,GO,NO>::solveTeko(Problem_Type* problem, BlockMultiVectorPtr_Type rhs ){
374
375// Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
376
377// bool verbose(problem->getVerbose());
378// int its=0;
379// if (problem->getParameterList()->get("Zero Initial Guess",true)) {
380// problem->getSolution()->putScalar(0.);
381// }
382// // typedef Teuchos::RCP< Thyra::ProductMultiVectorBase< double > > Teko::BlockedMultiVector
383// // convert them to teko compatible sub vectors
384// Teko::MultiVector x0_th = problem->getSolution()->getBlockNonConst(0)->getThyraMultiVector();
385// Teko::MultiVector x1_th = problem->getSolution()->getBlockNonConst(1)->getThyraMultiVector();
386
387// Teko::MultiVector b0_th;
388// Teko::MultiVector b1_th;
389// if (rhs.is_null()){
390// b0_th = problem->getRhs()->getBlockNonConst(0)->getThyraMultiVector();
391// b1_th = problem->getRhs()->getBlockNonConst(1)->getThyraMultiVector();
392// }
393// else{
394// b0_th = rhs->getBlockNonConst(0)->getThyraMultiVector();
395// b1_th = rhs->getBlockNonConst(1)->getThyraMultiVector();
396// }
397
398// std::vector<Teko::MultiVector> x_vec; x_vec.push_back(x0_th); x_vec.push_back(x1_th);
399// std::vector<Teko::MultiVector> b_vec; b_vec.push_back(b0_th); b_vec.push_back(b1_th);
400
401// Teko::MultiVector x = Teko::buildBlockedMultiVector(x_vec); // these will be used in the Teko solve
402// Teko::MultiVector b = Teko::buildBlockedMultiVector(b_vec);
403
404// ParameterListPtr_Type pListThyraSolver = sublist( problem->getParameterList(), "ThyraSolver" );
405
406// // pListThyraSolver->setParameters( problem->getParameterList()->sublist("ThyraPreconditioner") );
407
408// problem->getLinearSolverBuilder()->setParameterList(pListThyraSolver);
409// Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > lowsFactory = problem->getLinearSolverBuilder()->createLinearSolveStrategy("");
410
411// problem->setupPreconditioner( "Teko" );
412
413// // if (!pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive").compare("Multiplicative")) {
414// // pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
415// //
416// // Teuchos::RCP<const Thyra::LinearOpBase<SC> > thyra_linOp = problem->getPreconditioner()->getThyraPrec()->getUnspecifiedPrecOp();
417// // Thyra::apply( *thyra_linOp, Thyra::NOTRANS, *thyraB, thyraX.ptr() );
418// // pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",false);
419// // }
420
421
422// lowsFactory->setOStream(out);
423// // lowsFactory->setVerbLevel(Teuchos::VERB_EXTREME);
424
425// Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > solver = lowsFactory->createOp();
426
427// ThyraPrecPtr_Type thyraPrec = problem->getPreconditioner()->getThyraPrec();
428
429// ThyraLinOpConstPtr_Type thyraMatrix = problem->getPreconditioner()->getTekoOp();
430
431// Thyra::initializePreconditionedOp<SC>(*lowsFactory, thyraMatrix, thyraPrec.getConst(), solver.ptr());
432
433// {
434// Thyra::SolveStatus<SC> status = Thyra::solve<SC>(*solver, Thyra::NOTRANS, *b, x.ptr());
435// if (verbose)
436// std::cout << status << std::endl;
437
438// its = status.extraParameters->get("Belos/Iteration Count",0);
439// Teuchos::RCP< Thyra::ProductMultiVectorBase< double > > xBlocked = Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase< double > >( x );
440// problem->getSolution()->getBlockNonConst(0)->fromThyraMultiVector( xBlocked->getNonconstMultiVectorBlock( 0 ) );
441// problem->getSolution()->getBlockNonConst(1)->fromThyraMultiVector( xBlocked->getNonconstMultiVectorBlock( 1 ) );
442// }
443
444// return its;
445// }
446
447// template<class SC,class LO,class GO,class NO>
448// int LinearSolver<SC,LO,GO,NO>::solveTeko(TimeProblem_Type* timeProblem, BlockMultiVectorPtr_Type rhs ){
449
450// bool verbose(timeProblem->getVerbose());
451// int its=0;
452
453// ProblemPtr_Type problem = timeProblem->getUnderlyingProblem();
454// if (problem->getParameterList()->get("Zero Initial Guess",true)) {
455// problem->getSolution()->putScalar(0.);
456// }
457// // typedef Teuchos::RCP< Thyra::ProductMultiVectorBase< double > > Teko::BlockedMultiVector
458// // convert them to teko compatible sub vectors
459// Teko::MultiVector x0_th = problem->getSolution()->getBlockNonConst(0)->getThyraMultiVector();
460// Teko::MultiVector x1_th = problem->getSolution()->getBlockNonConst(1)->getThyraMultiVector();
461
462// Teko::MultiVector b0_th;
463// Teko::MultiVector b1_th;
464// if (rhs.is_null()){
465// b0_th = problem->getRhs()->getBlockNonConst(0)->getThyraMultiVector();
466// b1_th = problem->getRhs()->getBlockNonConst(1)->getThyraMultiVector();
467// }
468// else{
469// b0_th = rhs->getBlockNonConst(0)->getThyraMultiVector();
470// b1_th = rhs->getBlockNonConst(1)->getThyraMultiVector();
471// }
472
473// std::vector<Teko::MultiVector> x_vec; x_vec.push_back(x0_th); x_vec.push_back(x1_th);
474// std::vector<Teko::MultiVector> b_vec; b_vec.push_back(b0_th); b_vec.push_back(b1_th);
475
476// Teko::MultiVector x = Teko::buildBlockedMultiVector(x_vec); // these will be used in the Teko solve
477// Teko::MultiVector b = Teko::buildBlockedMultiVector(b_vec);
478
479// ParameterListPtr_Type pListThyraSolver = sublist( problem->getParameterList(), "ThyraSolver" );
480
481// // pListThyraSolver->setParameters( problem->getParameterList()->sublist("ThyraPreconditioner") );
482
483// problem->getLinearSolverBuilder()->setParameterList(pListThyraSolver);
484// Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > lowsFactory = problem->getLinearSolverBuilder()->createLinearSolveStrategy("");
485
486// problem->setupPreconditioner( "Teko" );
487
488// // if (!pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive").compare("Multiplicative")) {
489// // pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
490// //
491// // Teuchos::RCP<const Thyra::LinearOpBase<SC> > thyra_linOp = problem->getPreconditioner()->getThyraPrec()->getUnspecifiedPrecOp();
492// // Thyra::apply( *thyra_linOp, Thyra::NOTRANS, *thyraB, thyraX.ptr() );
493// // pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",false);
494// // }
495
496// Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
497
498// lowsFactory->setOStream(out);
499// lowsFactory->setVerbLevel(Teuchos::VERB_HIGH);
500
501// Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > solver = lowsFactory->createOp();
502
503// ThyraPrecPtr_Type thyraPrec = problem->getPreconditioner()->getThyraPrec();
504
505// ThyraLinOpConstPtr_Type thyraMatrix = problem->getPreconditioner()->getTekoOp();
506
507// Thyra::initializePreconditionedOp<SC>(*lowsFactory, thyraMatrix, thyraPrec.getConst(), solver.ptr());
508
509// {
510// Thyra::SolveStatus<SC> status = Thyra::solve<SC>(*solver, Thyra::NOTRANS, *b, x.ptr());
511// if (verbose)
512// std::cout << status << std::endl;
513
514// its = status.extraParameters->get("Belos/Iteration Count",0);
515// Teuchos::RCP< Thyra::ProductMultiVectorBase< double > > xBlocked = Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase< double > >( x );
516// problem->getSolution()->getBlockNonConst(0)->fromThyraMultiVector( xBlocked->getNonconstMultiVectorBlock( 0 ) );
517// problem->getSolution()->getBlockNonConst(1)->fromThyraMultiVector( xBlocked->getNonconstMultiVectorBlock( 1 ) );
518// }
519
520// return its;
521// }
522// #endif
523}
524#endif
FEDD::LinearSolver::solveMonolithic
int solveMonolithic(Problem_Type *problem, BlockMultiVectorPtr_Type rhs, std::string type)
Solve linear/linearized problem monolithicly.
Definition LinearSolver_def.hpp:71
FEDD::LinearSolver::solveBlock
int solveBlock(Problem_Type *problem, BlockMultiVectorPtr_Type rhs, std::string precType)
In case of a block system, solve block is called. It works also for teko.
Definition LinearSolver_def.hpp:207
FEDD::LinearSolver::solve
int solve(Problem_Type *problem, BlockMultiVectorPtr_Type rhs, std::string type="Monolithic")
Call to solve a linear/linearized problem with right-hand side rhs. Depending on 'type' solveMonolith...
Definition LinearSolver_def.hpp:23
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:33