Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
NonLinearProblem_def.hpp
1#ifndef NONLINEARPROBLEM_DEF_hpp
2#define NONLINEARPROBLEM_DEF_hpp
3
4#include "Problem.hpp"
5#include "feddlib/core/LinearAlgebra/BlockMap.hpp"
6#include "feddlib/core/LinearAlgebra/BlockMatrix.hpp"
7#include "feddlib/core/FE/Domain.hpp"
8#include "feddlib/problems/Solver/Preconditioner.hpp"
9
18
19namespace FEDD
20{
21
22 template <class SC, class LO, class GO, class NO>
23 NonLinearProblem<SC, LO, GO, NO>::NonLinearProblem(CommConstPtr_Type comm) : Problem<SC, LO, GO, NO>(comm),
24 previousSolution_(),
25 residualVec_(),
26 nonLinearTolerance_(1.e-6),
27 coeff_(0)
28 {
29 }
30
31 template <class SC, class LO, class GO, class NO>
32 NonLinearProblem<SC, LO, GO, NO>::NonLinearProblem(ParameterListPtr_Type &parameterList, CommConstPtr_Type comm) : Problem<SC, LO, GO, NO>(parameterList, comm),
33 previousSolution_(),
34 residualVec_(),
35 nonLinearTolerance_(1.e-6),
36 coeff_(0)
37 {
38 }
39 template <class SC, class LO, class GO, class NO>
40 NonLinearProblem<SC, LO, GO, NO>::~NonLinearProblem()
41 {
42 }
43
44 template <class SC, class LO, class GO, class NO>
45 void NonLinearProblem<SC, LO, GO, NO>::initializeProblem(int nmbVectors)
46 {
47 this->system_.reset(new BlockMatrix_Type(1));
48
49 this->initializeVectors(nmbVectors);
50
51 this->initializeVectorsNonLinear(nmbVectors);
52
53 // Init ThyraVectorSpcaes for NOX.
54 this->initVectorSpaces();
55
56 this->newtonStep_=0;
57 }
58
59 template <class SC, class LO, class GO, class NO>
60 void NonLinearProblem<SC, LO, GO, NO>::initializeVectorsNonLinear(int nmbVectors)
61 {
62
63 UN size = this->domainPtr_vec_.size();
64 this->previousSolution_.reset(new BlockMultiVector_Type(size));
65 this->residualVec_.reset(new BlockMultiVector_Type(size));
66
67 for (UN i = 0; i < size; i++)
68 {
69 if (this->dofsPerNode_vec_[i] > 1)
70 {
71 MapConstPtr_Type map = this->domainPtr_vec_[i]->getMapVecFieldUnique();
72 MultiVectorPtr_Type prevSolutionPart = Teuchos::rcp(new MultiVector_Type(map));
73 this->previousSolution_->addBlock(prevSolutionPart, i);
74 MultiVectorPtr_Type residualPart = Teuchos::rcp(new MultiVector_Type(map));
75 this->residualVec_->addBlock(residualPart, i);
76 }
77 else
78 {
79 MapConstPtr_Type map = this->domainPtr_vec_[i]->getMapUnique();
80 MultiVectorPtr_Type prevSolutionPart = Teuchos::rcp(new MultiVector_Type(map));
81 this->previousSolution_->addBlock(prevSolutionPart, i);
82 MultiVectorPtr_Type residualPart = Teuchos::rcp(new MultiVector_Type(map));
83 this->residualVec_->addBlock(residualPart, i);
84 }
85 }
86
87 this->residualVec_->putScalar(0.);
88 }
89
90 template <class SC, class LO, class GO, class NO>
91 void NonLinearProblem<SC, LO, GO, NO>::infoNonlinProblem()
92 {
93
94 bool verbose(this->comm_->getRank() == 0);
95 if (verbose)
96 {
97 std::cout << "\t ### ### ###" << std::endl;
98 std::cout << "\t ### Nonlinear Problem Information ###" << std::endl;
99 std::cout << "\t ### Linearization: " << this->parameterList_->sublist("General").get("Linearization", "default") << std::endl;
100 std::cout << "\t ### Relative tol: " << this->parameterList_->sublist("Parameter").get("relNonLinTol", 1.e-6) << "\t absolute tol: " << this->parameterList_->sublist("Parameter").get("absNonLinTol", 1.e-4) << "(not used for Newton or Fixed-Point)" << std::endl;
101 }
102 }
103
104 template <class SC, class LO, class GO, class NO>
105 void NonLinearProblem<SC, LO, GO, NO>::calculateNonLinResidualVec(SmallMatrix<double> &coeff, std::string type, double time, BlockMatrixPtr_Type systemMass)
106 {
107
108 coeff_ = coeff;
109 this->calculateNonLinResidualVec(type, time);
110 }
111
112 template <class SC, class LO, class GO, class NO>
113 void NonLinearProblem<SC, LO, GO, NO>::reAssembleAndFill(BlockMatrixPtr_Type bMat, std::string type)
114 {
115
116 this->assemble(type);
117 TEUCHOS_TEST_FOR_EXCEPTION(bMat->size() != this->system_->size(), std::logic_error, "Sizes of BlockMatrices are differen. reAssembleAndFill(...)");
118
119 for (int i = 0; i < this->system_->size(); i++)
120 {
121 for (int j = 0; j < this->system_->size(); j++)
122 {
123 if (this->system_->blockExists(i, j))
124 {
125 // MatrixPtr_Type mat = Teuchos::rcp(new Matrix_Type ( this->system_->getBlock( i, j ) ) );
126 bMat->addBlock(this->system_->getBlock(i, j), i, j);
127 }
128 }
129 }
130 }
131
132 template <class SC, class LO, class GO, class NO>
133 double NonLinearProblem<SC, LO, GO, NO>::calculateResidualNorm() const
134 {
135
136 Teuchos::Array<SC> residual(1);
137 residualVec_->norm2(residual());
138 TEUCHOS_TEST_FOR_EXCEPTION(residual.size() != 1, std::logic_error, "We need to change the code for numVectors>1.");
139 return residual[0];
140 }
141
142 template <class SC, class LO, class GO, class NO>
143 int NonLinearProblem<SC, LO, GO, NO>::solveUpdate()
144 {
145
146 // solution COPY!
147 *previousSolution_ = *this->solution_;
148
149 int its = this->solve(residualVec_);
150
151 return its;
152 }
153
154 template <class SC, class LO, class GO, class NO>
155 int NonLinearProblem<SC, LO, GO, NO>::solveAndUpdate(const std::string &criterion, double &criterionValue)
156 {
157 // BlockMatrixPtr_Type system
158 int its = solveUpdate();
159
160 if (criterion == "Update")
161 {
162 Teuchos::Array<SC> updateNorm(1);
163 this->solution_->norm2(updateNorm());
164 criterionValue = updateNorm[0];
165 }
166 this->solution_->update(1., *previousSolution_, 1.);
167
168 return its;
169 }
170
171 template <class SC, class LO, class GO, class NO>
172 typename NonLinearProblem<SC, LO, GO, NO>::BlockMultiVectorPtr_Type NonLinearProblem<SC, LO, GO, NO>::getResidualVector() const
173 {
174
175 return residualVec_;
176 }
177
178 template <class SC, class LO, class GO, class NO>
179 Thyra::ModelEvaluatorBase::InArgs<SC>
180 NonLinearProblem<SC, LO, GO, NO>::getNominalValues() const
181 {
182 return nominalValues_;
183 }
184
185 template <class SC, class LO, class GO, class NO>
186 Teuchos::RCP<const ::Thyra::VectorSpaceBase<SC>> NonLinearProblem<SC, LO, GO, NO>::get_x_space() const
187 {
188 return xSpace_;
189 }
190
191 template <class SC, class LO, class GO, class NO>
192 Teuchos::RCP<const ::Thyra::VectorSpaceBase<SC>> NonLinearProblem<SC, LO, GO, NO>::get_f_space() const
193 {
194 return fSpace_;
195 }
196
197 template <class SC, class LO, class GO, class NO>
198 Thyra::ModelEvaluatorBase::InArgs<SC>
199 NonLinearProblem<SC, LO, GO, NO>::createInArgs() const
200 {
201 return prototypeInArgs_;
202 }
203
204 // Private functions overridden from ModelEvaulatorDefaultBase
205
206 template <class SC, class LO, class GO, class NO>
207 Thyra::ModelEvaluatorBase::OutArgs<SC>
208 NonLinearProblem<SC, LO, GO, NO>::createOutArgsImpl() const
209 {
210 return prototypeOutArgs_;
211 }
212
213 template <class SC, class LO, class GO, class NO>
214 void NonLinearProblem<SC, LO, GO, NO>::initNOXParameters()
215 {
216
217 using Teuchos::RCP;
218 using Teuchos::rcp;
219 using ::Thyra::VectorBase;
220 typedef ::Thyra::ModelEvaluatorBase MEB;
221
222 MEB::InArgsSetup<SC> inArgs;
223 inArgs.setModelEvalDescription(this->description());
224 inArgs.setSupports(MEB::IN_ARG_x);
225 this->prototypeInArgs_ = inArgs;
226
227 MEB::OutArgsSetup<SC> outArgs;
228 outArgs.setModelEvalDescription(this->description());
229 outArgs.setSupports(MEB::OUT_ARG_f);
230 outArgs.setSupports(MEB::OUT_ARG_W_op);
231 outArgs.setSupports(MEB::OUT_ARG_W_prec);
232 this->prototypeOutArgs_ = outArgs;
233
234 this->nominalValues_ = inArgs;
235 }
236
237 template <class SC, class LO, class GO, class NO>
238 void NonLinearProblem<SC, LO, GO, NO>::initVectorSpaces()
239 {
240 this->initNOXParameters();
241
242 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method", "Monolithic");
243 if (!type.compare("Monolithic"))
244 initVectorSpacesMonolithic();
245 else if (!type.compare("Teko") || type == "FaCSI" || type == "FaCSI-Teko" || type == "FaCSI-Block" || type == "Diagonal" || type == "Triangular" || type == "PCD"|| type == "LSC")
246 initVectorSpacesBlock();
247 else
248 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Unkown preconditioner/solver type.");
249 }
250
251 template <class SC, class LO, class GO, class NO>
252 void NonLinearProblem<SC, LO, GO, NO>::initVectorSpacesMonolithic()
253 {
254
255 BlockMapPtr_Type map = Teuchos::rcp_const_cast<BlockMap_Type>(this->solution_->getMap());
256
257 // Teuchos::RCP<const XTpetra_Type> xTpetraMap = Teuchos::rcp_dynamic_cast<const XTpetra_Type>(map->getMergedMap()->getXpetraMap()->getMap());
258
259 TpetraMapConstPtr_Type tpetraMap = map->getMergedMap()->getTpetraMap();
260
261 this->xSpace_ = Thyra::createVectorSpace<SC, LO, GO, NO>(tpetraMap);
262 this->fSpace_ = Thyra::createVectorSpace<SC, LO, GO, NO>(tpetraMap);
263
264 typedef Teuchos::ScalarTraits<SC> ST;
265 x0_ = ::Thyra::createMember(this->xSpace_);
266 V_S(x0_.ptr(), ST::zero());
267
268 this->nominalValues_.set_x(x0_);
269 }
270
271 template <class SC, class LO, class GO, class NO>
272 void NonLinearProblem<SC, LO, GO, NO>::initVectorSpacesBlock()
273 {
274
275 BlockMapPtr_Type map = Teuchos::rcp_const_cast<BlockMap_Type>(this->solution_->getMap());
276
277 TpetraMap_Type tpetra_map;
278
279 Teuchos::Array<ThyraVecSpaceConstPtr_Type> vecSpaceArray(map->size());
280 for (int i = 0; i < map->size(); i++)
281 {
282 //Teuchos::RCP<const XTpetra_Type> xTpetraMap =
283 // Teuchos::rcp_dynamic_cast<const XTpetra_Type>(map->getBlock(i)->getXpetraMap()->getMap());
284 TpetraMapConstPtr_Type tpetraMap =map->getBlock(i)->getTpetraMap();
285 ThyraVecSpaceConstPtr_Type vecSpace = Thyra::createVectorSpace<SC, LO, GO, NO>(tpetraMap);
286 vecSpaceArray[i] = vecSpace;
287 }
288 this->xSpace_ = Teuchos::rcp(new Thyra::DefaultProductVectorSpace<SC>(vecSpaceArray()));
289 this->fSpace_ = Teuchos::rcp(new Thyra::DefaultProductVectorSpace<SC>(vecSpaceArray()));
290
291 typedef Teuchos::ScalarTraits<SC> ST;
292 x0_ = ::Thyra::createMember(this->xSpace_);
293 V_S(x0_.ptr(), ST::zero());
294
295 this->nominalValues_.set_x(x0_);
296 }
297
298
299 template<class SC,class LO,class GO,class NO>
300 void NonLinearProblem<SC,LO,GO,NO>::evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<SC> &inArgs,
301 const Thyra::ModelEvaluatorBase::OutArgs<SC> &outArgs
302 ) const
303 {
304 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
305 if ( !type.compare("Monolithic"))
306 evalModelImplMonolithic( inArgs, outArgs );
307 else if ( !type.compare("Teko")|| !type.compare("Diagonal") || type == "Triangular"|| !type.compare("PCD") || !type.compare("LSC")){
308 #ifdef FEDD_HAVE_TEKO
309 evalModelImplBlock( inArgs, outArgs );
310 #else
311 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Teko not found! Build Trilinos with Teko.");
312 #endif
313 }
314 else
315 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown preconditioner/solver type.");
316 }
317
322 template<class SC,class LO,class GO,class NO>
323 void NonLinearProblem<SC,LO,GO,NO>::evalModelImplMonolithic(const Thyra::ModelEvaluatorBase::InArgs<SC> &inArgs,
324 const Thyra::ModelEvaluatorBase::OutArgs<SC> &outArgs ) const
325 {
326 using Teuchos::RCP;
327 // using Teuchos::rcp;
328 using Teuchos::rcp_dynamic_cast;
329 using Teuchos::rcp_const_cast;
330 using Teuchos::ArrayView;
331 using Teuchos::Array;
332 // For output of xpetra/tpetra vectors
333 // RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
334 // RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
335
336 // Testing input arguments
337 TEUCHOS_TEST_FOR_EXCEPTION( inArgs.get_x().is_null(), std::logic_error, "inArgs.get_x() is null.");
338
339 RCP< const Thyra::VectorBase< SC > > vecThyra = inArgs.get_x();
340 RCP< Thyra::VectorBase< SC > > vecThyraNonConst = rcp_const_cast<Thyra::VectorBase< SC > >(vecThyra);
341
342 this->solution_->fromThyraMultiVector(vecThyraNonConst); // Setting solution to be the input vector inArgs
343
344 // Output arguments f and preconditioner
345 const RCP<Thyra::MultiVectorBase<SC> > f_out = outArgs.get_f();
346 const RCP<Thyra::LinearOpBase<SC> > W_out = outArgs.get_W_op();
347 const RCP<Thyra::PreconditionerBase<SC> > W_prec_out = outArgs.get_W_prec();
348
349 // Typedefs for Tpetra objects
350 typedef Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO> tpetra_extract;
351 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraMatrix_Type;
352 typedef RCP<TpetraMatrix_Type> TpetraMatrixPtr_Type;
353 typedef RCP<const TpetraMatrix_Type> TpetraMatrixConstPtr_Type;
354
355 // Determine operations to be done
356 const bool fill_f = nonnull(f_out);
357 const bool fill_W = nonnull(W_out);
358 const bool fill_W_prec = nonnull(W_prec_out);
359
360
361 if ( fill_f || fill_W || fill_W_prec ) {
362
363 // ****************
364 // Get the underlying xpetra objects
365 // ****************
366 // 1) Calculate residual vector f
367 if (fill_f) {
368
369 this->calculateNonLinResidualVec("standard"); // Calculating residual Vector
370
371 // Changing the residualVector into a ThyraMultivector
372
373 Teuchos::RCP<Thyra::MultiVectorBase<SC> > f_thyra = this->getResidualVector()->getThyraMultiVector();
374 f_out->assign(*f_thyra);
375 }
376
377 // 2) Calculate Newton tangent W
378 TpetraMatrixPtr_Type W;
379 if (fill_W) {
380 this->assemble("Newton"); // ReAssembling matrices with updated u in this class
381 this->setBoundariesSystem(); // setting boundaries to the system
382
383 // Cast reassembles matrix to tpetra operator/matrix
384 Teuchos::RCP<TpetraOp_Type> W_tpetra = tpetra_extract::getTpetraOperator(W_out);
385 Teuchos::RCP<TpetraMatrix_Type> W_tpetraMat = Teuchos::rcp_dynamic_cast<TpetraMatrix_Type>(W_tpetra);
386
387 TpetraMatrixConstPtr_Type W_systemTpetra = this->getSystem()->getMergedMatrix()->getTpetraMatrix();
388 TpetraMatrixPtr_Type W_systemTpetraNonConst = rcp_const_cast<TpetraMatrix_Type>(W_systemTpetra);
389
390 //Tpetra::CrsMatrixWrap<SC,LO,GO,NO>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO>&>(*W_systemXpetraNonConst);
391 //Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>& xTpetraMat = dynamic_cast<Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>&>(*crsOp.getCrsMatrix());
392
393 Teuchos::RCP<TpetraMatrix_Type> tpetraMatTpetra = W_systemTpetraNonConst; //xTpetraMat.getTpetra_CrsMatrixNonConst();
394
395 W_tpetraMat->resumeFill();
396
397 for (auto i=0; i<tpetraMatTpetra->getMap()->getLocalNumElements(); i++) {
398 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type indices; //ArrayView< const LO > indices
399 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type values;
400 tpetraMatTpetra->getLocalRowView( i, indices, values);
401 W_tpetraMat->replaceLocalValues( i, indices, values);
402 }
403 W_tpetraMat->fillComplete();
404
405 }
406 // 3) Compute preconditioner W_prec if needed
407 if (fill_W_prec) {
408
409 if (precInitOnly_){
410 // We have the option to reuse the preconditioner after the first X Newtonsteps
411 int newtonLimit = this->parameterList_->sublist("Parameter").get("newtonLimit",2);
412 if(this->newtonStep_ < newtonLimit || this->parameterList_->sublist("Parameter").get("Rebuild Preconditioner every Newton Iteration",true) )
413 {
414 this->setupPreconditioner( "Monolithic" ); // Rebuilding preconditioner
415 }
416 else{ // Reusing preconditioner
417 if (this->verbose_)
418 std::cout << " NonLinearProblem<SC,LO,GO,NO>::evalModelImplMonolithic:: Skipping preconditioner reconstruction" << std::endl;
419 }
420 }
421 else
422 precInitOnly_ = true;
423
424 // ch 26.04.19: After each setup of the preconditioner we check if we use a two-level precondtioner with multiplicative combination between the levels.
425 // If this is the case, we need to pre apply the coarse level to the residual(f_out).
426
427 std::string levelCombination = this->parameterList_->sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive");
428 if (!levelCombination.compare("Multiplicative")) {
429 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Multiplicative Level Combination is not supported for NOX.");
430 }
431
432 this->newtonStep_ ++;
433
434 }
435 }
436 }
441 #ifdef FEDD_HAVE_TEKO
442 template<class SC,class LO,class GO,class NO>
443 void NonLinearProblem<SC,LO,GO,NO>::evalModelImplBlock(const Thyra::ModelEvaluatorBase::InArgs<SC> &inArgs,
444 const Thyra::ModelEvaluatorBase::OutArgs<SC> &outArgs ) const
445 {
446 using Teuchos::RCP;
447 // using Teuchos::rcp;
448 using Teuchos::rcp_dynamic_cast;
449 using Teuchos::rcp_const_cast;
450 using Teuchos::ArrayView;
451 using Teuchos::Array;
452
453 // For output of xpetra/tpetra vectors
454 // RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
455 // RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
456
457 // Testing input arguments
458 TEUCHOS_TEST_FOR_EXCEPTION( inArgs.get_x().is_null(), std::logic_error, "inArgs.get_x() is null.");
459
460 RCP< const Thyra::VectorBase< SC > > vecThyra = inArgs.get_x();
461
462 RCP< Thyra::VectorBase< SC > > vecThyraNonConst = rcp_const_cast<Thyra::VectorBase< SC > >(vecThyra);
463
464 RCP< Thyra::ProductVectorBase< SC > > vecThyraBlock = rcp_dynamic_cast<Thyra::ProductVectorBase< SC > > (vecThyraNonConst);
465
466 this->solution_->getBlockNonConst(0)->fromThyraMultiVector( vecThyraBlock->getNonconstVectorBlock(0) );
467 this->solution_->getBlockNonConst(1)->fromThyraMultiVector( vecThyraBlock->getNonconstVectorBlock(1) );
468
469 // Output arguments f and preconditioner
470 const RCP<Thyra::MultiVectorBase<SC> > f_out = outArgs.get_f();
471 const RCP<Thyra::LinearOpBase<SC> > W_out = outArgs.get_W_op();
472 const RCP<Thyra::PreconditionerBase<SC> > W_prec_out = outArgs.get_W_prec();
473
474 // Typedefs for Tpetra objects
475 typedef Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO> tpetra_extract;
476 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraMatrix_Type;
477 typedef RCP<TpetraMatrix_Type> TpetraMatrixPtr_Type;
478 typedef RCP<const TpetraMatrix_Type> TpetraMatrixConstPtr_Type;
479
480 // Determine operations to be done
481 const bool fill_f = nonnull(f_out);
482 const bool fill_W = nonnull(W_out);
483 const bool fill_W_prec = nonnull(W_prec_out);
484
485 if ( fill_f || fill_W || fill_W_prec ) {
486
487 // ****************
488 // Get the underlying tpetra objects
489 // ****************
490 // 1) Calculate residual vector f
491 if (fill_f) {
492
493 this->calculateNonLinResidualVec("standard"); // Calculating residual Vector
494
495 Teko::MultiVector f0;
496 Teko::MultiVector f1;
497 f0 = this->getResidualVector()->getBlockNonConst(0)->getThyraMultiVector();
498 f1 = this->getResidualVector()->getBlockNonConst(1)->getThyraMultiVector();
499
500 std::vector<Teko::MultiVector> f_vec; f_vec.push_back(f0); f_vec.push_back(f1);
501
502 Teko::MultiVector f = Teko::buildBlockedMultiVector(f_vec);
503
504 f_out->assign(*f);
505 }
506
507 // 2) Calculate Newton tangent W
508 TpetraMatrixPtr_Type W;
509 if (fill_W) {
510
511 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMatrix;
512
513 this->assemble("Newton"); // ReAssembling matrices with updated u in this class
514 this->setBoundariesSystem(); // setting boundaries to the system
515
516 // Cast reassembles matrix to tpetra operator/matrix
517 RCP<ThyraBlockOp_Type> W_blocks = rcp_dynamic_cast<ThyraBlockOp_Type>(W_out);
518 RCP<const ThyraOp_Type> W_block00 = W_blocks->getBlock(0,0);
519 RCP<ThyraOp_Type> W_block00NonConst = rcp_const_cast<ThyraOp_Type>( W_block00 );
520 RCP<TpetraOp_Type> W_tpetra = tpetra_extract::getTpetraOperator( W_block00NonConst );
521
522 RCP<TpetraMatrix_Type> W_tpetraMat = Teuchos::rcp_dynamic_cast<TpetraMatrix_Type>(W_tpetra);
523
524 TpetraMatrixConstPtr_Type W_matrixTpetra = this->getSystem()->getBlock(0,0)->getTpetraMatrix();
525 TpetraMatrixPtr_Type W_matrixTpetraNonConst = rcp_const_cast<TpetraMatrix_Type>(W_matrixTpetra);
526 RCP<TpetraMatrix_Type> tpetraMatTpetra = W_matrixTpetraNonConst;
527
528 W_tpetraMat->resumeFill();
529
530 for (auto i=0; i<tpetraMatTpetra->getMap()->getLocalNumElements(); i++) {
531 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type indices; //ArrayView< const LO > indices
532 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type values;
533 tpetraMatTpetra->getLocalRowView( i, indices, values);
534 W_tpetraMat->replaceLocalValues( i, indices, values);
535 }
536 W_tpetraMat->fillComplete();
537
538 }
539 // 3) Compute preconditioner W_prec if needed
540 if (fill_W_prec) {
541 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
542 if (precInitOnly_){
543 // We have the option to reuse the preconditioner after the first X Newtonsteps
544 int newtonLimit = this->parameterList_->sublist("Parameter").get("newtonLimit",2);
545 if(this->newtonStep_ < newtonLimit || this->parameterList_->sublist("Parameter").get("Rebuild Preconditioner every Newton Iteration",true) )
546 {
547 this->setupPreconditioner( type );
548 }
549 else{
550 if (this->verbose_)
551 std::cout << " \n NonLinearProblem<SC,LO,GO,NO>::evalModelImplBlock:: Skipping preconditioner reconstruction \n " << std::endl;
552 }
553 }
554 else
555 precInitOnly_ = true;
556 // ch 26.04.19: After each setup of the preconditioner we check if we use a two-level precondtioner with multiplicative combination between the levels.
557 // If this is the case, we need to pre apply the coarse level to the residual(f_out).
558
559 ParameterListPtr_Type tmpSubList = sublist( sublist( sublist( sublist( this->parameterList_, "Teko Parameters" ) , "Preconditioner Types" ) , "Teko" ) , "Inverse Factory Library" );
560
561 std::string levelCombination1 = tmpSubList->sublist( "FROSch-Velocity" ).get("Level Combination","Additive");
562 std::string levelCombination2 = tmpSubList->sublist( "FROSch-Pressure" ).get("Level Combination","Additive");
563
564 if ( !levelCombination1.compare("Multiplicative") || !levelCombination2.compare("Multiplicative") ) {
565
566 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Multiplicative Level Combination is not supported for NOX.");
567 ParameterListPtr_Type solverPList = this->getLinearSolverBuilder()->getNonconstParameterList();
568 }
569 this->newtonStep_ ++;
570
571 }
572 }
573 }
574 #endif
575
576
577 template<class SC,class LO,class GO,class NO>
578 Teuchos::RCP<Thyra::LinearOpBase<SC> > NonLinearProblem<SC,LO,GO,NO>::create_W_op() const
579 {
580 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
581 if ( !type.compare("Monolithic"))
582 return create_W_op_Monolithic( );
583 else if ( !type.compare("Teko") || !type.compare("Diagonal") || !type.compare("Triangular") || !type.compare("PCD") || !type.compare("LSC") ){
584 #ifdef FEDD_HAVE_TEKO
585 return create_W_op_Block( );
586 #else
587 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Teko not found! Build Trilinos with Teko.");
588 #endif
589 }
590 else
591 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown preconditioner/solver type.");
592 }
593
594 template<class SC,class LO,class GO,class NO>
595 Teuchos::RCP<Thyra::LinearOpBase<SC> > NonLinearProblem<SC,LO,GO,NO>::create_W_op_Monolithic() const
596 {
597 Teuchos::RCP<const Thyra::LinearOpBase<SC> > W_opConst = this->system_->getThyraLinOp();
598 Teuchos::RCP<Thyra::LinearOpBase<SC> > W_op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<SC> >(W_opConst);
599 return W_op;
600 }
601
602 #ifdef FEDD_HAVE_TEKO
603 template<class SC,class LO,class GO,class NO>
604 Teuchos::RCP<Thyra::LinearOpBase<SC> > NonLinearProblem<SC,LO,GO,NO>::create_W_op_Block() const
605 {
606
607 Teko::LinearOp thyraF = this->system_->getBlock(0,0)->getThyraLinOp();
608 Teko::LinearOp thyraBT = this->system_->getBlock(0,1)->getThyraLinOp();
609 Teko::LinearOp thyraB = this->system_->getBlock(1,0)->getThyraLinOp();
610
611 if (!this->system_->blockExists(1,1)){
612 MatrixPtr_Type dummy = Teuchos::rcp( new Matrix_Type( this->system_->getBlock(1,0)->getMap(), 1 ) );
613 dummy->fillComplete();
614 this->system_->addBlock( dummy, 1, 1 );
615 }
616
617 Teko::LinearOp thyraC = this->system_->getBlock(1,1)->getThyraLinOp();
618
619 Teuchos::RCP<const Thyra::LinearOpBase<SC> > W_opConst = Thyra::block2x2(thyraF,thyraBT,thyraB,thyraC);
620 Teuchos::RCP<Thyra::LinearOpBase<SC> > W_op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<SC> >(W_opConst);
621 return W_op;
622 }
623 #endif
624
625 template<class SC,class LO,class GO,class NO>
626 Teuchos::RCP<Thyra::PreconditionerBase<SC> > NonLinearProblem<SC,LO,GO,NO>::create_W_prec() const
627 {
628 this->initializeSolverBuilder();
629
630 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
631 this->setBoundariesSystem();
632
633 if (!type.compare("Teko") || !type.compare("Diagonal") || !type.compare("Triangular") || !type.compare("PCD") || !type.compare("LSC")) { //
634 this->setupPreconditioner( type );
635 precInitOnly_ = false;
636 }
637 else{
638 this->setupPreconditioner( type ); // initializePreconditioner( type );
639 precInitOnly_ = false;
640 }
641
642
643 Teuchos::RCP<const Thyra::PreconditionerBase<SC> > thyraPrec = this->getPreconditionerConst()->getThyraPrecConst();
644 Teuchos::RCP<Thyra::PreconditionerBase<SC> > thyraPrecNonConst = Teuchos::rcp_const_cast<Thyra::PreconditionerBase<SC> >(thyraPrec);
645
646 return thyraPrecNonConst;
647
648 }
649
650}
651
652#endif
FEDD::NonLinearProblem
Definition NonLinearProblem_decl.hpp:28
FEDD::NonLinearProblem::initializeVectorsNonLinear
void initializeVectorsNonLinear(int nmbVectors=1)
Initialisation of the non-linear vectors like residual and previous solution.
Definition NonLinearProblem_def.hpp:60
FEDD::NonLinearProblem::getResidualVector
BlockMultiVectorPtr_Type getResidualVector() const
Get the residual vector.
Definition NonLinearProblem_def.hpp:172
FEDD::NonLinearProblem::initializeProblem
void initializeProblem(int nmbVectors=1)
Initialisation of the non-linear problem with system, vectors, and Thyra vector spaces for NOX.
Definition NonLinearProblem_def.hpp:45
FEDD::NonLinearProblem::solveUpdate
int solveUpdate()
This is where the linear solve specifically happens.
Definition NonLinearProblem_def.hpp:143
FEDD::NonLinearProblem::create_W_op
Teuchos::RCP< Thyra::LinearOpBase< SC > > create_W_op() const
Block Approach for Nonlinear Solver NOX. Input. Includes calculation of the residual vector and updat...
Definition NonLinearProblem_def.hpp:578
FEDD::NonLinearProblem::calculateResidualNorm
double calculateResidualNorm() const
Calculate the 2-norm of the residual vector.
Definition NonLinearProblem_def.hpp:133
FEDD::NonLinearProblem::NonLinearProblem
NonLinearProblem(CommConstPtr_Type comm)
Constructor.
Definition NonLinearProblem_def.hpp:23
FEDD::NonLinearProblem::solveAndUpdate
int solveAndUpdate(const std::string &criterion, double &criterionValue)
Solving the non-linear problem and updating the solution.
Definition NonLinearProblem_def.hpp:155
FEDD::NonLinearProblem::calculateNonLinResidualVec
virtual void calculateNonLinResidualVec(std::string type="standard", double time=0.) const =0
Virtual function which is implemented in the specific non-linear problem classes to calculate the non...
FEDD::NonLinearProblem::infoNonlinProblem
void infoNonlinProblem()
Information about the non-linear problem.
Definition NonLinearProblem_def.hpp:91
FEDD::SmallMatrix
This class represents a templated small Matrix of type T. Primarily created for 2x2 and 3x3 matrices....
Definition SmallMatrix.hpp:20
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36