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