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