Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
Problem_def.hpp
1#ifndef PROBLEM_DEF_hpp
2#define PROBLEM_DEF_hpp
3
4#include "feddlib/core/LinearAlgebra/BlockMatrix.hpp"
5#include "feddlib/core/FE/Domain.hpp"
6#include "feddlib/core/FE/FE.hpp"
7#include "feddlib/problems/Solver/Preconditioner.hpp"
8#include "feddlib/problems/Solver/LinearSolver.hpp"
9#include "feddlib/core/General/BCBuilder.hpp"
10#include "feddlib/core/LinearAlgebra/BlockMultiVector.hpp"
11
19
20
21namespace FEDD
22{
23 template <class SC, class LO, class GO, class NO>
24 Problem<SC, LO, GO, NO>::Problem(CommConstPtr_Type comm) : dim_(-1),
25 comm_(comm),
26 system_(),
27 rhs_(),
28 solution_(),
29 preconditioner_(),
30 linearSolverBuilder_(),
31 verbose_(comm->getRank() == 0),
32 parameterList_(),
33 domainPtr_vec_(),
34 domain_FEType_vec_(),
35 variableName_vec_(),
36 bcFactory_(),
37 feFactory_(),
38 dofsPerNode_vec_(),
39 sourceTerm_(),
40 rhsFuncVec_(),
41 parasSourceFunc_(0)
42#ifdef FEDD_TIMER
43 ,
44 solveProblemTimer_(Teuchos::TimeMonitor::getNewCounter("FEDD - Problem - Solve")), bcMatrixTimer_(Teuchos::TimeMonitor::getNewCounter("FEDD - Problem - Set Boundaries Matrix")), bcRHSTimer_(Teuchos::TimeMonitor::getNewCounter("FEDD - Problem - Set Boundaries RHS"))
45#endif
46 {
47 linearSolverBuilder_.reset(new Stratimikos::DefaultLinearSolverBuilder());
48 preconditioner_ = Teuchos::rcp(new Preconditioner_Type(this));
49 feFactory_.reset(new FEFac_Type());
50 }
51
52 template <class SC, class LO, class GO, class NO>
53 Problem<SC, LO, GO, NO>::Problem(ParameterListPtr_Type &parameterList, CommConstPtr_Type comm) : dim_(-1),
54 comm_(comm),
55 system_(),
56 rhs_(),
57 solution_(),
58 preconditioner_(),
59 linearSolverBuilder_(),
60 verbose_(comm->getRank() == 0),
61 parameterList_(parameterList),
62 domainPtr_vec_(),
63 domain_FEType_vec_(),
64 variableName_vec_(),
65 bcFactory_(),
66 feFactory_(),
67 dofsPerNode_vec_(),
68 sourceTerm_(),
69 rhsFuncVec_(),
70 parasSourceFunc_(0)
71#ifdef FEDD_TIMER
72 ,
73 solveProblemTimer_(Teuchos::TimeMonitor::getNewCounter("FEDD - Problem - Solve")), bcMatrixTimer_(Teuchos::TimeMonitor::getNewCounter("FEDD - Problem - Set Boundaries Matrix")), bcRHSTimer_(Teuchos::TimeMonitor::getNewCounter("FEDD - Problem - Set Boundaries RHS"))
74#endif
75 {
76 linearSolverBuilder_.reset(new Stratimikos::DefaultLinearSolverBuilder());
77 preconditioner_ = Teuchos::rcp(new Preconditioner_Type(this));
78 feFactory_.reset(new FEFac_Type());
79 }
80
81 template <class SC, class LO, class GO, class NO>
82 void Problem<SC, LO, GO, NO>::infoProblem()
83 {
84 bool verbose(comm_->getRank() == 0);
85 if (verbose)
86 {
87 std::cout << "\t ### Problem Information ###" << std::endl;
88 std::cout << "\t ### Dimension: " << dim_ << std::endl;
89 std::cout << "\t ### Linear problem: "
90 << "??" << std::endl;
91 std::cout << "\t ### Number of blocks/equations/variables: " << domainPtr_vec_.size() << std::endl;
92 for (int i = 0; i < domainPtr_vec_.size(); i++)
93 {
94 std::cout << "\t \t # Block " << i + 1 << "\t name: " << variableName_vec_.at(i) << "\t d.o.f.s: " << dofsPerNode_vec_.at(i) << "\t FE type: " << domain_FEType_vec_.at(i) << std::endl;
95 }
96
97 // ch 15.04.19: Hier ggf. unterscheiden zwischen Monolithic und Teko bzw. anderen Block-Precs.
98 ParameterListPtr_Type pListThyraPrec = sublist(parameterList_, "ThyraPreconditioner");
99 std::cout << "\t ### ### ###" << std::endl;
100 std::cout << "\t ### Preconditioner Information ###" << std::endl;
101 std::cout << "\t ### Type: " << parameterList_->sublist("General").get("Preconditioner Method", "Monolithic") << std::endl;
102 std::cout << "\t ### Prec.: " << pListThyraPrec->get("Preconditioner Type", "FROSch") << std::endl;
103
104 if (!pListThyraPrec->get("Preconditioner Type", "FROSch").compare("FROSch") && parameterList_->sublist("General").get("Preconditioner Method", "Monolithic") == "Monolithic")
105 {
106 std::cout << "\t ### Variant: " << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("FROSch Preconditioner Type", "TwoLevelBlockPreconditioner") << std::endl;
107 std::cout << "\t ### Two Level: "
108 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("TwoLevel", false)
109 << "\t Overlap: "
110 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("Overlap", 0)
111 << "\t Level Combination: "
112 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("Level Combination", "Additive") << std::endl;
113
114 std::cout << "\t OverlappingOperator Type: "
115 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("OverlappingOperator Type", "AlgebraicOverlappingOperator") << std::endl;
116
117 std::cout << "\t CoarseOperator Type: "
118 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("CoarseOperator Type", "GDSWCoarseOperator") << std::endl;
119
120 for (int i = 0; i < this->parameterList_->get("Number of blocks", 2); i++)
121 {
122 std::cout << "\t \t # Block " << i + 1 << "\t d.o.f.s: "
123 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("DofsPerNode" + std::to_string(i + 1), 1)
124 << "\t d.o.f. ordering: "
125 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("DofOrdering" + std::to_string(i + 1), "NodeWise") << std::endl;
126 }
127 }
128 else
129 {
130 std::cout << "\t ### Full preconditioner information only available for Monolithic preconditioner type ###" << std::endl;
131 }
132 }
133 }
134
135 template <class SC, class LO, class GO, class NO>
136 void Problem<SC, LO, GO, NO>::initializeProblem(int nmbVectors)
137 {
138
139 this->system_.reset(new BlockMatrix_Type(1));
140 this->initializeVectors(nmbVectors);
141 }
142
143 template <class SC, class LO, class GO, class NO>
144 void Problem<SC, LO, GO, NO>::addRhsFunction(RhsFunc_Type func)
145 {
146 this->rhsFuncVec_.push_back(func);
147 }
148 template <class SC, class LO, class GO, class NO>
149 void Problem<SC, LO, GO, NO>::addRhsFunction(RhsFunc_Type func, int i)
150 {
151 if (this->rhsFuncVec_.size() <= i + 1)
152 this->rhsFuncVec_[i] = func;
153 else if (this->rhsFuncVec_.size() == i)
154 this->rhsFuncVec_.push_back(func);
155 else
156 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Insertion Index smaller then size of vector");
157 }
158
159 /*template<class SC,class LO,class GO,class NO>
160 void Problem<SC,LO,GO,NO>::addRhsFunctionAndFlag(RhsFunc_Type func, int i, int flag){
161
162 TEUCHOS_TEST_FOR_EXCEPTION(!parameterList_->sublist("Parameter").get("Source Type","volume").compare("volume") , std::runtime_error, "Source Type Error: RHS Flags only make sence with surface source type.");
163
164 if(this->rhsFuncVec_.size() <= i+1){
165 this->rhsFuncVec_[i] = func;
166 this->rhsFuncFlagVec_[i] = flag;
167 }
168 else if(this->rhsFuncVec_.size() == i){
169 this->rhsFuncVec_.push_back(func);
170 this->rhsFuncFlagVec_.push_back(flag);
171 }
172 else
173 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Insertion Index smaller then size of vector");
174
175
176 }*/
177
178 template <class SC, class LO, class GO, class NO>
179 RhsFunc_Type &Problem<SC, LO, GO, NO>::getRhsFunction(int i)
180 {
181 return rhsFuncVec_[i];
182 }
183
184 template <class SC, class LO, class GO, class NO>
185 void Problem<SC, LO, GO, NO>::addVariable(const DomainConstPtr_Type &domain, std::string FEType, std::string name, int dofsPerNode)
186 {
187
188 domainPtr_vec_.push_back(domain);
189 domain_FEType_vec_.push_back(FEType);
190 variableName_vec_.push_back(name);
191 feFactory_->addFE(domain);
192 dofsPerNode_vec_.push_back(dofsPerNode);
193 }
194
195 // template<class SC,class LO,class GO,class NO>
196 // void Problem<SC,LO,GO,NO>::reAssemble(){
197 // if (verbose_) {
198 // cout << "Nothing to reassemble for linear problem." << endl;
199 // }
200 // }
201
202 template <class SC, class LO, class GO, class NO>
203 void Problem<SC, LO, GO, NO>::assembleSourceTerm(double time) const
204 {
205
206 TEUCHOS_TEST_FOR_EXCEPTION(sourceTerm_.is_null(), std::runtime_error, "Initialize source term before you assemble it - sourceTerm pointer is null");
207 this->sourceTerm_->putScalar(0.);
208 std::string sourceType = parameterList_->sublist("Parameter").get("Source Type", "volume");
209
210 if (sourceType == "volume")
211 assembleVolumeTerm(time);
212 else if (sourceType == "surface")
213 assembleSurfaceTerm(time);
214 }
215
216 template <class SC, class LO, class GO, class NO>
217 void Problem<SC, LO, GO, NO>::assembleVolumeTerm(double time) const
218 {
219 for (UN i = 0; i < sourceTerm_->size(); i++)
220 {
221 if (this->rhsFuncVec_.size() > i)
222 {
223 if (!this->rhsFuncVec_[i].empty())
224 {
225 MultiVectorPtr_Type FERhs;
226 // funcParameter[0] is always the time
227 vec_dbl_Type funcParameter(1, 0.);
228 funcParameter[0] = time;
229
230 // how can we use different parameters for different blocks here?
231 for (int j = 0; j < parasSourceFunc_.size(); j++)
232 funcParameter.push_back(parasSourceFunc_[j]);
233
234 std::string type;
235 if (this->getDofsPerNode(i) > 1)
236 {
237 FERhs = Teuchos::rcp(new MultiVector_Type(this->domainPtr_vec_.at(i)->getMapVecFieldRepeated()));
238 type = "Vector";
239 }
240 else
241 {
242 FERhs = Teuchos::rcp(new MultiVector_Type(this->domainPtr_vec_.at(i)->getMapRepeated()));
243 type = "Scalar";
244 }
245
246 this->feFactory_->assemblyRHS(this->dim_,
247 this->domain_FEType_vec_.at(i),
248 FERhs,
249 type,
250 this->rhsFuncVec_[i],
251 funcParameter);
252
253 this->sourceTerm_->getBlockNonConst(i)->exportFromVector(FERhs, false, "Add");
254 }
255 }
256 }
257 }
258
259 template <class SC, class LO, class GO, class NO>
260 void Problem<SC, LO, GO, NO>::assembleSurfaceTerm(double time) const
261 {
262 for (UN i = 0; i < sourceTerm_->size(); i++)
263 {
264 if (this->rhsFuncVec_.size() > i)
265 {
266
267 if (!this->rhsFuncVec_[i].empty())
268 {
269 MultiVectorPtr_Type FERhs;
270 // funcParameter[0] is always the time
271 vec_dbl_Type funcParameter(1, 0.);
272 funcParameter[0] = time;
273
274 // how can we use different parameters for different blocks here?
275 for (int j = 0; j < parasSourceFunc_.size(); j++)
276 funcParameter.push_back(parasSourceFunc_[j]);
277
278 // we add an additional parameter to place the surface flag of the element there during the assembly and shift the degree of the function to the last place now
279 funcParameter.push_back(funcParameter[funcParameter.size() - 1]);
280 std::string type;
281 if (this->getDofsPerNode(i) > 1)
282 {
283 FERhs = Teuchos::rcp(new MultiVector_Type(this->domainPtr_vec_.at(i)->getMapVecFieldRepeated()));
284 type = "Vector";
285 }
286 else
287 {
288 FERhs = Teuchos::rcp(new MultiVector_Type(this->domainPtr_vec_.at(i)->getMapRepeated()));
289 type = "Scalar";
290 }
291
292 this->feFactory_->assemblySurfaceIntegral(this->getDomain(i)->getDimension(),
293 this->getDomain(i)->getFEType(),
294 FERhs,
295 "Vector",
296 this->rhsFuncVec_[i],
297 funcParameter);
298
299 this->sourceTerm_->getBlockNonConst(i)->exportFromVector(FERhs, false, "Add");
300 }
301 }
302 }
303 // this->sourceTerm_->scale(-1.); // this scaling is needed for TPM problem
304 }
305
306 template <class SC, class LO, class GO, class NO>
307 int Problem<SC, LO, GO, NO>::solve(BlockMultiVectorPtr_Type rhs)
308 {
309 int its;
310 if (verbose_)
311 std::cout << "-- Solve System ..." << std::endl;
312 {
313
314#ifdef FEDD_TIMER
315 TimeMonitor_Type solveTM(*solveProblemTimer_);
316#endif
317 LinearSolver<SC, LO, GO, NO> linSolver;
318 std::string type = parameterList_->sublist("General").get("Preconditioner Method", "Monolithic");
319 its = linSolver.solve(this, rhs, type); // if rhs is null. Then the rhs_ of this is used in the linear solver
320 }
321
322 if (verbose_)
323 std::cout << " done -- " << std::endl;
324
325 return its;
326 }
327
328 template <class SC, class LO, class GO, class NO>
329 void Problem<SC, LO, GO, NO>::setupPreconditioner(std::string type) const
330 {
331
332 preconditioner_->buildPreconditioner(type);
333 }
334
335 template <class SC, class LO, class GO, class NO>
336 void Problem<SC, LO, GO, NO>::initializePreconditioner(std::string type) const
337 {
338
339 preconditioner_->initializePreconditioner(type);
340 }
341
342 template <class SC, class LO, class GO, class NO>
343 void Problem<SC, LO, GO, NO>::addBoundaries(const BCConstPtr_Type &bcFactory)
344 {
345
346 bcFactory_ = bcFactory;
347 }
348
349 template <class SC, class LO, class GO, class NO>
350 void Problem<SC, LO, GO, NO>::setBoundaries(double time) const
351 {
352#ifdef FEDD_TIMER
353 TimeMonitor_Type bcMatTM(*bcMatrixTimer_);
354 TimeMonitor_Type bcRHSTM(*bcRHSTimer_);
355#endif
356 bcFactory_->set(system_, rhs_, time);
357 }
358
359 template <class SC, class LO, class GO, class NO>
360 void Problem<SC, LO, GO, NO>::setBoundariesRHS(double time) const
361 {
362#ifdef FEDD_TIMER
363 TimeMonitor_Type bcRHSTM(*bcRHSTimer_);
364#endif
365 bcFactory_->setRHS(rhs_, time);
366 }
367
368 template <class SC, class LO, class GO, class NO>
369 void Problem<SC, LO, GO, NO>::setAllDirichletZero(BlockMultiVectorPtr_Type u) const
370 {
371#ifdef FEDD_TIMER
372 TimeMonitor_Type bcRHSTM(*bcRHSTimer_);
373#endif
374 bcFactory_->setAllDirichletZero(u);
375 }
376
377 template <class SC, class LO, class GO, class NO>
378 void Problem<SC, LO, GO, NO>::setBoundariesSystem() const
379 {
380#ifdef FEDD_TIMER
381 TimeMonitor_Type bcMatTM(*bcMatrixTimer_);
382#endif
383 bcFactory_->setSystem(system_);
384 }
385
386 template <class SC, class LO, class GO, class NO>
387 void Problem<SC, LO, GO, NO>::initializeVectors(int nmbVectors)
388 {
389
390 UN size = domainPtr_vec_.size();
391 solution_.reset(new BlockMultiVector_Type(size));
392 rhs_.reset(new BlockMultiVector_Type(size));
393 sourceTerm_.reset(new BlockMultiVector_Type(size));
394 rhsFuncVec_.resize(size);
395
396 for (UN i = 0; i < size; i++)
397 {
398 if (dofsPerNode_vec_[i] > 1)
399 {
400 MapConstPtr_Type map = domainPtr_vec_[i]->getMapVecFieldUnique();
401 MultiVectorPtr_Type solutionPart = Teuchos::rcp(new MultiVector_Type(map));
402 solution_->addBlock(solutionPart, i);
403 MultiVectorPtr_Type rhsPart = Teuchos::rcp(new MultiVector_Type(map));
404 rhs_->addBlock(rhsPart, i);
405 MultiVectorPtr_Type sourceTermPart = Teuchos::rcp(new MultiVector_Type(map));
406 sourceTerm_->addBlock(sourceTermPart, i);
407 }
408 else
409 {
410 MapConstPtr_Type map = domainPtr_vec_[i]->getMapUnique();
411 MultiVectorPtr_Type solutionPart = Teuchos::rcp(new MultiVector_Type(map));
412 solution_->addBlock(solutionPart, i);
413 MultiVectorPtr_Type rhsPart = Teuchos::rcp(new MultiVector_Type(map));
414 rhs_->addBlock(rhsPart, i);
415 MultiVectorPtr_Type sourceTermPart = Teuchos::rcp(new MultiVector_Type(map));
416 sourceTerm_->addBlock(sourceTermPart, i);
417 }
418 }
419 }
420
421 template <class SC, class LO, class GO, class NO>
422 typename Problem<SC, LO, GO, NO>::BlockMultiVectorPtr_Type Problem<SC, LO, GO, NO>::getRhs()
423 {
424
425 return rhs_;
426 }
427
428 template <class SC, class LO, class GO, class NO>
429 typename Problem<SC, LO, GO, NO>::BlockMultiVectorPtr_Type Problem<SC, LO, GO, NO>::getRhs() const
430 {
431
432 return rhs_;
433 }
434
435 template <class SC, class LO, class GO, class NO>
436 typename Problem<SC, LO, GO, NO>::BlockMultiVectorPtr_Type Problem<SC, LO, GO, NO>::getSolution()
437 {
438
439 return solution_;
440 }
441
442 template <class SC, class LO, class GO, class NO>
443 typename Problem<SC, LO, GO, NO>::BlockMatrixPtr_Type Problem<SC, LO, GO, NO>::getSystem() const
444 {
445
446 return system_;
447 }
448
449 template <class SC, class LO, class GO, class NO>
450 typename Problem<SC, LO, GO, NO>::PreconditionerPtr_Type Problem<SC, LO, GO, NO>::getPreconditioner()
451 {
452 return preconditioner_;
453 }
454
455 template <class SC, class LO, class GO, class NO>
456 typename Problem<SC, LO, GO, NO>::PreconditionerConstPtr_Type Problem<SC, LO, GO, NO>::getPreconditionerConst() const
457 {
458 return preconditioner_;
459 }
460
461 template <class SC, class LO, class GO, class NO>
462 void Problem<SC, LO, GO, NO>::setPreconditionerThyraFromLinOp(ThyraLinOpPtr_Type precLinOp)
463 {
464 preconditioner_->setPreconditionerThyraFromLinOp(precLinOp);
465 }
466
467 template <class SC, class LO, class GO, class NO>
468 bool Problem<SC, LO, GO, NO>::getVerbose() const
469 {
470
471 return verbose_;
472 }
473
474 template <class SC, class LO, class GO, class NO>
475 typename Problem<SC, LO, GO, NO>::FEFacConstPtr_Type Problem<SC, LO, GO, NO>::getFEFactory()
476 {
477
478 return feFactory_;
479 }
480
481 template <class SC, class LO, class GO, class NO>
482 typename Problem<SC, LO, GO, NO>::BCConstPtr_Type Problem<SC, LO, GO, NO>::getBCFactory()
483 {
484
485 return bcFactory_;
486 }
487
488 template <class SC, class LO, class GO, class NO>
489 typename Problem<SC, LO, GO, NO>::DomainConstPtr_Type Problem<SC, LO, GO, NO>::getDomain(int i) const
490 {
491
492 return domainPtr_vec_.at(i);
493 }
494
495 template <class SC, class LO, class GO, class NO>
496 std::string Problem<SC, LO, GO, NO>::getFEType(int i) const
497 {
498
499 return domain_FEType_vec_.at(i);
500 }
501
502 template <class SC, class LO, class GO, class NO>
503 std::string Problem<SC, LO, GO, NO>::getVariableName(int i) const
504 {
505
506 return variableName_vec_.at(i);
507 }
508
509 template <class SC, class LO, class GO, class NO>
510 int Problem<SC, LO, GO, NO>::getDofsPerNode(int i) const
511 {
512
513 return dofsPerNode_vec_.at(i);
514 }
515
516 template <class SC, class LO, class GO, class NO>
517 typename Problem<SC, LO, GO, NO>::ParameterListPtr_Type Problem<SC, LO, GO, NO>::getParameterList() const
518 {
519
520 return parameterList_;
521 }
522
523 template <class SC, class LO, class GO, class NO>
524 void Problem<SC, LO, GO, NO>::addToRhs(BlockMultiVectorPtr_Type x) const
525 {
526
527 rhs_->update(1., *x, 1.);
528 }
529
530 template <class SC, class LO, class GO, class NO>
531 typename Problem<SC, LO, GO, NO>::BlockMultiVectorPtr_Type Problem<SC, LO, GO, NO>::getSourceTerm()
532 {
533 TEUCHOS_TEST_FOR_EXCEPTION(sourceTerm_.is_null(), std::runtime_error, "Problem has no source term to return - source term pointer is null");
534 return sourceTerm_;
535 }
536
537 template <class SC, class LO, class GO, class NO>
538 bool Problem<SC, LO, GO, NO>::hasSourceTerm() const
539 {
540 for (int i = 0; i < rhsFuncVec_.size(); i++)
541 {
542 if (!rhsFuncVec_[i].empty())
543 return true;
544 }
545 return false;
546 }
547
548 template <class SC, class LO, class GO, class NO>
549 void Problem<SC, LO, GO, NO>::initSolutionWithVector(MultiVector_Type &mv)
550 {
551 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "initSolutionWithVector not implemented. DO we need this?");
552
553 // *solution_ = mv;
554 }
555
556 template <class SC, class LO, class GO, class NO>
557 void Problem<SC, LO, GO, NO>::initializeSolverBuilder() const
558 {
559
560 ParameterListPtr_Type pListThyraPrec = sublist(parameterList_, "ThyraPreconditioner");
561
562 linearSolverBuilder_->setParameterList(pListThyraPrec);
563 }
564
565
566
567 // Calls the FE Factory to change the linearization of all FE objects
568 template <class SC, class LO, class GO, class NO>
569 void Problem<SC, LO, GO, NO>::changeAssFELinearization(std::string linearization)
570 {
571 this->feFactory_->changeLinearizationFE(linearization);
572 }
573
574
575 // Functions that return the H1 and L2 Norm of a given Vector. The Norms being defined as:
576 // || mv ||_L2 = mv^T * M * mv
577 // || mv ||_H1 = mv^T * K * mv
578 // with M beeing the mass matrix and K beeing the stiffness matrix
579
580 template <class SC, class LO, class GO, class NO>
581 double Problem<SC, LO, GO, NO>::calculateL2Norm(MultiVectorConstPtr_Type mv, int domainInd)
582 {
583
584 MatrixPtr_Type M;
585 M = Teuchos::rcp(new Matrix_Type(this->domainPtr_vec_.at(domainInd)->getMapUnique(), this->getDomain(domainInd)->getApproxEntriesPerRow()));
586 this->feFactory_->assemblyMass(this->dim_, this->domain_FEType_vec_.at(domainInd), "Scalar", M);
587
588 Teuchos::RCP<MultiVector<SC, LO, GO, NO>> mvOutput = Teuchos::rcp(new MultiVector_Type(this->domainPtr_vec_.at(domainInd)->getMapUnique()));
589
590 M->apply(*mv, *mvOutput);
591
592 Teuchos::ArrayRCP<SC> vector = mv->getDataNonConst(0);
593 Teuchos::ArrayRCP<SC> outputVector = mvOutput->getDataNonConst(0);
594
595 double result = 0;
596 for (int i = 0; i < vector.size(); i++)
597 {
598 result += vector[i] * outputVector[i];
599 }
600 Teuchos::reduceAll<int, double>(*comm_, Teuchos::REDUCE_SUM, result, Teuchos::outArg(result));
601
602 return result;
603 }
604 template <class SC, class LO, class GO, class NO>
605 double Problem<SC, LO, GO, NO>::calculateH1Norm(MultiVectorConstPtr_Type mv, int blockId1, int blockId2, int domainInd)
606 {
607
608 MatrixPtr_Type K = this->getSystem()->getBlock(blockId1, blockId2);
609
610 Teuchos::RCP<MultiVector<SC, LO, GO, NO>> mvOutput = Teuchos::rcp(new MultiVector_Type(K->getMap()));
611
612 K->apply(*mv, *mvOutput); // this represents mvOutput = K * mv ;
613
614 Teuchos::ArrayRCP<SC> vector = mv->getDataNonConst(0);
615 Teuchos::ArrayRCP<SC> outputVector = mvOutput->getDataNonConst(0);
616
617 double result = 0;
618 for (int i = 0; i < vector.size(); i++)
619 {
620 result += vector[i] * outputVector[i];
621 }
622 Teuchos::reduceAll<int, double>(*comm_, Teuchos::REDUCE_SUM, result, Teuchos::outArg(result));
623
624 return result;
625 }
626
627 template <class SC, class LO, class GO, class NO>
628 void Problem<SC, LO, GO, NO>::infoParameter(bool full,std::string name)
629 {
630 bool verbose(comm_->getRank() == 0);
631 if (verbose)
632 {
633 std::cout << "#####################################" << std::endl;
634 if(name != "empty")
635 std::cout << " ### " << name << " Problem Information ###" << std::endl;
636 else
637 std::cout << " ### Problem Information ###" << std::endl;
638 std::cout << " ### Dimension: " << dim_ << std::endl;
639 std::cout << " ### Number of blocks/equations/variables: " << domainPtr_vec_.size() << std::endl;
640 for (int i = 0; i < domainPtr_vec_.size(); i++)
641 {
642 std::cout << " ## Block " << i + 1 << " name: " << variableName_vec_.at(i) << " d.o.f.s: " << dofsPerNode_vec_.at(i) << " FE type: " << domain_FEType_vec_.at(i) << std::endl;
643 }
644 if(full)
645 std::cout << " ####################################" << std::endl;
646 ParameterListPtr_Type parameterlist = sublist(parameterList_, "Parameter");
647 std::cout << " ### Parameter Information ###" << std::endl;
648 if(parameterlist->get("Mesh Type", "???") != "???"){
649 std::cout << " ### Mesh Type: " << parameterlist->get("Mesh Type", "???") << std::endl;
650 if(parameterlist->get("Mesh Type", "???") == "structured" || parameterlist->get("Mesh Type", "???") == "structured_bfs" )
651 std::cout << " ### H/h= " << parameterlist->get("H/h", 0) << std::endl;
652 else
653 std::cout << " ### Mesh File Name 1: " << parameterList_->sublist("Mesh Partitioner").get("Mesh 1 Name", "???") << std::endl;
654 }
655
656 // Geometry Problem
657 if(parameterlist->get("Model Geometry", "None") != "None" )
658 std::cout << " ### Model Geometry: " << parameterlist->get("Model Geometry", "None") << std::endl;
659 if(parameterlist->get("Model", "None") != "None" )
660 std::cout << " ### Model: " << parameterlist->get("Model", "None") << std::endl;
661 if(parameterlist->get("Distance Laplace", 0.) > 0 )
662 std::cout << " ### Distance Laplace: " << parameterlist->get("Distance Laplace", 0.) << std::endl;
663 if(parameterlist->get("Coefficient Laplace", 0.) > 0 )
664 std::cout << " ### Coefficient Laplace: " << parameterlist->get("Coefficient Laplace", 0.) << std::endl;
665 // Solid Problem
666 if(parameterlist->get("Material model", "None") != "None" )
667 std::cout << " ### Material model: " << parameterlist->get("Material model", "None") << std::endl;
668
669 if(parameterlist->get("Density", 0.) > 0 )
670 std::cout << " ### Density: " << parameterlist->get("Density", 0.) << std::endl;
671
672 // Fluid Problem
673 if(parameterlist->get("Viscosity", 0.) > 0)
674 std::cout << " ### Viscosity: " << parameterlist->get("Viscosity", 0.) << std::endl ;
675 if(abs(parameterlist->get("MaxVelocity", 0.)) > 0 )
676 std::cout << " ### Maximum Velocity: " << parameterlist->get("MaxVelocity", 0.) << std::endl;
677 else if(abs(parameterlist->get("Max Velocity", 0.)) > 0 )
678 std::cout << " ### or Maximum Velocity: " << parameterlist->get("Max Velocity", 0.) << std::endl;
679
680 if(abs(parameterlist->get("Flowrate", 0.)) > 0 )
681 std::cout << " ### Flowrate: " << parameterlist->get("Flowrate", 0.) << std::endl;
682
683 if((parameterlist->get("Pressure Boundary Condition", "None")) != "None" ){
684 std::cout << " ### Pressure Boundary Condition: " << parameterlist->get("Pressure Boundary Condition", "None") << std::endl;
685 std::cout << " ### Reference fluid pressure: " << parameterlist->get("Reference fluid pressure", 0.) << std::endl;
686 }
687 if(parameterlist->get("E", 0.) > 0 )
688 std::cout << " ### E Modul: " << parameterlist->get("E", 0.) << std::endl;
689 if(parameterlist->get("Mu", 0.) > 0 )
690 std::cout << " ### Mu: " << parameterlist->get("Mu", 0.) << std::endl;
691 if(parameterlist->get("Poisson Ratio", 0.) > 0 )
692 std::cout << " ### Poisson Ratio: " << parameterlist->get("Poisson Ratio", 0.) << std::endl;
693
694 if(full)
695 std::cout << " ### Rel. Tol.: " << parameterlist->get("relNonLinTol", 0.)
696 << " ### Abs. Tol.: " << parameterlist->get("absNonLinTol", 0.) << std::endl;
697
698
699 std::cout << "####################################" << std::endl;
700 // ch 15.04.19: Hier ggf. unterscheiden zwischen Monolithic und Teko bzw. anderen Block-Precs.
701 if(full){
702 ParameterListPtr_Type pListThyraPrec = sublist(parameterList_, "ThyraPreconditioner");
703 std::cout << " ### Preconditioner Information ###" << std::endl;
704 std::cout << " ## Type: " << parameterList_->sublist("General").get("Preconditioner Method", "Monolithic") << std::endl;
705
706 if (!pListThyraPrec->get("Preconditioner Type", "FROSch").compare("FROSch") && parameterList_->sublist("General").get("Preconditioner Method", "Monolithic") == "Monolithic")
707 {
708 std::cout << " ## Variant: " << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("FROSch Preconditioner Type", "TwoLevelBlockPreconditioner") << std::endl;
709 std::cout << " ## Overlap: "
710 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("Overlap", 0) << std::endl;
711 std::cout << " ## OverlappingOperator Type: "
712 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("OverlappingOperator Type", "AlgebraicOverlappingOperator") << std::endl;
713
714 std::cout << " ## CoarseOperator Type: "
715 << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("CoarseOperator Type", "GDSWCoarseOperator") << std::endl;
716
717 if(pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("CoarseOperator Type", "GDSWCoarseOperator") == "IPOUHarmonicCoarseOperator"){
718 for (int i = 0; i < this->parameterList_->get("Number of blocks", 2); i++)
719 {
720 std::cout << " # IPOU Block "<< std::to_string(i + 1) <<": " << pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").sublist("IPOUHarmonicCoarseOperator").sublist("Blocks").sublist(std::to_string(i + 1)).sublist("InterfacePartitionOfUnity").get("Type","NOTFOUND") << std::endl;
721 }
722 }
723
724 }
725 else if (!parameterList_->sublist("General").get("Preconditioner Method", "Monolithic").compare("Teko")){
726 std::cout << " ### Block Preconditioner Type: " << parameterList_->sublist("Teko Parameters").sublist("Preconditioner Types").sublist("Teko").get("Inverse Type","SIMPLE") << std::endl;
727 std::cout << " ### Velocity Preconditioner: " << parameterList_->sublist("Teko Parameters").sublist("Preconditioner Types").sublist("Teko").sublist("Inverse Factory Library").sublist("FROSch-Velocity").get("CoarseOperator Type","GDSW#") << std::endl;
728 std::cout << " ### Pressure Preconditioner: " << parameterList_->sublist("Teko Parameters").sublist("Preconditioner Types").sublist("Teko").sublist("Inverse Factory Library").sublist("FROSch-Pressure").get("CoarseOperator Type","GDSW#") << std::endl;
729
730 }
731 else if (!parameterList_->sublist("General").get("Preconditioner Method", "Monolithic").compare("Diagonal") ||
732 !parameterList_->sublist("General").get("Preconditioner Method", "Monolithic").compare("Triangular") ||
733 !parameterList_->sublist("General").get("Preconditioner Method", "Monolithic").compare("PCD") ||
734 !parameterList_->sublist("General").get("Preconditioner Method", "Monolithic").compare("LSC"))
735 {
736 // ParameterListPtr_Type pListThyraPrecBlock = sublist(parameterList_, "Block Preconditioner");
737
738 std::cout << " ### Block Preconditioner Type: " << parameterList_->sublist("General").get("Preconditioner Method", "Monolithic") << std::endl;
739 std::cout << " ## Velocity Preconditioner: " << parameterList_->sublist("Velocity preconditioner").sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("CoarseOperator Type", "Nada") << std::endl;
740 std::cout << " # Variant: " << parameterList_->sublist("Velocity preconditioner").sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("FROSch Preconditioner Type", "TwoLevelBlockPreconditioner") << std::endl;
741 std::cout << " # Overlap: "
742 << parameterList_->sublist("Velocity preconditioner").sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("Overlap", 0) << std::endl;
743
744 std::cout << " ## Schur Complement Preconditioner: " << parameterList_->sublist("Schur complement preconditioner").sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("CoarseOperator Type", "Nada") << std::endl;
745 std::cout << " # Variant: " << parameterList_->sublist("Schur complement preconditioner").sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("FROSch Preconditioner Type", "TwoLevelBlockPreconditioner") << std::endl;
746 std::cout << " # Overlap: "
747 << parameterList_->sublist("Schur complement preconditioner").sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("Overlap", 0) << std::endl;
748 }
749 else
750 {
751 std::cout << " ### Full preconditioner information only available for Monolithic/Teko preconditioner type ###" << std::endl;
752 }
753
754 std::cout << " ####################################" << std::endl;
755
756 std::cout << " ### Git commit hash: " << GIT_COMMIT_HASH << std::endl;
757 std::cout << " ### Git branch name: " << GIT_BRANCH_NAME << std::endl;
758
759 std::cout << "#####################################" << std::endl;
760 }
761 }
762 }
763
764
765}
766#endif
FEDD::Problem::addRhsFunction
void addRhsFunction(RhsFunc_Type func)
Definition Problem_def.hpp:144
FEDD::Problem::getRhsFunction
RhsFunc_Type & getRhsFunction(int i)
Definition Problem_def.hpp:179
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36