Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
Preconditioner_def.hpp
1#ifndef PRECONDITIONER_START
2#define PRECONDITIONER_START(A,S) Teuchos::RCP<Teuchos::TimeMonitor> A = Teuchos::rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(std::string("Preconditioner: ") + std::string(S))));
3#endif
4
5#ifndef PRECONDITIONER_STOP
6#define PRECONDITIONER_STOP(A) A.reset();
7#endif
8
9#ifndef Preconditioner_DEF_hpp
10#define Preconditioner_DEF_hpp
11#include <Thyra_DefaultZeroLinearOp_decl.hpp>
12#ifdef FEDD_HAVE_IFPACK2
13#include <Thyra_Ifpack2PreconditionerFactory_def.hpp>
14#endif
15
16#include "feddlib/core/FE/Domain.hpp"
17#include "feddlib/core/General/BCBuilder.hpp"
18
19#include "feddlib/problems/abstract/MinPrecProblem.hpp"
20#include "feddlib/problems/abstract/Problem.hpp"
21#include "feddlib/problems/abstract/TimeProblem.hpp"
22#include "feddlib/problems/specific/NavierStokes.hpp"
23#include "feddlib/problems/specific/NonLinElasticity.hpp"
24#include "feddlib/problems/specific/LinElas.hpp"
25#include "feddlib/problems/specific/Geometry.hpp"
26#include "feddlib/problems/specific/FSI.hpp"
27#include "feddlib/problems/Solver/PrecOpFaCSI.hpp"
28#include "feddlib/problems/Solver/PrecBlock2x2.hpp"
29
38
39namespace FEDD {
40template <class SC,class LO,class GO,class NO>
41Preconditioner<SC,LO,GO,NO>::Preconditioner(Problem_Type* problem):
42thyraPrec_(),
43precondtionerIsBuilt_(false),
44problem_(),
45timeProblem_(),
46precFactory_()
47#ifdef FEDD_HAVE_TEKO
48,tekoLinOp_()
49,velocityMassMatrix_()
50,rh_()
51#endif
52,fsiLinOp_()
53,precFluid_()
54,precStruct_()
55,precGeo_()
56,probFluid_()
57,probSolid_()
58,probGeo_()
59#ifdef PRECONDITIONER_TIMER
60,preconditionerTimer_ (Teuchos::TimeMonitor::getNewCounter("Preconditioner: Setup"))
61#endif
62{
63 problem_.reset( problem, false );
64
65 Stratimikos::enableFROSch<LO,GO,NO>( *problem_->getLinearSolverBuilder() );
66
67#ifdef FEDD_HAVE_IFPACK2
68 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> CRSMat;
69 typedef Tpetra::Operator<SC,LO,GO,NO> TP_op;
70 typedef Thyra::PreconditionerFactoryBase<SC> Base;
71 typedef Thyra::Ifpack2PreconditionerFactory<CRSMat> Impl;
72
73 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMatrix_Type;
74 typedef Teuchos::RCP<TpetraCrsMatrix_Type> TpetraCrsMatrixPtr_Type;
75
76 problem_->getLinearSolverBuilder()->setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
77#endif
78#ifdef FEDD_HAVE_TEKO
79 Teko::addTekoToStratimikosBuilder( *problem_->getLinearSolverBuilder() );
80#endif
81
82
83
84}
85
86template <class SC,class LO,class GO,class NO>
87Preconditioner<SC,LO,GO,NO>::Preconditioner(TimeProblem_Type* problem):
88thyraPrec_(),
89precondtionerIsBuilt_(false),
90problem_(),
91timeProblem_(),
92precFactory_()
93#ifdef FEDD_HAVE_TEKO
94,tekoLinOp_()
95,velocityMassMatrix_()
96,rh_()
97#endif
98,fsiLinOp_()
99,precFluid_()
100,precStruct_()
101,precGeo_()
102,probFluid_()
103,probSolid_()
104,probGeo_()
105#ifdef PRECONDITIONER_TIMER
106,preconditionerTimer_ (Teuchos::TimeMonitor::getNewCounter("Preconditioner: Setup"))
107#endif
108{
109 timeProblem_.reset( problem, false );
110
111 // We need to ensure that already set information is upheld
112 if(!problem->getUnderlyingProblem()->preconditioner_->getPressureProjection().is_null()){
113 setPressureProjection(problem->getUnderlyingProblem()->preconditioner_->getPressureProjection());
114 }
115
116 if(!problem->getUnderlyingProblem()->preconditioner_->getVelocityMassMatrix().is_null()){
117 setVelocityMassMatrix(problem->getUnderlyingProblem()->preconditioner_->getVelocityMassMatrix());
118 }
119
120 if(!problem->getUnderlyingProblem()->preconditioner_->getPressureLaplaceMatrix().is_null()){
121 setPressureLaplaceMatrix(problem->getUnderlyingProblem()->preconditioner_->getPressureLaplaceMatrix());
122 }
123 if(!problem->getUnderlyingProblem()->preconditioner_->getPressureMassMatrix().is_null()){
124 setPressureMassMatrix(problem->getUnderlyingProblem()->preconditioner_->getPressureMassMatrix());
125
126 }
127 if(!problem->getUnderlyingProblem()->preconditioner_->getPCDOperatorMatrix().is_null()){
128 setPCDOperator(problem->getUnderlyingProblem()->preconditioner_->getPCDOperatorMatrix());
129 }
130
131}
132
133template <class SC,class LO,class GO,class NO>
134Preconditioner<SC,LO,GO,NO>::~Preconditioner()
135{
136}
137
138template<class SC,class LO,class GO,class NO>
139void Preconditioner<SC,LO,GO,NO>::setPreconditionerThyraFromLinOp( ThyraLinOpPtr_Type precLinOp ){
140 TEUCHOS_TEST_FOR_EXCEPTION( thyraPrec_.is_null(), std::runtime_error, "thyraPrec_ is null." );
141 Teuchos::RCP<Thyra::DefaultPreconditioner<SC> > defaultPrec = Teuchos::rcp_dynamic_cast<Thyra::DefaultPreconditioner<SC> > (thyraPrec_);
142 TEUCHOS_TEST_FOR_EXCEPTION( defaultPrec.is_null(), std::runtime_error, "Cast to DefaultPrecondtioner failed." );
143 TEUCHOS_TEST_FOR_EXCEPTION( precLinOp.is_null(), std::runtime_error, "precLinOp is null." );
144 defaultPrec->initializeUnspecified( precLinOp );
145}
146
147template <class SC,class LO,class GO,class NO>
148typename Preconditioner<SC,LO,GO,NO>::ThyraPrecPtr_Type Preconditioner<SC,LO,GO,NO>::getThyraPrec(){
149 return thyraPrec_;
150}
151
152template <class SC,class LO,class GO,class NO>
153typename Preconditioner<SC,LO,GO,NO>::ThyraPrecConstPtr_Type Preconditioner<SC,LO,GO,NO>::getThyraPrecConst() const{
154 return thyraPrec_;
155}
156
157template <class SC,class LO,class GO,class NO>
158void Preconditioner<SC,LO,GO,NO>::initializePreconditioner( std::string type )
159{
160 if ( type == "Monolithic" || type == "FaCSI" || type == "Diagonal" || type == "Triangular"){
161 if (type == "Monolithic")
162 initPreconditionerMonolithic( );
163 else if (type == "FaCSI" || type == "Diagonal" || type == "Triangular" || type == "PCD" || type == "LSC")
164 initPreconditionerBlock( );
165
166 }
167 else if (type == "Teko" || type == "FaCSI-Teko" || type == "FaCSI-Block"){
168 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Please construct the Teko precondtioner completely.");
169 }
170 else
171 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown preconditioner type for initialization.");
172}
173
174template <class SC,class LO,class GO,class NO>
175void Preconditioner<SC,LO,GO,NO>::setPressureProjection(BlockMultiVectorPtr_Type pressureProjection) const{
176 pressureProjection_ = pressureProjection;
177}
178
179template <class SC,class LO,class GO,class NO>
180void Preconditioner<SC,LO,GO,NO>::initPreconditionerMonolithic( )
181{
182
183 LinSolverBuilderPtr_Type solverBuilder;
184 Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > thyraRangeSpace;
185 Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > thyraDomainSpace;
186
187 if (!problem_.is_null()){
188 solverBuilder = problem_->getLinearSolverBuilder();
189 thyraRangeSpace = Thyra::tpetraVectorSpace<SC,LO,GO,NO>( problem_->getSystem()->getMap()->getMergedMap()->getTpetraMap()); //Tpetra::ThyraUtils<SC,LO,GO,NO>::toThyra( problem_->getSystem()->getMap()->getMergedMap()->getTpetraMap() );
190 thyraDomainSpace = Thyra::tpetraVectorSpace<SC,LO,GO,NO>( problem_->getSystem()->getMap()->getMergedMap()->getTpetraMap()); //Tpetra::ThyraUtils<SC,LO,GO,NO>::toThyra( problem_->getSystem()->getMap()->getMergedMap()->getTpetraMap() );
191 }
192 else if(!timeProblem_.is_null()){
193 solverBuilder = timeProblem_->getUnderlyingProblem()->getLinearSolverBuilder();
194 thyraRangeSpace = Thyra::tpetraVectorSpace<SC,LO,GO,NO>( timeProblem_->getSystem()->getMap()->getMergedMap()->getTpetraMap() );
195 thyraDomainSpace = Thyra::tpetraVectorSpace<SC,LO,GO,NO>( timeProblem_->getSystem()->getMap()->getMergedMap()->getTpetraMap() );
196
197 }
198
199 Teuchos::RCP<Thyra::LinearOpBase<SC> > thyraLinOp = Teuchos::rcp( new Thyra::DefaultZeroLinearOp<SC>( thyraRangeSpace, thyraDomainSpace ) );
200
201 Teuchos::RCP<Thyra::PreconditionerFactoryBase<SC> > precFactory = solverBuilder->createPreconditioningStrategy("");
202 thyraPrec_ = precFactory->createPrec();
203
204 Teuchos::RCP<Thyra::DefaultPreconditioner<SC> > defaultPrec = Teuchos::rcp_dynamic_cast<Thyra::DefaultPreconditioner<SC> >(thyraPrec_);
205
206 defaultPrec->initializeUnspecified(thyraLinOp);
207}
208
209template <class SC,class LO,class GO,class NO>
210void Preconditioner<SC,LO,GO,NO>::initPreconditionerBlock( )
211{
212 using Teuchos::Array;
213 using Teuchos::RCP;
214 using Teuchos::rcp;
215 BlockMultiVectorPtr_Type sol;
216 BlockMultiVectorPtr_Type rhs;
217 LinSolverBuilderPtr_Type solverBuilder;
218 int size = 0;
219 if (!problem_.is_null()){
220 solverBuilder = problem_->getLinearSolverBuilder();
221 size = problem_->getSystem()->size();
222 sol = problem_->getSolution();
223 rhs = problem_->getRhs();
224 }
225 else if(!timeProblem_.is_null()){
226 solverBuilder = timeProblem_->getUnderlyingProblem()->getLinearSolverBuilder();
227 size = timeProblem_->getSystem()->size();
228 sol = timeProblem_->getSolution();
229 rhs = timeProblem_->getRhs();
230 }
231
232 Array< RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesRange( size );
233 Array< RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesDomain( size );
234
235 for (int i=0; i<size; i++) {
236 vectorSpacesRange[i] = rhs->getBlock(i)->getMap()->getThyraVectorSpaceBase();
237 vectorSpacesDomain[i] = sol->getBlock(i)->getMap()->getThyraVectorSpaceBase();
238 }
239
240 RCP<const Thyra::DefaultProductVectorSpace<SC> > pR = Thyra::productVectorSpace<SC>( vectorSpacesRange );
241 RCP<const Thyra::DefaultProductVectorSpace<SC> > pD = Thyra::productVectorSpace<SC>( vectorSpacesDomain );
242
243 RCP<Thyra::LinearOpBase<SC> > thyraLinOp = rcp( new Thyra::DefaultZeroLinearOp<SC>( pR, pD ) );
244
245 Teuchos::RCP<Thyra::PreconditionerFactoryBase<SC> > precFactory = solverBuilder->createPreconditioningStrategy("");
246 thyraPrec_ = precFactory->createPrec();
247
248 Teuchos::RCP<Thyra::DefaultPreconditioner<SC> > defaultPrec = Teuchos::rcp_dynamic_cast<Thyra::DefaultPreconditioner<SC> >(thyraPrec_);
249
250 defaultPrec->initializeUnspecified(thyraLinOp);
251}
252
253template <class SC,class LO,class GO,class NO>
254void Preconditioner<SC,LO,GO,NO>::buildPreconditioner( std::string type )
255{
256
257#ifdef PRECONDITIONER_TIMER
258 CommConstPtr_Type comm;
259 if (!problem_.is_null())
260 comm = problem_->getComm();
261 else if(!timeProblem_.is_null())
262 comm = timeProblem_->getComm();
263 comm->barrier();
264 Teuchos::TimeMonitor preconditionerTimeMonitor(*preconditionerTimer_);
265#endif
266
267 if (!type.compare("Monolithic")){
268 buildPreconditionerMonolithic( );
269 }
270 else if (!type.compare("Teko")){
271#ifdef FEDD_HAVE_TEKO
272 buildPreconditionerTeko( );
273#else
274 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Teko not found! Build Trilinos with Teko.");
275#endif
276 }
277 else if( type == "FaCSI" || type == "FaCSI-Teko" || type == "FaCSI-Block" ){
278 buildPreconditionerFaCSI( type );
279 }
280 else if(type == "Triangular" || type == "Diagonal" || type == "PCD" || type == "LSC"){
281 buildPreconditionerBlock2x2( );
282 }
283 else
284 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown preconditioner type.");
285
286#ifdef PRECONDITIONER_TIMER
287 comm->barrier();
288#endif
289}
290
291template <class SC,class LO,class GO,class NO>
292void Preconditioner<SC,LO,GO,NO>::buildPreconditionerMonolithic( )
293{
294
295 CommConstPtr_Type comm;
296 if (!problem_.is_null())
297 comm = problem_->getComm();
298 else if(!timeProblem_.is_null())
299 comm = timeProblem_->getComm();
300 else
301 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Preconditioner can not be used without a problem.");
302
303 bool verbose ( comm->getRank() == 0 );
304
305 ParameterListPtr_Type parameterList;
306
307 if (!problem_.is_null())
308 parameterList = problem_->getParameterList();
309 else if(!timeProblem_.is_null())
310 parameterList = timeProblem_->getParameterList();
311
312 std::string precType = parameterList->sublist("ThyraPreconditioner").get("Preconditioner Type", "FROSch");
313
314 bool useRepeatedMaps = parameterList->get( "Use repeated maps", true );
315 bool useNodeLists = parameterList->get( "Use node lists", true );
316 ParameterListPtr_Type pListThyraPrec = sublist( parameterList, "ThyraPreconditioner" );
317 ParameterListPtr_Type plFrosch = sublist( sublist( pListThyraPrec, "Preconditioner Types" ), "FROSch");
318
319 ThyraLinOpConstPtr_Type thyraMatrix;
320 if (!problem_.is_null())
321 thyraMatrix = problem_->getSystem()->getThyraLinOp();
322 else if(!timeProblem_.is_null())
323 thyraMatrix = timeProblem_->getSystemCombined()->getThyraLinOp();
324
325 UN numberOfBlocks = parameterList->get("Number of blocks",1);
326
327 typedef Tpetra::MultiVector<SC,LO,GO,NO> TMultiVector;
328 typedef Teuchos::RCP<TMultiVector> TMultiVectorPtr;
329 typedef Teuchos::ArrayRCP<TMultiVectorPtr> TMultiVectorPtrVecPtr;
330
331 // ------------------------
332 // Defs to cast back from tpetra to xpetra
333 typedef Xpetra::MultiVector<SC,LO,GO,NO> XMultiVector;
334 typedef Teuchos::RCP<XMultiVector> XMultiVectorPtr;
335 typedef Teuchos::ArrayRCP<XMultiVectorPtr> XMultiVectorPtrVecPtr;
336
337 typedef Xpetra::Map<LO,GO,NO> XpetraMap_Type;
338 typedef Teuchos::RCP<XpetraMap_Type> XpetraMapPtr_Type;
339 typedef Teuchos::RCP<const XpetraMap_Type> XpetraMapConstPtr_Type;
340 typedef const XpetraMapConstPtr_Type XpetraMapConstPtrConst_Type;
341 // ---------
342 // XMapVecPtrVecPtr
343 //Teuchos::ArrayRCP<Teuchos::RCP<Tpetra::Map<LO,GO,NO> > > repeatedMaps(numberOfBlocks);
344 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::Map<LO,GO,NO> > > repeatedMaps(numberOfBlocks);
345 // --------
346
347 XMultiVectorPtrVecPtr nodeListVec( numberOfBlocks );
348 if (!useNodeLists)
349 nodeListVec = Teuchos::null;
350
351 // timeProblem_->getSystemCombined()->print();
352 //Set Precondtioner lists
353 if (!precondtionerIsBuilt_) {
354 if ( !precType.compare("FROSch") ){
355 Teuchos::ArrayRCP<FROSch::DofOrdering> dofOrderings(numberOfBlocks);
356 Teuchos::ArrayRCP<UN> dofsPerNodeVector(numberOfBlocks);
357 for (UN i = 0; i < numberOfBlocks; i++) {
358 dofsPerNodeVector[i] = (UN) pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get( "DofsPerNode" + std::to_string(i+1), 1);
359
360 if (useRepeatedMaps) {
361 if (dofsPerNodeVector[i] > 1){
362
363 if (!problem_.is_null()) {
364 if (problem_->getDomain(i)->getFEType() == "P0") {
365 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Vector field map not implemented for P0 elements.");
366 }
367
368 //Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = problem_->getDomain(i)->getMapVecFieldRepeated()->getTpetraMap();
369 //Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
370
371 MapConstPtr_Type mapConstTmp = problem_->getDomain(i)->getMapVecFieldRepeated();//->getTpetraMap();
372 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
373 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
374
375 repeatedMaps[i] = mapX;
376
377
378 }
379 else if(!timeProblem_.is_null()){
380 if (timeProblem_->getDomain(i)->getFEType() == "P0") {
381 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Vector field map not implemented for P0 elements.");
382 }
383 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = timeProblem_->getDomain(i)->getMapVecFieldRepeated()->getTpetraMap();
384 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
385
386 MapConstPtr_Type mapConstTmp = timeProblem_->getDomain(i)->getMapVecFieldRepeated();//->getTpetraMap();
387 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
388 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
389
390 repeatedMaps[i] = mapX;
391 }
392 }
393 else{
394 if (!problem_.is_null()){
395 if (problem_->getDomain(i)->getFEType() == "P0") {
396 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = problem_->getDomain(i)->getElementMap()->getTpetraMap();
397 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
398 MapConstPtr_Type mapConstTmp = problem_->getDomain(i)->getElementMap();//->getTpetraMap();
399 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
400 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
401
402 repeatedMaps[i] = mapX;
403 }
404 else{
405 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = problem_->getDomain(i)->getMapRepeated()->getTpetraMap();
406 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
407 MapConstPtr_Type mapConstTmp = problem_->getDomain(i)->getMapRepeated();//->getTpetraMap();
408 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
409 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
410
411 repeatedMaps[i] = mapX;
412 }
413 }
414 else if (!timeProblem_.is_null()){
415 if (timeProblem_->getDomain(i)->getFEType() == "P0") {
416 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = timeProblem_->getDomain(i)->getElementMap()->getTpetraMap();
417 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
418 MapConstPtr_Type mapConstTmp = timeProblem_->getDomain(i)->getElementMap();//->getTpetraMap();
419 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
420 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
421
422 repeatedMaps[i] = mapX;
423 }
424 else{
425 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = timeProblem_->getDomain(i)->getMapRepeated()->getTpetraMap();
426 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
427 MapConstPtr_Type mapConstTmp = timeProblem_->getDomain(i)->getMapRepeated();//->getTpetraMap();
428 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
429 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
430
431 repeatedMaps[i] = mapX;
432 }
433 }
434 }
435 }
436
437 if (useNodeLists) {
438 if (!problem_.is_null()){
439 TEUCHOS_TEST_FOR_EXCEPTION( problem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Node lists cannot be used for P0 elements." );
440 Teuchos::RCP< Tpetra::MultiVector<SC,LO,GO,NO> > nodeListTpetra = problem_->getDomain(i)->getNodeListMV()->getTpetraMultiVectorNonConst();
441 Teuchos::RCP< Xpetra::TpetraMultiVector<SC,LO,GO,NO> > nodeListXpetraTpetra = Teuchos::rcp(new Xpetra::TpetraMultiVector<SC,LO,GO,NO>(nodeListTpetra));
442 Teuchos::RCP< Xpetra::MultiVector<SC,LO,GO,NO> > nodeListXpetra = Teuchos::rcp_dynamic_cast<Xpetra::MultiVector<SC,LO,GO,NO>>(nodeListXpetraTpetra);
443
444 nodeListVec[i] = nodeListXpetra; //problem_->getDomain(i)->getNodeListMV()->getTpetraMultiVectorNonConst();
445 }
446 else if (!timeProblem_.is_null()){
447 TEUCHOS_TEST_FOR_EXCEPTION( timeProblem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Node lists cannot be used for P0 elements." );
448 Teuchos::RCP< Tpetra::MultiVector<SC,LO,GO,NO> > nodeListTpetra = timeProblem_->getDomain(i)->getNodeListMV()->getTpetraMultiVectorNonConst();
449 Teuchos::RCP< Xpetra::TpetraMultiVector<SC,LO,GO,NO> > nodeListXpetraTpetra = Teuchos::rcp(new Xpetra::TpetraMultiVector<SC,LO,GO,NO>(nodeListTpetra));
450 Teuchos::RCP< Xpetra::MultiVector<SC,LO,GO,NO> > nodeListXpetra = Teuchos::rcp_dynamic_cast<Xpetra::MultiVector<SC,LO,GO,NO>>(nodeListXpetraTpetra);
451
452 nodeListVec[i] = nodeListXpetra;//timeProblem_->getDomain(i)->getNodeListMV()->getTpetraMultiVectorNonConst();
453 }
454
455 }
456
457
458 if (!pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get( "DofOrdering" + std::to_string(i+1), "NodeWise" ).compare("DimensionWise"))
459 dofOrderings[i] = FROSch::DimensionWise;
460 else if (!pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get( "DofOrdering" + std::to_string(i+1), "NodeWise" ).compare("NodeWise"))
461 dofOrderings[i] = FROSch::NodeWise;
462 else
463 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Please chose a valid DOF ordering for ThyraPreconditioner(FROSch).");
464 }
465 int dim;
466 if (!problem_.is_null())
467 dim = problem_->getDomain(0)->getDimension();
468 else if(!timeProblem_.is_null())
469 dim = timeProblem_->getDomain(0)->getDimension();
470
471 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Dimension", dim);
472 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Repeated Map Vector",repeatedMaps);
473 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Coordinates List Vector",nodeListVec);
474 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("DofOrdering Vector",dofOrderings);
475 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("DofsPerNode Vector",dofsPerNodeVector);
476 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Mpi Ranks Coarse",parameterList->sublist("General").get("Mpi Ranks Coarse",0) );
477
478 // This a pressure projection is only used for saddle point problems. We check here if we have a pressure projection set and if we have more than one block or one block with dim dof per node (i.e. fluid problem)
479 // This is unfortunately called pressure correction in FROSch, but is a projection!
480 if(!pressureProjection_.is_null() && ( dofsPerNodeVector.size() > 1 || dofsPerNodeVector[0] == 1) ){
481 pressureProjection_->merge(); // We merge the projection vector, as FROSch does not distinguish between blocks
482
483 Teuchos::RCP< Tpetra::MultiVector<SC,LO,GO,NO> > vectorTpetra = pressureProjection_->getMergedVectorNonConst()->getTpetraMultiVectorNonConst();
484 Teuchos::RCP< Xpetra::TpetraMultiVector<SC,LO,GO,NO> > vectorXpetraTpetra = Teuchos::rcp(new Xpetra::TpetraMultiVector<SC,LO,GO,NO>(vectorTpetra));
485 Teuchos::RCP< Xpetra::MultiVector<SC,LO,GO,NO> > vectorXpetra = Teuchos::rcp_dynamic_cast<Xpetra::MultiVector<SC,LO,GO,NO>>(vectorXpetraTpetra);
486
487 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").sublist("AlgebraicOverlappingOperator").set("Projection",vectorXpetra);
488 // In case of pressure correction we set the parameter in the parameterlist to true
489 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").sublist("AlgebraicOverlappingOperator").set("Use Pressure Correction", true);
490 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").sublist("AlgebraicOverlappingOperator").set("Use Local Pressure Correction", true);
491 }
492
493 /* We need to set the ranges of local problems and the coarse problem here.
494 When using an unstructured decomposition of, e.g., FSI, with 2 domains, which might be on a different set of ranks, we need to set the following parameters for FROSch here. Similarly we need to set a coarse rank problem range. For now, we use extra coarse ranks only for structured decompositions
495 */
496
497 int lowerBound = 100000000;
498 int upperBound = -1;
499 for (UN i = 0; i < numberOfBlocks; i++) {
500 tuple_intint_Type rankRange;
501 if (!problem_.is_null())
502 rankRange = problem_->getDomain(i)->getMesh()->getRankRange();
503 else if(!timeProblem_.is_null())
504 rankRange = timeProblem_->getDomain(i)->getMesh()->getRankRange();
505
506 if (std::get<0>(rankRange) < lowerBound)
507 lowerBound = std::get<0>(rankRange);
508 if (std::get<1>(rankRange) > upperBound)
509 upperBound = std::get<1>(rankRange);
510 }
511
512 int lowerBoundCoarse = lowerBound;
513 int upperBoundCoarse = upperBound;
514 // For now only for structured decompositions
515 if ( parameterList->sublist("General").get("Mpi Ranks Coarse",0) > 0 ){
516 lowerBoundCoarse = upperBound + 1;
517 upperBoundCoarse = comm->getSize() - 1;
518 }
519 if (verbose) {
520 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
521 std::cout << "\t --- Range for local problems of preconditioner from " << lowerBound << " to " << upperBound << std::endl;
522 std::cout << "\t --- Range for coarse problem of preconditioner from " << lowerBoundCoarse << " to " << upperBoundCoarse << std::endl;
523 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
524 }
525
526 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Local problem ranks lower bound", lowerBound );
527 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Local problem ranks upper bound", upperBound );
528
529 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Coarse problem ranks lower bound", lowerBoundCoarse );
530 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Coarse problem ranks upper bound", upperBoundCoarse );
531
532 // if ( !pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive").compare("Multiplicative") ){
533 // pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
534 // }
535
536 // Here we set the parameterlist which is a member variable of class Problem
537 ParameterListPtr_Type pListThyraSolver = sublist( parameterList, "ThyraSolver" ); //ch 11.02.19 To we need full Solver details
538 pListThyraSolver->setParameters( *pListThyraPrec );
539
540 LinSolverBuilderPtr_Type solverBuilder;
541 if (!problem_.is_null())
542 solverBuilder = problem_->getLinearSolverBuilder();
543 else if(!timeProblem_.is_null())
544 solverBuilder = timeProblem_->getLinearSolverBuilder();
545
546 // We save the pointer to the coarse matrix in this parameter list inside FROSch
547 pListPhiExport_ = pListThyraSolver;
548
549 solverBuilder->setParameterList(pListThyraSolver);
550 precFactory_ = solverBuilder->createPreconditioningStrategy("");
551
552 if ( thyraPrec_.is_null() ){
553 thyraPrec_ = precFactory_->createPrec();
554 }
555
556 Thyra::initializePrec<SC>(*precFactory_, thyraMatrix, thyraPrec_.ptr()); // (precfactory, fwdOp, prec) Problem: PreconditionerBase<SC>* thyraPrec_
557 precondtionerIsBuilt_ = true;
558
559 }
560 else
561 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unknown preconditioner type; You can only compute a FROSch preconditioner here." );
562 }
563 else {
564 TEUCHOS_TEST_FOR_EXCEPTION( precFactory_.is_null(), std::runtime_error, "precFactory_ is null.");
565 Thyra::initializePrec<SC>(*precFactory_, thyraMatrix, thyraPrec_.ptr());
566 }
567
568 if (parameterList->sublist( "Exporter" ).get("Export coarse functions",false))
569 exportCoarseBasis();
570}
571
572template <class SC,class LO,class GO,class NO>
573void Preconditioner<SC,LO,GO,NO>::buildPreconditionerMonolithicFSI( )
574{
575
576 CommConstPtr_Type comm;
577 if (!problem_.is_null())
578 comm = problem_->getComm();
579 else if(!timeProblem_.is_null())
580 comm = timeProblem_->getComm();
581 else
582 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Preconditioner can not be used without a problem.");
583
584 bool verbose ( comm->getRank() == 0 );
585
586 ParameterListPtr_Type parameterList;
587
588 if (!problem_.is_null())
589 parameterList = problem_->getParameterList();
590 else if(!timeProblem_.is_null())
591 parameterList = timeProblem_->getParameterList();
592
593 std::string precType = parameterList->sublist("ThyraPreconditioner").get("Preconditioner Type", "FROSch");
594
595 bool useRepeatedMaps = parameterList->get( "Use repeated maps", true );
596 bool useNodeLists = parameterList->get( "Use node lists", true );
597 ParameterListPtr_Type pListThyraPrec = sublist( parameterList, "ThyraPreconditioner" );
598 ParameterListPtr_Type plFrosch = sublist( sublist( pListThyraPrec, "Preconditioner Types" ), "FROSch");
599// setRecyclingParameter( plFrosch );
600
601
602 ThyraLinOpConstPtr_Type thyraMatrix;
603 if (!problem_.is_null())
604 thyraMatrix = problem_->getSystem()->getThyraLinOp();
605 else if(!timeProblem_.is_null())
606 thyraMatrix = timeProblem_->getSystemCombined()->getThyraLinOp();
607
608 UN numberOfBlocks = parameterList->get("Number of blocks",1);
609 TEUCHOS_TEST_FOR_EXCEPTION( numberOfBlocks<4 || numberOfBlocks>5, std::logic_error, "Unknown FSI size." );
610
611
612 typedef Tpetra::MultiVector<SC,LO,GO,NO> TMultiVector;
613 typedef Teuchos::RCP<TMultiVector> TMultiVectorPtr;
614 typedef Teuchos::ArrayRCP<TMultiVectorPtr> TMultiVectorPtrVecPtr;
615
616 // ------------------------
617 // Defs to cast back from tpetra to xpetra
618 typedef Xpetra::MultiVector<SC,LO,GO,NO> XMultiVector;
619 typedef Teuchos::RCP<XMultiVector> XMultiVectorPtr;
620 typedef Teuchos::ArrayRCP<XMultiVectorPtr> XMultiVectorPtrVecPtr;
621
622 typedef Xpetra::Map<LO,GO,NO> XpetraMap_Type;
623 typedef Teuchos::RCP<XpetraMap_Type> XpetraMapPtr_Type;
624 typedef Teuchos::RCP<const XpetraMap_Type> XpetraMapConstPtr_Type;
625 typedef const XpetraMapConstPtr_Type XpetraMapConstPtrConst_Type;
626
627 // XMapPtrVecPtr
628 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::Map<LO,GO,NO> > > repeatedMaps(numberOfBlocks);
629 // ------------------------
630
631 TMultiVectorPtrVecPtr nodeListVec( numberOfBlocks );
632 nodeListVec = Teuchos::null;
633
634 //Set Precondtioner lists
635 if (!precondtionerIsBuilt_) {
636 if ( !precType.compare("FROSch") ){
637 Teuchos::ArrayRCP<FROSch::DofOrdering> dofOrderings(numberOfBlocks);
638 Teuchos::ArrayRCP<UN> dofsPerNodeVector(numberOfBlocks);
639 for (UN i = 0; i < numberOfBlocks; i++) {
640 dofsPerNodeVector[i] = (UN) pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get( "DofsPerNode" + std::to_string(i+1), 1);
641 if (i==3) { //interface coupling
642 TEUCHOS_TEST_FOR_EXCEPTION( timeProblem_.is_null(), std::logic_error, "FSI time problem is null!" );
643 TEUCHOS_TEST_FOR_EXCEPTION( timeProblem_->getDomain(i)->getFEType() == "P0", std::logic_error, "We should not be able to use P0 for interface coupling." );
644 //Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp timeProblem_->getDomain(i)->getInterfaceMapUnique()->getTpetraMap();
645 //Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
646 MapConstPtr_Type mapConstTmp =timeProblem_->getDomain(i)->getInterfaceMapUnique();//->getTpetraMap();
647 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
648 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
649
650 repeatedMaps[i] = mapX;
651 }
652 else {
653 if (useRepeatedMaps) {
654 if (dofsPerNodeVector[i] > 1){
655
656 if (!problem_.is_null()) {
657 if (problem_->getDomain(i)->getFEType() == "P0") {
658 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Vector field map not implemented for P0 elements.");
659 }
660
661 //Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = problem_->getDomain(i)->getMapVecFieldRepeated()->getTpetraMap();
662 //Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
663
664 MapConstPtr_Type mapConstTmp = problem_->getDomain(i)->getMapVecFieldRepeated();//->getTpetraMap();
665 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
666 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
667
668 repeatedMaps[i] = mapX;
669
670 }
671 else if(!timeProblem_.is_null()){
672 if (timeProblem_->getDomain(i)->getFEType() == "P0") {
673 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Vector field map not implemented for P0 elements.");
674 }
675 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = timeProblem_->getDomain(i)->getMapVecFieldRepeated()->getTpetraMap();
676 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
677
678 MapConstPtr_Type mapConstTmp = timeProblem_->getDomain(i)->getMapVecFieldRepeated();//->getTpetraMap();
679 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
680 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
681
682 repeatedMaps[i] = mapX;
683 }
684 }
685 else{
686 if (!problem_.is_null()){
687 if (problem_->getDomain(i)->getFEType() == "P0") {
688 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = problem_->getDomain(i)->getElementMap()->getTpetraMap();
689 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
690 MapConstPtr_Type mapConstTmp = problem_->getDomain(i)->getElementMap();//->getTpetraMap();
691 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
692 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
693
694
695
696 repeatedMaps[i] = mapX;
697 }
698 else{
699 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = problem_->getDomain(i)->getMapRepeated()->getTpetraMap();
700 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
701 MapConstPtr_Type mapConstTmp = problem_->getDomain(i)->getMapRepeated();//->getTpetraMap();
702 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
703 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
704
705 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
706 mapX->describe(*out,Teuchos::VERB_EXTREME);
707 repeatedMaps[i] = mapX;
708 }
709 }
710 else if (!timeProblem_.is_null()){
711 if (timeProblem_->getDomain(i)->getFEType() == "P0") {
712 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = timeProblem_->getDomain(i)->getElementMap()->getTpetraMap();
713 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
714 MapConstPtr_Type mapConstTmp = timeProblem_->getDomain(i)->getElementMap();//->getTpetraMap();
715 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
716 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
717
718 repeatedMaps[i] = mapX;
719 }
720 else{
721 // Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp = timeProblem_->getDomain(i)->getMapRepeated()->getTpetraMap();
722 // Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
723 MapConstPtr_Type mapConstTmp = timeProblem_->getDomain(i)->getMapRepeated();//->getTpetraMap();
724 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
725 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
726
727 repeatedMaps[i] = mapX;
728 }
729 }
730 }
731 }
732
733 if (!pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get( "DofOrdering" + std::to_string(i+1), "NodeWise" ).compare("DimensionWise"))
734 dofOrderings[i] = FROSch::DimensionWise;
735 else if (!pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get( "DofOrdering" + std::to_string(i+1), "NodeWise" ).compare("NodeWise"))
736 dofOrderings[i] = FROSch::NodeWise;
737 else
738 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Please chose a valid DOF ordering for ThyraPreconditioner(FROSch).");
739
740 }
741 }
742 int dim;
743 if (!problem_.is_null())
744 dim = problem_->getDomain(0)->getDimension();
745 else if(!timeProblem_.is_null())
746 dim = timeProblem_->getDomain(0)->getDimension();
747
748 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Dimension", dim);
749 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Repeated Map Vector",repeatedMaps);
750 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Coordinates List Vector",nodeListVec);
751 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("DofOrdering Vector",dofOrderings);
752 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("DofsPerNode Vector",dofsPerNodeVector);
753 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Mpi Ranks Coarse",parameterList->sublist("General").get("Mpi Ranks Coarse",0) );
754
755 /* We need to set the ranges of local problems and the coarse problem here.
756 When using an unstructured decomposition of, e.g., FSI, with 2 domains, which might be on a different set of ranks, we need to set the following parameters for FROSch here. Similarly we need to set a coarse rank problem range. For now, we use extra coarse ranks only for structured decompositions
757 */
758
759 int lowerBound = 100000000;
760 int upperBound = -1;
761 for (UN i = 0; i < numberOfBlocks; i++) {
762 tuple_intint_Type rankRange;
763 if (!problem_.is_null())
764 rankRange = problem_->getDomain(i)->getMesh()->getRankRange();
765 else if(!timeProblem_.is_null())
766 rankRange = timeProblem_->getDomain(i)->getMesh()->getRankRange();
767
768 if (std::get<0>(rankRange) < lowerBound)
769 lowerBound = std::get<0>(rankRange);
770 if (std::get<1>(rankRange) > upperBound)
771 upperBound = std::get<1>(rankRange);
772 }
773
774 int lowerBoundCoarse = lowerBound;
775 int upperBoundCoarse = upperBound;
776 // For now only for structured decompositions
777 if ( parameterList->sublist("General").get("Mpi Ranks Coarse",0) > 0 ){
778 lowerBoundCoarse = upperBound + 1;
779 upperBoundCoarse = comm->getSize() - 1;
780 }
781 if (verbose) {
782 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
783 std::cout << "\t --- Range for local problems of preconditioner from " << lowerBound << " to " << upperBound << std::endl;
784 std::cout << "\t --- Range for coarse problem of preconditioner from " << lowerBoundCoarse << " to " << upperBoundCoarse << std::endl;
785 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
786 }
787
788 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Local problem ranks lower bound", lowerBound );
789 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Local problem ranks upper bound", upperBound );
790
791 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Coarse problem ranks lower bound", lowerBoundCoarse );
792 pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set( "Coarse problem ranks upper bound", upperBoundCoarse );
793
794 // if ( !pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive").compare("Multiplicative") ){
795 // pListThyraPrec->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
796 // }
797
798 // Here we set the parameterlist which is a member variable of class Problem
799 ParameterListPtr_Type pListThyraSolver = sublist( parameterList, "ThyraSolver" ); //ch 11.02.19 To we need full Solver details
800 pListThyraSolver->setParameters( *pListThyraPrec );
801
802 LinSolverBuilderPtr_Type solverBuilder;
803 if (!problem_.is_null())
804 solverBuilder = problem_->getLinearSolverBuilder();
805 else if(!timeProblem_.is_null())
806 solverBuilder = timeProblem_->getUnderlyingProblem()->getLinearSolverBuilder();
807
808 // We save the pointer to the coarse matrix in this parameter list inside FROSch
809 pListPhiExport_ = pListThyraSolver;
810
811 solverBuilder->setParameterList(pListThyraSolver);
812 precFactory_ = solverBuilder->createPreconditioningStrategy("");
813
814 if ( thyraPrec_.is_null() )
815 thyraPrec_ = precFactory_->createPrec();
816
817
818 Thyra::initializePrec<SC>(*precFactory_, thyraMatrix, thyraPrec_.ptr());
819 precondtionerIsBuilt_ = true;
820
821 }
822 else
823 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unknown preconditioner type; You can only compute a FROSch preconditioner here." );
824 }
825 else {
826 TEUCHOS_TEST_FOR_EXCEPTION( precFactory_.is_null(), std::runtime_error, "precFactory_ is null.");
827 Thyra::initializePrec<SC>(*precFactory_, thyraMatrix, thyraPrec_.ptr());
828 }
829
830 if (parameterList->sublist( "Exporter" ).get("Export coarse functions",false))
831 exportCoarseBasisFSI();
832
833}
834
835
836#ifdef FEDD_HAVE_TEKO
837template <class SC,class LO,class GO,class NO>
838void Preconditioner<SC,LO,GO,NO>::buildPreconditionerTeko( )
839{
840
841 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
842
843 ParameterListPtr_Type parameterList;
844 if (!problem_.is_null())
845 parameterList = problem_->getParameterList();
846 else if(!timeProblem_.is_null())
847 parameterList = timeProblem_->getParameterList();
848 else
849 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Preconditioner can not be used without a problem.");
850
851
852 LinSolverBuilderPtr_Type solverBuilder;
853 if (!problem_.is_null())
854 solverBuilder = problem_->getLinearSolverBuilder();
855 else if(!timeProblem_.is_null())
856 solverBuilder = timeProblem_->getUnderlyingProblem()->getLinearSolverBuilder();
857
858 ParameterListPtr_Type tekoPList= sublist( parameterList, "Teko Parameters" );
859
860 if (precFactory_.is_null()) {
861 ParameterListPtr_Type tmpSubList = sublist( sublist( sublist( sublist( parameterList, "Teko Parameters" ) , "Preconditioner Types" ) , "Teko" ) , "Inverse Factory Library" );
862
863 //only sets repeated maps in parameterlist if FROSch is used for both block
864 if ( !tmpSubList->sublist("FROSch-Pressure").get("Type","FROSch").compare("FROSch") &&
865 !tmpSubList->sublist("FROSch-Velocity").get("Type","FROSch").compare("FROSch") ) {
866 setVelocityParameters( tekoPList, parameterList->sublist("General").get("Mpi Ranks Coarse",0) );
867 setPressureParameters( tekoPList, parameterList->sublist("General").get("Mpi Ranks Coarse",0) );
868 }
869 }
870
871 BlockMatrixPtr_Type system;
872 if (!problem_.is_null())
873 system = problem_->getSystem();
874 else if(!timeProblem_.is_null())
875 system = timeProblem_->getSystemCombined();
876
877 TEUCHOS_TEST_FOR_EXCEPTION( system->size()!=2, std::logic_error, "Wrong size of system for Teko-Block-Preconditioners.");
878
879 Teko::LinearOp thyraF = system->getBlock(0,0)->getThyraLinOp();
880 Teko::LinearOp thyraBT = system->getBlock(0,1)->getThyraLinOp();
881 Teko::LinearOp thyraB = system->getBlock(1,0)->getThyraLinOp();
882
883 if (!system->blockExists(1,1)){
884 MatrixPtr_Type dummy;
885 dummy.reset( new Matrix_Type( system->getBlock(1,0)->getMap(), 1 ) );
886 dummy->fillComplete();
887 system->addBlock( dummy, 1, 1 );
888 }
889 Teko::LinearOp thyraC = system->getBlock(1,1)->getThyraLinOp();
890
891 tekoLinOp_ = Thyra::block2x2(thyraF,thyraBT,thyraB,thyraC);
892
893 if (!precondtionerIsBuilt_) {
894 if ( precFactory_.is_null() ){
895 ParameterListPtr_Type pListThyraSolver = sublist( parameterList, "ThyraSolver" );
896
897 pListThyraSolver->setParameters( *tekoPList );
898
899 solverBuilder->setParameterList( pListThyraSolver );
900 precFactory_ = solverBuilder->createPreconditioningStrategy("");//createPreconditioningStrategy(*solverBuilder);
901
902 rh_.reset(new Teko::RequestHandler());
903
904 if(!tekoPList->sublist("Preconditioner Types").sublist("Teko").get("Inverse Type", "SIMPLE").compare("LSC") ||
905 !tekoPList->sublist("Preconditioner Types").sublist("Teko").get("Inverse Type", "SIMPLE").compare("LSC-Pressure-Laplace") ||
906 !tekoPList->sublist("Preconditioner Types").sublist("Teko").get("Inverse Type", "SIMPLE").compare("SIMPLE")){
907 Teko::LinearOp thyraMass = velocityMassMatrix_;
908 Teuchos::RCP< Teko::StaticRequestCallback<Teko::LinearOp> > callbackMass = Teuchos::rcp(new Teko::StaticRequestCallback<Teko::LinearOp> ( "Velocity Mass Matrix", thyraMass ) );
909 rh_->addRequestCallback( callbackMass );
910
911 Teko::LinearOp thyraLaplace = pressureLaplace_;
912
913 Teuchos::RCP< Teko::StaticRequestCallback<Teko::LinearOp> > callbackLaplace = Teuchos::rcp(new Teko::StaticRequestCallback<Teko::LinearOp> ( "Pressure Laplace Operator", thyraLaplace ) );
914 rh_->addRequestCallback( callbackLaplace );
915 }
916 else if(!tekoPList->sublist("Preconditioner Types").sublist("Teko").get("Inverse Type", "SIMPLE").compare("PCD")){
917
918 // Velocity Mass Matrix
919 Teko::LinearOp thyraMass = velocityMassMatrix_;
920 Teuchos::RCP< Teko::StaticRequestCallback<Teko::LinearOp> > callbackMass = Teuchos::rcp(new Teko::StaticRequestCallback<Teko::LinearOp> ( "Velocity Mass Matrix", thyraMass ) );
921 rh_->addRequestCallback( callbackMass );
922
923 // Pressure Laplace
924 Teko::LinearOp thyraLaplace = pressureLaplace_;
925 Teuchos::RCP< Teko::StaticRequestCallback<Teko::LinearOp> > callbackLaplace = Teuchos::rcp(new Teko::StaticRequestCallback<Teko::LinearOp> ( "Pressure Laplace Operator", thyraLaplace ) );
926 rh_->addRequestCallback( callbackLaplace );
927
928 // Pressure Mass
929 Teko::LinearOp thyraPressureMass = pressureMass_;
930 Teuchos::RCP< Teko::StaticRequestCallback<Teko::LinearOp> > callbackPressureMass = Teuchos::rcp(new Teko::StaticRequestCallback<Teko::LinearOp> ( "Pressure Mass Matrix", thyraPressureMass ) );
931 rh_->addRequestCallback( callbackPressureMass );
932
933 // PCD
934 if (!timeProblem_.is_null()){
935 timeProblem_->assemble("UpdateConvectionDiffusionOperator");
936 }
937 else{
938 problem_->assemble("UpdateConvectionDiffusionOperator");
939 }
940 Teko::LinearOp thyraPCD = pcdOperator_;
941 callbackPCD_ = Teuchos::rcp(new Teko::StaticRequestCallback<Teko::LinearOp> ( "PCD Operator", thyraPCD ) );
942 rh_->addRequestCallback( callbackPCD_ );
943
944 }
945 Teuchos::RCP< Teko::StratimikosFactory > tekoFactory = Teuchos::rcp_dynamic_cast<Teko::StratimikosFactory>(precFactory_);
946 tekoFactory->setRequestHandler( rh_ );
947 }
948
949
950 if ( thyraPrec_.is_null() ){
951 thyraPrec_ = precFactory_->createPrec();
952 }
953
954 Teuchos::RCP< const Thyra::DefaultLinearOpSource< SC > > thyraMatrixSourceOp = defaultLinearOpSource (tekoLinOp_);
955 // Thyra::initializePrec<SC>(*precFactory, thyraMatrixSourceOp, thyraPrec_.ptr());
956
957 precFactory_->initializePrec(thyraMatrixSourceOp, thyraPrec_.get());
958 precondtionerIsBuilt_ = true;
959
960 }
961 else{
962 if(!tekoPList->sublist("Preconditioner Types").sublist("Teko").get("Inverse Type", "SIMPLE").compare("PCD")){
963 // PCD: As part of the pcd preconditioner depends on the current velocity, we need to update it in each iteration.
964 //pcdOperatorMatrixPtr_->print();
965 // PCD
966 if (!timeProblem_.is_null()){
967 timeProblem_->assemble("UpdateConvectionDiffusionOperator");
968 }
969 else{
970 problem_->assemble("UpdateConvectionDiffusionOperator");
971 }
972 Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler());
973
974 // PCD
975 Teko::LinearOp thyraPCD;
976
977 // When we deal with a time problem only the underlying problem holds the updated pcd matrix. Why: When we assemble the PCD operator it within i.e. the Navier Stokes class.
978 // The Navier stokes class is derived from a nonlinear problem originally derived from a problem which has a precondidioner object to which the PCD operator is added in each reassembly.
979 // Since the Navier Stokes problem does not know the timeProblem, we need to extract the information from the problem which is added to the time problem / which the timeProblem is based on.
980 if(!timeProblem_.is_null())
981 {
982 // timeProblem_->getUnderlyingProblem()->preconditioner_->getPCDOperatorMatrix()->print();
983 pcdOperator_ = timeProblem_->getUnderlyingProblem()->preconditioner_->getPCDOperatorMatrix()->getThyraLinOp();
984 thyraPCD = timeProblem_->getUnderlyingProblem()->preconditioner_->getPCDOperatorMatrix()->getThyraLinOp();
985 }
986 else
987 thyraPCD= pcdOperator_;
988
989 // Updating matrix in the pointer
990 Teuchos::RCP< Teko::StaticRequestCallback<Teko::LinearOp> > callbackTmp = Teuchos::rcp(new Teko::StaticRequestCallback<Teko::LinearOp> ( "PCD Operator", thyraPCD ) );
991 *callbackPCD_ = *callbackTmp;
992
993 }
994
995 Teuchos::RCP< const Thyra::DefaultLinearOpSource< SC > > thyraMatrixSourceOp = defaultLinearOpSource (tekoLinOp_);
996 // Thyra::initializePrec<SC>(*precFactory, thyraMatrixSourceOp, thyraPrec_.ptr());
997
998 precFactory_->initializePrec(thyraMatrixSourceOp, thyraPrec_.get());
999 }
1000
1001}
1002
1003template <class SC,class LO,class GO,class NO>
1004void Preconditioner<SC,LO,GO,NO>::buildPreconditionerFaCSI( std::string type )
1005{
1006 typedef Domain<SC,LO,GO,NO> Domain_Type;
1007 typedef Teuchos::RCP<const Domain_Type> DomainConstPtr_Type;
1008 typedef std::vector<DomainConstPtr_Type> DomainConstPtr_vec_Type;
1009
1010 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
1011 // here we assume that FSI is always a time problem, this can be
1012 ParameterListPtr_Type parameterList;
1013 if (!timeProblem_.is_null())
1014 parameterList = timeProblem_->getParameterList();
1015 else
1016 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Preconditioner can not be used without a time problem.");
1017
1018 // Get FSI problem
1019 ProblemPtr_Type steadyProblem = timeProblem_->getUnderlyingProblem();
1020 Teuchos::RCP< FSI<SC,LO,GO,NO> > steadyFSI = Teuchos::rcp_dynamic_cast<FSI<SC,LO,GO,NO> >(steadyProblem);
1021 BlockMatrixPtr_Type fsiSystem = timeProblem_->getSystemCombined();
1022
1023 ParameterListPtr_Type pLFluid = steadyFSI->getFluidProblem()->getParameterList();
1024
1025
1026 std::string precTypeFluid;
1027 if (type == "FaCSI")
1028 precTypeFluid = "Monolithic";
1029 else if (type == "FaCSI-Teko"){
1030 precTypeFluid = "Teko";
1031 }
1032 else if (type == "FaCSI-Block"){
1033 precTypeFluid = parameterList->sublist("Parameter Fluid").get("Preconditioner Type", "PCD");
1034 }
1035
1036 CommConstPtr_Type comm = timeProblem_->getComm();
1037 bool useFluidPreconditioner = parameterList->sublist("General").get("Use Fluid Preconditioner", true);
1038 bool useSolidPreconditioner = parameterList->sublist("General").get("Use Solid Preconditioner", true);
1039 bool onlyDiagonal = parameterList->sublist("General").get("Only Diagonal", false);
1040 Teuchos::RCP< PrecOpFaCSI<SC,LO,GO,NO> > facsi
1041 = Teuchos::rcp(new PrecOpFaCSI<SC,LO,GO,NO> ( comm, precTypeFluid == "Monolithic", useFluidPreconditioner, useSolidPreconditioner, onlyDiagonal) );
1042
1043 if (comm->getRank() == 0) {
1044 if (onlyDiagonal)
1045 std::cout << "\t### No preconditioner will be used! ###" << std::endl;
1046 else
1047 std::cout << "\t### FaCSI standard ###" << std::endl;
1048 }
1049
1050 BlockMatrixPtr_Type fluidSystem = Teuchos::rcp( new BlockMatrix_Type(2) );
1051
1052 // We build copies of the fluid system with homogenous Dirichlet boundary conditions on the interface
1053 // MatrixPtr_Type f = Teuchos::rcp(new Matrix_Type( fsiSystem->getBlock(0,0) ) );
1054
1055 // MatrixPtr_Type bt = Teuchos::rcp(new Matrix_Type( fsiSystem->getBlock(0,1) ) );
1056 // MatrixPtr_Type b = Teuchos::rcp(new Matrix_Type( fsiSystem->getBlock(1,0) ) );
1057 // MatrixPtr_Type c;
1058 // if ( fsiSystem->blockExists(1,1) )
1059 // c = Teuchos::rcp(new Matrix_Type( fsiSystem->getBlock(1,1) ) );
1060 // fluidSystem->addBlock( f, 0, 0 );
1061 // fluidSystem->addBlock( bt, 0, 1 );
1062 // fluidSystem->addBlock( b, 1, 0 );
1063 // if ( fsiSystem->blockExists(1,1) )
1064 // fluidSystem->addBlock( c, 1, 1 );
1065
1066
1067 // We want to use the underlying Navier-Stokes Fluid Problem to build the preconditioner
1068 // We start with the fluid time problem
1069 Teuchos::RCP< TimeProblem<SC,LO,GO,NO> > fluidProblem = steadyFSI->problemTimeFluid_;
1070 fluidProblem->combineSystems(); // Build combined system || check if even is neccesary
1071 fluidProblem->setBoundariesSystem(); // Set boundaries || might also need fsi bc
1072 // The we cast the timeproblem to original Navier-Stokes problem and use it to build preconditioner
1073 Teuchos::RCP< NavierStokes<SC,LO,GO,NO> > fluidProblemSteady = Teuchos::rcp_dynamic_cast<NavierStokes<SC,LO,GO,NO> >(fluidProblem->getUnderlyingProblem());
1074
1075 faCSIBCFactory_->setSystem( fluidProblem->getSystemCombined() );
1076
1077 fluidProblemSteady->setupPreconditioner( precTypeFluid );
1078 precFluid_ = fluidProblemSteady->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1079
1080 //Setup structure problem
1081 bool nonlinearStructure = false;
1082 if (steadyFSI->getStructureProblem().is_null())
1083 nonlinearStructure = true;
1084
1085 ParameterListPtr_Type pLStructure;
1086 if (nonlinearStructure)
1087 pLStructure = steadyFSI->getNonLinStructureProblem()->getParameterList();
1088 else
1089 pLStructure = steadyFSI->getStructureProblem()->getParameterList();
1090
1091 if (probSolid_.is_null()){
1092 probSolid_ = Teuchos::rcp( new MinPrecProblem_Type( pLStructure, timeProblem_->getComm() ) );
1093 DomainConstPtr_vec_Type structDomain;
1094 if (nonlinearStructure)
1095 structDomain = steadyFSI->getNonLinStructureProblem()->getDomainVector();
1096 else
1097 structDomain = steadyFSI->getStructureProblem()->getDomainVector();
1098
1099 probSolid_->initializeDomains( structDomain );
1100 probSolid_->initializeLinSolverBuilder( timeProblem_->getLinearSolverBuilder() );
1101 }
1102
1103 BlockMatrixPtr_Type structSystem = Teuchos::rcp( new BlockMatrix_Type(1) );
1104
1105 structSystem->addBlock( fsiSystem->getBlock(2,2), 0, 0 );
1106
1107 probSolid_->initializeSystem( structSystem );
1108
1109 probSolid_->setupPreconditioner( );
1110
1111 precStruct_ = probSolid_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1112
1113
1114 //Setup geometry problem
1115
1116 if (timeProblem_->getSystem()->size()>4) {
1117 ParameterListPtr_Type pLGeometry = steadyFSI->getGeometryProblem()->getParameterList();
1118 if (probGeo_.is_null()) {
1119 probGeo_ = Teuchos::rcp( new MinPrecProblem_Type( pLGeometry, timeProblem_->getComm() ) );
1120 DomainConstPtr_vec_Type geoDomain = steadyFSI->getGeometryProblem()->getDomainVector();
1121 probGeo_->initializeDomains( geoDomain );
1122 probGeo_->initializeLinSolverBuilder( timeProblem_->getLinearSolverBuilder() );
1123 }
1124
1125 BlockMatrixPtr_Type geoSystem = Teuchos::rcp( new BlockMatrix_Type(1) );
1126
1127 geoSystem->addBlock( fsiSystem->getBlock(4,4), 0, 0 );
1128
1129 probGeo_->initializeSystem( geoSystem );
1130
1131 probGeo_->setupPreconditioner( );
1132
1133 precGeo_ = probGeo_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1134 }
1135 bool shape = false;
1136 if (fsiSystem->size()>4) {
1137 if (shape){
1138 facsi->setGIShape( fsiSystem->getBlock(3,0)->getThyraLinOpNonConst()/*C1*/,
1139 fsiSystem->getBlock(0,3)->getThyraLinOpNonConst()/*C1T*/,
1140 fsiSystem->getBlock(3,2)->getThyraLinOpNonConst()/*C2*/,
1141 fsiSystem->getBlock(4,2)->getThyraLinOpNonConst()/*C4*/,
1142 precStruct_,
1143 precFluid_,
1144 fsiSystem->getBlock(0,0)->getThyraLinOpNonConst()/*fF*/,
1145 fsiSystem->getBlock(0,1)->getThyraLinOpNonConst()/*fBT*/,
1146 precGeo_,
1147 fsiSystem->getBlock(0,4)->getThyraLinOpNonConst()/*shape v*/,
1148 fsiSystem->getBlock(1,4)->getThyraLinOpNonConst()/*shape p*/ );
1149 }
1150 else {
1151 facsi->setGI( fsiSystem->getBlock(3,0)->getThyraLinOpNonConst()/*C1*/,
1152 fsiSystem->getBlock(0,3)->getThyraLinOpNonConst()/*C1T*/,
1153 fsiSystem->getBlock(3,2)->getThyraLinOpNonConst()/*C2*/,
1154 fsiSystem->getBlock(4,2)->getThyraLinOpNonConst()/*C4*/,
1155 precStruct_,
1156 precFluid_,
1157 fsiSystem->getBlock(0,0)->getThyraLinOpNonConst()/*fF*/,
1158 fsiSystem->getBlock(0,1)->getThyraLinOpNonConst()/*fBT*/,
1159 precGeo_ );
1160 }
1161 }
1162 else {
1163 facsi->setGE( fsiSystem->getBlock(3,0)->getThyraLinOpNonConst()/*C1*/,
1164 fsiSystem->getBlock(0,3)->getThyraLinOpNonConst()/*C1T*/,
1165 fsiSystem->getBlock(3,2)->getThyraLinOpNonConst()/*C2*/,
1166 precStruct_,
1167 precFluid_,
1168 fsiSystem->getBlock(0,0)->getThyraLinOpNonConst()/*fF*/,
1169 fsiSystem->getBlock(0,1)->getThyraLinOpNonConst()/*fBT*/ );
1170 }
1171
1172 LinSolverBuilderPtr_Type solverBuilder = timeProblem_->getUnderlyingProblem()->getLinearSolverBuilder();
1173
1174 if ( precFactory_.is_null() )
1175 precFactory_ = solverBuilder->createPreconditioningStrategy("");
1176
1177 if ( thyraPrec_.is_null() )
1178 thyraPrec_ = precFactory_->createPrec();
1179
1180 Teuchos::RCP< Thyra::DefaultPreconditioner<SC> > defaultPrec =
1181 Teuchos::rcp_dynamic_cast< Thyra::DefaultPreconditioner<SC> > (thyraPrec_);
1182 ThyraLinOpPtr_Type linOp =
1183 Teuchos::rcp_dynamic_cast< Thyra::LinearOpBase<SC> > (facsi);
1184
1185 defaultPrec->initializeUnspecified( linOp );
1186
1187 precondtionerIsBuilt_ = true;
1188
1189}
1190
1191template <class SC,class LO,class GO,class NO>
1192void Preconditioner<SC,LO,GO,NO>::setPressureMassMatrix(MatrixPtr_Type massMatrix) const{
1193 pressureMassMatrixPtr_ = massMatrix;
1194 pressureMass_= massMatrix->getThyraLinOp();
1195}
1196
1197template <class SC,class LO,class GO,class NO>
1198void Preconditioner<SC,LO,GO,NO>::setPressureLaplaceMatrix(MatrixPtr_Type matrix) const{
1199 pressureLaplace_ =matrix->getThyraLinOp();
1200 pressureLaplaceMatrixPtr_ = matrix;
1201}
1202
1203// template <class SC,class LO,class GO,class NO>
1204// void Preconditioner<SC,LO,GO,NO>::setPressureMass(MatrixPtr_Type matrix) const{
1205// pressureMass_ = matrix->getThyraLinOp();
1206// pressureMassMatrixPtr_ = matrix;
1207// }
1208
1209template <class SC,class LO,class GO,class NO>
1210void Preconditioner<SC,LO,GO,NO>::setPCDOperator(MatrixPtr_Type matrix) const{
1211 pcdOperator_ = matrix->getThyraLinOp();
1212 pcdOperatorMatrixPtr_ = matrix;
1213}
1214
1215// Function to build a general 2 x 2 Block preconditioner
1216// Currently only used for (Navier-)Stokes type problems
1217// This includes the
1218// - Diagonal Prec, where the Schur complement is replaced by - 1/nu M_p
1219// - Triangular Prec, where the Schur complement is replaced by - 1/nu M_p
1220// - PCD Prec, where the Schur complement is replaced by -M_p F_p^-1 A_p
1221// - LSC Prec, where the Schur complement is replaced by -A_p^-1 (B (M_v^-1) F (M_v^-1) B^T ) A_p^-1
1222
1223template <class SC,class LO,class GO,class NO>
1224void Preconditioner<SC,LO,GO,NO>::buildPreconditionerBlock2x2( )
1225{
1226 PRECONDITIONER_START(buildPreconditionerBlock2x2, " buildPreconditionerBlock2x2");
1227
1228 typedef Domain<SC,LO,GO,NO> Domain_Type;
1229 typedef Teuchos::RCP<const Domain_Type> DomainConstPtr_Type;
1230 typedef std::vector<DomainConstPtr_Type> DomainConstPtr_vec_Type;
1231
1232 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
1233
1234 ParameterListPtr_Type parameterList;
1235 BlockMatrixPtr_Type system;
1236 CommConstPtr_Type comm;
1237 ProblemPtr_Type steadyProblem;
1238 if (!timeProblem_.is_null()){
1239 std::cout << "buildPreconditionerBlock2x2 timeProblem_ " << std::endl;
1240 parameterList = timeProblem_->getParameterList();
1241 system = timeProblem_->getSystemCombined();
1242 comm = timeProblem_->getComm();
1243 steadyProblem = timeProblem_->getUnderlyingProblem();
1244 }
1245 else{
1246 parameterList = problem_->getParameterList();
1247 system = problem_->getSystem();
1248 comm = problem_->getComm();
1249 steadyProblem = problem_;
1250 }
1251
1252 bool verbose( comm->getRank() == 0 );
1253
1254 if(verbose){
1255 std::cout << " ######################## " << std::endl;
1256 std::cout << " Build 2x2 Preconditioner " << std::endl;
1257 std::cout << " ######################## " << std::endl;
1258 }
1259 ParameterListPtr_Type plVelocity( new Teuchos::ParameterList( parameterList->sublist("Velocity preconditioner") ) );
1260 ParameterListPtr_Type plSchur( new Teuchos::ParameterList( parameterList->sublist("Schur complement preconditioner") ) );
1261
1262 //Addig General information to both lists
1263 plVelocity->sublist("General").setParameters(parameterList->sublist("General"));
1264 plSchur->sublist("General").setParameters(parameterList->sublist("General"));
1265
1266 Teuchos::RCP< PrecBlock2x2<SC,LO,GO,NO> > blockPrec2x2
1267 = Teuchos::rcp(new PrecBlock2x2<SC,LO,GO,NO> ( comm ) );
1268
1269 // The velocity problem is always treated the same
1270 if (probVelocity_.is_null()){
1271 probVelocity_ = Teuchos::rcp( new MinPrecProblem_Type( plVelocity, comm ) );
1272 DomainConstPtr_vec_Type domain1(0);
1273 domain1.push_back( steadyProblem->getDomain(0) );
1274 probVelocity_->initializeDomains( domain1 );
1275 probVelocity_->initializeLinSolverBuilder( steadyProblem->getLinearSolverBuilder() );
1276 }
1277
1278 BlockMatrixPtr_Type system1 = Teuchos::rcp( new BlockMatrix_Type(1) );
1279
1280 MatrixPtr_Type F = system->getBlock(0,0);//Teuchos::rcp(new Matrix_Type( system->getBlock(0,0) ) );
1281
1282 system1->addBlock( F, 0, 0 );
1283
1284 probVelocity_->initializeSystem( system1 );
1285
1286 PRECONDITIONER_START(setupFInv, " Setup Preconditioner for F");
1287 probVelocity_->setupPreconditioner( "Monolithic" ); // single matrix
1288 PRECONDITIONER_STOP(setupFInv);
1289
1290 precVelocity_ = probVelocity_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1291
1292 if (probSchur_.is_null()) {
1293 if(!timeProblem_.is_null()){
1294 probSchur_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1295 DomainConstPtr_vec_Type domain2(0);
1296 domain2.push_back( timeProblem_->getDomain(1) );
1297 probSchur_->initializeDomains( domain2 );
1298 probSchur_->initializeLinSolverBuilder( timeProblem_->getLinearSolverBuilder() );
1299 }
1300 else
1301 {
1302 probSchur_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1303 DomainConstPtr_vec_Type domain2(0);
1304 domain2.push_back( problem_->getDomain(1) );
1305 probSchur_->initializeDomains( domain2 );
1306 probSchur_->initializeLinSolverBuilder( problem_->getLinearSolverBuilder() );
1307
1308 }
1309 }
1310
1311 std::string type = parameterList->sublist("General").get("Preconditioner Method","Diagonal");
1312
1313 // We distinguish for the Schur complement component
1314 // Setup additional things
1315 PRECONDITIONER_START(setupSInv, " Setup Preconditioner for S");
1316
1317 if(type == "Diagonal" || type == "Triangular"){
1318 BlockMatrixPtr_Type Mp = Teuchos::rcp( new BlockMatrix_Type(1) );
1319
1320 Mp->addBlock( pressureMassMatrixPtr_, 0, 0 );
1321
1322 probSchur_->initializeSystem( Mp );
1323
1324 probSchur_->setupPreconditioner( "Monolithic" ); // single matrix
1325
1326 precSchur_ = probSchur_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1327 }
1328 else if (type == "PCD") {
1329 // For PCD we additionally need to setup the monolithic preconditioner for the
1330 // Laplace operator and the pressure mass matrix
1331 // We include a diagonal inverse approximation of the mass matrix
1332 if (verbose) {
1333 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1334 std::cout << "\t --- Building PCD Operator Components " << std::endl;
1335 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1336 }
1337
1338 if (probLaplace_.is_null()) {
1339 if (!timeProblem_.is_null()){
1340 probLaplace_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1341 DomainConstPtr_vec_Type domain2(0);
1342 domain2.push_back( timeProblem_->getDomain(1) );
1343 probLaplace_->initializeDomains( domain2 );
1344 probLaplace_->initializeLinSolverBuilder( timeProblem_->getLinearSolverBuilder() );
1345 }
1346 else{
1347 probLaplace_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1348 DomainConstPtr_vec_Type domain2(0);
1349 domain2.push_back( problem_->getDomain(1) );
1350 probLaplace_->initializeDomains( domain2 );
1351 probLaplace_->initializeLinSolverBuilder( problem_->getLinearSolverBuilder() );
1352 }
1353 }
1354
1355 if(laplaceInverse_.is_null()){// The Schwarz approximation of Laplace operator only needs to be build once
1356 if (verbose) {
1357 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1358 std::cout << "\t --- PCD: Setup A_p Schwarz Approximation " << std::endl;
1359 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1360 }
1361 BlockMatrixPtr_Type Ap = Teuchos::rcp( new BlockMatrix_Type(1) );
1362 Ap->addBlock(pressureLaplaceMatrixPtr_,0,0);
1363
1364 probLaplace_->initializeSystem( Ap );
1365
1366 probLaplace_->setupPreconditioner( "Monolithic" ); // single matrix
1367 laplaceInverse_ = probLaplace_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1368 }
1369
1370 if (probMass_.is_null()) {
1371 if (!timeProblem_.is_null()){
1372 probMass_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1373 DomainConstPtr_vec_Type domain2(0);
1374 domain2.push_back( timeProblem_->getDomain(1) );
1375 probMass_->initializeDomains( domain2 );
1376 probMass_->initializeLinSolverBuilder( timeProblem_->getLinearSolverBuilder() );
1377 }
1378 else{
1379 probMass_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1380 DomainConstPtr_vec_Type domain2(0);
1381 domain2.push_back( problem_->getDomain(1) );
1382 probMass_->initializeDomains( domain2 );
1383 probMass_->initializeLinSolverBuilder( problem_->getLinearSolverBuilder() );
1384 }
1385 }
1386
1387 if(massMatrixInverse_.is_null()){// The Schwarz approximation of pressure mass matrix only needs to be build once
1388 if (verbose) {
1389 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1390 std::cout << "\t --- PCD: Setup M_p Schwarz Approximation " << std::endl;
1391 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1392 }
1393 BlockMatrixPtr_Type Qp = Teuchos::rcp( new BlockMatrix_Type(1) );
1394 Qp->addBlock(pressureMassMatrixPtr_,0,0);
1395
1396 // Approximation of Mp is either done by Monolithic preconditioner or by a diagonal
1397 // Approximation with 'Diagonal' or TODO: AbsRowSum
1398 bool explicitInverse = parameterList->sublist("General").get("Mu Explicit Inverse",true);
1399 std::string typeDiag = parameterList->sublist("General").get("Diagonal Approximation","Diagonal");
1400
1401 if(explicitInverse)
1402 {
1403 probMass_->initializeSystem( Qp );
1404 probMass_->setupPreconditioner( "Monolithic" ); // single matrix
1405 massMatrixInverse_ = probMass_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1406 }
1407 else
1408 {
1409 massMatrixInverse_ = pressureMassMatrixPtr_->buildDiagonalInverse(typeDiag)->getThyraLinOpNonConst() ;
1410 }
1411 }
1412
1413 }
1414 else if (type == "LSC") {
1415
1416 if (verbose) {
1417 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1418 std::cout << "\t --- Building LSC Operator Components " << std::endl;
1419 std::cout << "\t --- -------------------------------------------------------- ---"<< std::endl;
1420 }
1421
1422 if (probLaplace_.is_null()) {
1423 if (!timeProblem_.is_null()){
1424 probLaplace_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1425 DomainConstPtr_vec_Type domain2(0);
1426 domain2.push_back( timeProblem_->getDomain(1) );
1427 probLaplace_->initializeDomains( domain2 );
1428 probLaplace_->initializeLinSolverBuilder( timeProblem_->getLinearSolverBuilder() );
1429 }
1430 else{
1431 probLaplace_ = Teuchos::rcp( new MinPrecProblem_Type( plSchur, comm ) );
1432 DomainConstPtr_vec_Type domain2(0);
1433 domain2.push_back( problem_->getDomain(1) );
1434 probLaplace_->initializeDomains( domain2 );
1435 probLaplace_->initializeLinSolverBuilder( problem_->getLinearSolverBuilder() );
1436 }
1437 }
1438
1439 BlockMatrixPtr_Type Ap = Teuchos::rcp( new BlockMatrix_Type(1) );
1440 Ap->addBlock(pressureLaplaceMatrixPtr_,0,0);
1441
1442 probLaplace_->initializeSystem( Ap );
1443
1444 probLaplace_->setupPreconditioner( "Monolithic" ); // single matrix
1445 laplaceInverse_ = probLaplace_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1446
1447 if (probVMass_.is_null()) {
1448 if (!timeProblem_.is_null()){
1449 probVMass_ = Teuchos::rcp( new MinPrecProblem_Type( plVelocity, comm ) );
1450 DomainConstPtr_vec_Type domain1(0);
1451 domain1.push_back( timeProblem_->getDomain(0) );
1452 probVMass_->initializeDomains( domain1 );
1453 probVMass_->initializeLinSolverBuilder( timeProblem_->getLinearSolverBuilder() );
1454 }
1455 else{
1456 probVMass_ = Teuchos::rcp( new MinPrecProblem_Type( plVelocity, comm ) );
1457 DomainConstPtr_vec_Type domain1(0);
1458 domain1.push_back( problem_->getDomain(0) );
1459 probVMass_->initializeDomains( domain1 );
1460 probVMass_->initializeLinSolverBuilder( problem_->getLinearSolverBuilder() );
1461 }
1462 }
1463 BlockMatrixPtr_Type Qv = Teuchos::rcp( new BlockMatrix_Type(1) );
1464 Qv->addBlock(velocityMassMatrixMatrixPtr_,0,0);
1465
1466 // Approximation of Mp is either done by Monolithic preconditioner or by a diagonal
1467 // Approximation with 'Diagonal' or TODO: AbsRowSum
1468 bool explicitInverse = parameterList->sublist("General").get("Mu Explicit Inverse",true);
1469 std::string typeDiag = parameterList->sublist("General").get("Diagonal Approximation","Diagonal");
1470
1471 if(explicitInverse)
1472 {
1473 probVMass_->initializeSystem( Qv );
1474 probVMass_->setupPreconditioner( "Monolithic" ); // single matrix
1475 massMatrixVInverse_ = probVMass_->getPreconditioner()->getThyraPrec()->getNonconstUnspecifiedPrecOp();
1476 }
1477 else
1478 {
1479 massMatrixVInverse_ = velocityMassMatrixMatrixPtr_->buildDiagonalInverse(typeDiag)->getThyraLinOpNonConst() ;
1480 }
1481
1482 }
1483 PRECONDITIONER_STOP(setupSInv);
1484
1485 // Building block Prec and passing along the different operators
1486 // that are required to build the respective preconditioners
1487 if (type == "Diagonal") {
1488 blockPrec2x2->setDiagonal(precVelocity_,
1489 precSchur_);
1490 }
1491 else if (type == "Triangular") {
1492 ThyraLinOpPtr_Type BT = system->getBlock(0,1)->getThyraLinOpNonConst();
1493 blockPrec2x2->setTriangular(precVelocity_,
1494 precSchur_,
1495 BT);
1496 }
1497 else if (type == "PCD") {
1498 if (!timeProblem_.is_null()){
1499 timeProblem_->assemble("UpdateConvectionDiffusionOperator");
1500 }
1501 else{
1502 problem_->assemble("UpdateConvectionDiffusionOperator");
1503 }
1504 MatrixPtr_Type pcdOperatorScaled = Teuchos::rcp( new Matrix_Type( pcdOperatorMatrixPtr_ ) );
1505 pcdOperatorScaled->resumeFill();
1506 pcdOperatorScaled->scale(-1.0);
1507 pcdOperatorScaled->fillComplete();
1508 ThyraLinOpPtr_Type BT = system->getBlock(0,1)->getThyraLinOpNonConst();
1509 blockPrec2x2->setTriangular(precVelocity_,
1510 laplaceInverse_,
1511 pcdOperatorScaled->getThyraLinOpNonConst(),
1512 massMatrixInverse_,
1513 massMatrixVInverse_,
1514 BT);
1515 ThyraLinOpPtr_Type B = system->getBlock(1,0)->getThyraLinOpNonConst();
1516 blockPrec2x2->setB(B);
1517 }
1518 else if (type == "LSC") {
1519 ThyraLinOpPtr_Type BT = system->getBlock(0,1)->getThyraLinOpNonConst();
1520 blockPrec2x2->setTriangular(precVelocity_,
1521 laplaceInverse_,
1522 massMatrixVInverse_,
1523 BT);
1524
1525
1526 ThyraLinOpPtr_Type B = system->getBlock(1,0)->getThyraLinOpNonConst();
1527 blockPrec2x2->setB(B);
1528
1529 ThyraLinOpPtr_Type F = system->getBlock(0,0)->getThyraLinOpNonConst();
1530 blockPrec2x2->setF(F);
1531
1532
1533 }
1534
1535 LinSolverBuilderPtr_Type solverBuilder;
1536 if (!timeProblem_.is_null())
1537 solverBuilder = timeProblem_->getLinearSolverBuilder();
1538 else
1539 solverBuilder = problem_->getLinearSolverBuilder();
1540
1541 if ( precFactory_.is_null() )
1542 precFactory_ = solverBuilder->createPreconditioningStrategy("");
1543
1544 if ( thyraPrec_.is_null() )
1545 thyraPrec_ = precFactory_->createPrec();
1546
1547 Teuchos::RCP< Thyra::DefaultPreconditioner<SC> > defaultPrec =
1548 Teuchos::rcp_dynamic_cast< Thyra::DefaultPreconditioner<SC> > (thyraPrec_);
1549 ThyraLinOpPtr_Type linOp =
1550 Teuchos::rcp_dynamic_cast< Thyra::LinearOpBase<SC> > (blockPrec2x2);
1551
1552 defaultPrec->initializeUnspecified( linOp );
1553
1554 precondtionerIsBuilt_ = true;
1555
1556 PRECONDITIONER_STOP(buildPreconditionerBlock2x2);
1557
1558
1559}
1560
1561
1562
1563template <class SC,class LO,class GO,class NO>
1564void Preconditioner<SC,LO,GO,NO>::setVelocityParameters( ParameterListPtr_Type parameterList, int coarseRanks )
1565{
1566 CommConstPtr_Type comm;
1567 if (!problem_.is_null())
1568 comm = problem_->getComm();
1569 else if(!timeProblem_.is_null())
1570 comm = timeProblem_->getComm();
1571
1572 bool verbose( comm->getRank() == 0 );
1573
1574 // Xpetra for now
1575 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::Map<LO,GO,NO> > > repeatedMaps(1);
1576 typedef Xpetra::Map<LO,GO,NO> XpetraMap_Type;
1577 typedef Teuchos::RCP<XpetraMap_Type> XpetraMapPtr_Type;
1578 typedef Teuchos::RCP<const XpetraMap_Type> XpetraMapConstPtr_Type;
1579 typedef const XpetraMapConstPtr_Type XpetraMapConstPtrConst_Type;
1580
1581 Teuchos::ArrayRCP<FROSch::DofOrdering> dofOrderings(1);
1582 Teuchos::ArrayRCP<UN> dofsPerNodeVector(1);
1583 ParameterListPtr_Type velocitySubList = sublist( sublist( sublist( sublist( parameterList, "Preconditioner Types" ) , "Teko" ) , "Inverse Factory Library" ) , "FROSch-Velocity" );
1584 dofsPerNodeVector[0] = (UN) velocitySubList->get( "DofsPerNode", 2);
1585 TEUCHOS_TEST_FOR_EXCEPTION(dofsPerNodeVector[0]<2, std::logic_error, "DofsPerNode for velocity must be atleast 2.");
1586
1587 /*Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp;
1588 if (!problem_.is_null())
1589 mapConstTmp = problem_->getDomain(0)->getMapVecFieldRepeated()->getTpetraMap();
1590 else if(!timeProblem_.is_null())
1591 mapConstTmp = timeProblem_->getDomain(0)->getMapVecFieldRepeated()->getTpetraMap();
1592
1593 Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);*/
1594
1595 MapConstPtr_Type mapConstTmp;
1596 if (!problem_.is_null())
1597 mapConstTmp = problem_->getDomain(0)->getMapVecFieldRepeated();
1598 else if(!timeProblem_.is_null())
1599 mapConstTmp = timeProblem_->getDomain(0)->getMapVecFieldRepeated();
1600
1601 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
1602 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
1603
1604 repeatedMaps[0] = mapX;
1605
1606 if (!velocitySubList->get( "DofOrdering", "NodeWise" ).compare("DimensionWise"))
1607 dofOrderings[0] = FROSch::DimensionWise;
1608 else if (!velocitySubList->get( "DofOrdering", "NodeWise" ).compare("NodeWise"))
1609 dofOrderings[0] = FROSch::NodeWise;
1610 else
1611 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Please chose a valid DOF ordering for ThyraPreconditioner(FROSch).");
1612
1613 velocitySubList->set("Repeated Map Vector",repeatedMaps);
1614 velocitySubList->set("DofOrdering Vector",dofOrderings);
1615 velocitySubList->set("DofsPerNode Vector",dofsPerNodeVector);
1616 velocitySubList->set( "Mpi Ranks Coarse", coarseRanks );
1617
1618 int lowerBound = 100000000;
1619 int upperBound = -1;
1620
1621 tuple_intint_Type rankRange;
1622 if (!problem_.is_null())
1623 rankRange = problem_->getDomain(0)->getMesh()->getRankRange();
1624 else if(!timeProblem_.is_null())
1625 rankRange = timeProblem_->getDomain(0)->getMesh()->getRankRange();
1626
1627 if (std::get<0>(rankRange) < lowerBound)
1628 lowerBound = std::get<0>(rankRange);
1629 if (std::get<1>(rankRange) > upperBound)
1630 upperBound = std::get<1>(rankRange);
1631
1632 int lowerBoundCoarse = lowerBound;
1633 int upperBoundCoarse = upperBound;
1634 // For now only for structured decompositions
1635 if ( coarseRanks > 0 ){
1636 lowerBoundCoarse = upperBound + 1;
1637 upperBoundCoarse = comm->getSize() - 1;
1638 }
1639 if (verbose) {
1640 std::cout << "\t --- ------------------------------------------------------------------- ---"<< std::endl;
1641 std::cout << "\t --- Range for local problems of preconditioner Teko velocity from " << lowerBound << " to " << upperBound << std::endl;
1642 std::cout << "\t --- Range for coarse problem of preconditioner Teko velocity from " << lowerBoundCoarse << " to " << upperBoundCoarse << std::endl;
1643 std::cout << "\t --- ------------------------------------------------------------------- ---"<< std::endl;
1644 }
1645
1646 velocitySubList->set( "Local problem ranks lower bound", lowerBound );
1647 velocitySubList->set( "Local problem ranks upper bound", upperBound );
1648
1649 velocitySubList->set( "Coarse problem ranks lower bound", lowerBoundCoarse );
1650 velocitySubList->set( "Coarse problem ranks upper bound", upperBoundCoarse );
1651
1652
1653
1654}
1655
1656
1657template <class SC,class LO,class GO,class NO>
1658void Preconditioner<SC,LO,GO,NO>::setPressureParameters( ParameterListPtr_Type parameterList, int coarseRanks )
1659{
1660 CommConstPtr_Type comm;
1661 if (!problem_.is_null())
1662 comm = problem_->getComm();
1663 else if(!timeProblem_.is_null())
1664 comm = timeProblem_->getComm();
1665
1666 bool verbose( comm->getRank() == 0 );
1667
1668 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::Map<LO,GO,NO> > > repeatedMaps(1);
1669 typedef Xpetra::Map<LO,GO,NO> XpetraMap_Type;
1670 typedef Teuchos::RCP<XpetraMap_Type> XpetraMapPtr_Type;
1671 typedef Teuchos::RCP<const XpetraMap_Type> XpetraMapConstPtr_Type;
1672 typedef const XpetraMapConstPtr_Type XpetraMapConstPtrConst_Type;
1673
1674 Teuchos::ArrayRCP<FROSch::DofOrdering> dofOrderings(1);
1675 Teuchos::ArrayRCP<UN> dofsPerNodeVector(1);
1676 ParameterListPtr_Type pressureSubList = sublist( sublist( sublist( sublist( parameterList, "Preconditioner Types" ) , "Teko" ) , "Inverse Factory Library" ) , "FROSch-Pressure" );
1677 dofsPerNodeVector[0] = (UN) pressureSubList->get( "DofsPerNode", 1);
1678 TEUCHOS_TEST_FOR_EXCEPTION(dofsPerNodeVector[0]!=1, std::logic_error, "DofsPerNode for pressure must be 1.");
1679
1680 /*Teuchos::RCP<const Tpetra::Map<LO,GO,NO> > mapConstTmp;
1681
1682 if (!problem_.is_null()){
1683 if ( problem_->getDomain(1)->getFEType()=="P0" )
1684 mapConstTmp = problem_->getDomain(1)->getElementMap()->getTpetraMap();
1685 else
1686 mapConstTmp = problem_->getDomain(1)->getMapRepeated()->getTpetraMap();
1687 }
1688 else if(!timeProblem_.is_null()){
1689 if ( timeProblem_->getDomain(1)->getFEType()=="P0" )
1690 mapConstTmp = timeProblem_->getDomain(1)->getElementMap()->getTpetraMap();
1691 else
1692 mapConstTmp = timeProblem_->getDomain(1)->getMapRepeated()->getTpetraMap();
1693 }*/
1694 MapConstPtr_Type mapConstTmp;
1695 if (!problem_.is_null()){
1696 if ( problem_->getDomain(1)->getFEType()=="P0" )
1697 mapConstTmp = problem_->getDomain(1)->getElementMap();
1698 else
1699 mapConstTmp = problem_->getDomain(1)->getMapRepeated();
1700 }
1701 else if(!timeProblem_.is_null()){
1702 if ( timeProblem_->getDomain(1)->getFEType()=="P0" )
1703 mapConstTmp = timeProblem_->getDomain(1)->getElementMap();
1704 else
1705 mapConstTmp = timeProblem_->getDomain(1)->getMapRepeated();
1706 }
1707
1708 //Teuchos::RCP<Tpetra::Map<LO,GO,NO> > mapTmp = Teuchos::rcp_const_cast<Tpetra::Map<LO,GO,NO> > (mapConstTmp);
1709
1710 XpetraMapConstPtr_Type mapConstX = Xpetra::MapFactory<LO,GO,NO>::Build( Xpetra::UseTpetra, mapConstTmp->getGlobalNumElements(), mapConstTmp->getNodeElementList(), mapConstTmp->getIndexBase(), mapConstTmp->getComm() );
1711 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapX= Teuchos::rcp_const_cast<Xpetra::Map<LO,GO,NO> > (mapConstX);
1712
1713 repeatedMaps[0] = mapX;
1714
1715 if (!pressureSubList->get( "DofOrdering", "NodeWise" ).compare("DimensionWise"))
1716 dofOrderings[0] = FROSch::DimensionWise;
1717 else if (!pressureSubList->get( "DofOrdering", "NodeWise" ).compare("NodeWise"))
1718 dofOrderings[0] = FROSch::NodeWise;
1719 else
1720 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Please chose a valid DOF ordering for ThyraPreconditioner(FROSch).");
1721
1722 pressureSubList->set("Repeated Map Vector",repeatedMaps);
1723 pressureSubList->set("DofOrdering Vector",dofOrderings);
1724 pressureSubList->set("DofsPerNode Vector",dofsPerNodeVector);
1725 pressureSubList->set( "Mpi Ranks Coarse", coarseRanks );
1726
1727 int lowerBound = 100000000;
1728 int upperBound = -1;
1729
1730
1731 tuple_intint_Type rankRange;
1732 if (!problem_.is_null())
1733 rankRange = problem_->getDomain(1)->getMesh()->getRankRange();
1734 else if(!timeProblem_.is_null())
1735 rankRange = timeProblem_->getDomain(1)->getMesh()->getRankRange();
1736
1737 if (std::get<0>(rankRange) < lowerBound)
1738 lowerBound = std::get<0>(rankRange);
1739 if (std::get<1>(rankRange) > upperBound)
1740 upperBound = std::get<1>(rankRange);
1741
1742 int lowerBoundCoarse = lowerBound;
1743 int upperBoundCoarse = upperBound;
1744 // For now only for structured decompositions
1745 if ( coarseRanks > 0 ){
1746 lowerBoundCoarse = upperBound + 1;
1747 upperBoundCoarse = comm->getSize() - 1;
1748 }
1749 if (verbose) {
1750 std::cout << "\t --- ------------------------------------------------------------------- ---"<< std::endl;
1751 std::cout << "\t --- Range for local problems of preconditioner Teko pressure from " << lowerBound << " to " << upperBound << std::endl;
1752 std::cout << "\t --- Range for coarse problem of preconditioner Teko pressure from " << lowerBoundCoarse << " to " << upperBoundCoarse << std::endl;
1753 std::cout << "\t --- ------------------------------------------------------------------- ---"<< std::endl;
1754 }
1755
1756 pressureSubList->set( "Local problem ranks lower bound", lowerBound );
1757 pressureSubList->set( "Local problem ranks upper bound", upperBound );
1758
1759 pressureSubList->set( "Coarse problem ranks lower bound", lowerBoundCoarse );
1760 pressureSubList->set( "Coarse problem ranks upper bound", upperBoundCoarse );
1761
1762
1763}
1764
1765template <class SC,class LO,class GO,class NO>
1766typename Preconditioner<SC,LO,GO,NO>::ThyraLinOpConstPtr_Type Preconditioner<SC,LO,GO,NO>::getTekoOp(){
1767 return tekoLinOp_;
1768}
1769
1770template <class SC,class LO,class GO,class NO>
1771void Preconditioner<SC,LO,GO,NO>::setVelocityMassMatrix(MatrixPtr_Type massMatrix) const{
1772 velocityMassMatrix_ = massMatrix->getThyraLinOp();
1773 velocityMassMatrixMatrixPtr_ = massMatrix;
1774}
1775#endif
1776
1777template <class SC,class LO,class GO,class NO>
1778void Preconditioner<SC,LO,GO,NO>::exportCoarseBasis( ){
1779
1780 CommConstPtr_Type comm;
1781 if (!problem_.is_null())
1782 comm = problem_->getComm();
1783 else if(!timeProblem_.is_null())
1784 comm = timeProblem_->getComm();
1785
1786 TEUCHOS_TEST_FOR_EXCEPTION( pListPhiExport_.is_null(), std::runtime_error, "No parameterlist to extract Phi pointer.");
1787 ParameterListPtr_Type pLPrec = sublist( sublist( pListPhiExport_, "Preconditioner Types" ) ,"FROSch" );
1788 std::string coarseType = pLPrec->get( "CoarseOperator Type", "RGDSWCoarseOperator" );
1789 ParameterListPtr_Type pLCoarse = sublist( pLPrec, coarseType );
1790// TEUCHOS_TEST_FOR_EXCEPTION( coarseType!="RGDSWCoarseOperator" && coarseType!="GDSWCoarseOperator", std::runtime_error, "Export Phi only for GDSWCoarseOperator and RGDSWCoarseOperator.");
1791
1792 TEUCHOS_TEST_FOR_EXCEPTION( !pLCoarse->isParameter("RCP(Phi)"), std::runtime_error, "No parameter to extract Phi pointer.");
1793
1794 Teuchos::RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > phiTpetra;
1795
1796 TEUCHOS_TEST_FOR_EXCEPTION( !pLCoarse->isType<decltype(phiTpetra)>("RCP(Phi)"), std::runtime_error, "Wrong type of pointer to extract Phi.");
1797
1798 phiTpetra = pLCoarse->get<decltype(phiTpetra)>("RCP(Phi)");
1799
1800 MatrixPtr_Type phiMatrix = Teuchos::rcp( new Matrix_Type( phiTpetra ) );
1801 int numberOfBlocks;
1802 {
1803 ParameterListPtr_Type parameterList;
1804 if (!problem_.is_null())
1805 parameterList = problem_->getParameterList();
1806 else if(!timeProblem_.is_null())
1807 parameterList = timeProblem_->getParameterList();
1808 numberOfBlocks = parameterList->get("Number of blocks",1);
1809 }
1810
1811 std::vector<MapConstPtr_Type> blockMaps(numberOfBlocks);
1812 MultiVectorPtr_Type phiMV;
1813 phiMatrix->toMV( phiMV );
1814
1815 for (UN i = 0; i < numberOfBlocks; i++) {
1816 int dofsPerNode = pLPrec->get( "DofsPerNode" + std::to_string(i+1), 1);
1817 MapConstPtr_Type map;
1818 if (dofsPerNode>1) {
1819 if (!problem_.is_null()) {
1820 TEUCHOS_TEST_FOR_EXCEPTION( problem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1821 map = problem_->getDomain(i)->getMapVecFieldUnique();
1822 }
1823 else if(!timeProblem_.is_null()){
1824 TEUCHOS_TEST_FOR_EXCEPTION( timeProblem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1825 map = timeProblem_->getDomain(i)->getMapVecFieldUnique();
1826 }
1827 }
1828 else{
1829 if (!problem_.is_null()) {
1830 TEUCHOS_TEST_FOR_EXCEPTION( problem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1831 map = problem_->getDomain(i)->getMapUnique();
1832 }
1833 else if(!timeProblem_.is_null()){
1834 TEUCHOS_TEST_FOR_EXCEPTION( timeProblem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1835 map = timeProblem_->getDomain(i)->getMapUnique();
1836 }
1837 }
1838 blockMaps[i] = map;
1839 }
1840
1841
1842 BlockMultiVectorPtr_Type phiBlockMV = Teuchos::rcp( new BlockMultiVector_Type ( blockMaps, phiMV->getNumVectors() ) );
1843 phiBlockMV->setMergedVector( phiMV );
1844 phiBlockMV->split();
1845
1846 ParameterListPtr_Type pLProblem;
1847 if (!problem_.is_null())
1848 pLProblem = problem_->getParameterList();
1849 else if(!timeProblem_.is_null())
1850 pLProblem = timeProblem_->getParameterList();
1851
1852 BlockMultiVectorPtr_Type phiSumBlockMV = phiBlockMV->sumColumns();
1853
1854 for (int i=0; i<phiBlockMV->size(); i++) {
1855 bool exportThisBlock = !pLProblem->sublist("Exporter").get( "Exclude coarse functions block" + std::to_string(i+1), false);
1856
1857 if (exportThisBlock){
1858 Teuchos::RCP<ExporterParaView<SC,LO,GO,NO> > exporter(new ExporterParaView<SC,LO,GO,NO>());
1859
1860 DomainConstPtr_Type domain;
1861 if (!problem_.is_null())
1862 domain = problem_->getDomain(i);
1863 else if(!timeProblem_.is_null())
1864 domain = timeProblem_->getDomain(i);
1865
1866 int dofsPerNode = domain->getDimension();
1867 std::string fileName = pLProblem->sublist("Exporter").get( "Name coarse functions block" + std::to_string(i+1), "phi");
1868
1869 MeshPtr_Type meshNonConst = Teuchos::rcp_const_cast<Mesh_Type>( domain->getMesh() );
1870 exporter->setup(fileName, meshNonConst, domain->getFEType());
1871
1872 for (int j=0; j<phiBlockMV->getNumVectors(); j++) {
1873
1874 MultiVectorConstPtr_Type exportPhi = phiBlockMV->getBlock( i )->getVector( j );
1875
1876 std::string varName = fileName + std::to_string(j);
1877 if ( pLPrec->get( "DofsPerNode" + std::to_string(i+1), 1) > 1 )
1878 exporter->addVariable( exportPhi, varName, "Vector", dofsPerNode, domain->getMapUnique() );
1879 else
1880 exporter->addVariable( exportPhi, varName, "Scalar", 1, domain->getMapUnique() );
1881
1882 }
1883 MultiVectorConstPtr_Type exportSumPhi = phiSumBlockMV->getBlock(i);
1884 std::string varName = fileName + "Sum";
1885 if ( pLPrec->get( "DofsPerNode" + std::to_string(i+1), 1) > 1 )
1886 exporter->addVariable( exportSumPhi, varName, "Vector", dofsPerNode, domain->getMapUnique() );
1887 else
1888 exporter->addVariable( exportSumPhi, varName, "Scalar", 1, domain->getMapUnique() );
1889
1890 exporter->save(0.0);
1891 }
1892 }
1893
1894}
1895
1896template <class SC,class LO,class GO,class NO>
1897void Preconditioner<SC,LO,GO,NO>::exportCoarseBasisFSI( ){
1898
1899 CommConstPtr_Type comm;
1900 if (!problem_.is_null())
1901 comm = problem_->getComm();
1902 else if(!timeProblem_.is_null())
1903 comm = timeProblem_->getComm();
1904
1905 TEUCHOS_TEST_FOR_EXCEPTION( pListPhiExport_.is_null(), std::runtime_error, "No parameterlist to extract Phi pointer.");
1906 ParameterListPtr_Type pLPrec = sublist( sublist( pListPhiExport_, "Preconditioner Types" ) ,"FROSch" );
1907 std::string coarseType = pLPrec->get( "CoarseOperator Type", "RGDSWCoarseOperator" );
1908 ParameterListPtr_Type pLCoarse = sublist( pLPrec, coarseType );
1909// TEUCHOS_TEST_FOR_EXCEPTION( coarseType!="RGDSWCoarseOperator" && coarseType!="GDSWCoarseOperator", std::runtime_error, "Export Phi only for GDSWCoarseOperator and RGDSWCoarseOperator.");
1910
1911 TEUCHOS_TEST_FOR_EXCEPTION( !pLCoarse->isParameter("Phi Pointer"), std::runtime_error, "No parameter to extract Phi pointer.");
1912
1913 Teuchos::RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > phiTpetra;
1914
1915 TEUCHOS_TEST_FOR_EXCEPTION( !pLCoarse->isType<decltype(phiTpetra)>("Phi Pointer"), std::runtime_error, "Wrong type of pointer to extract Phi.");
1916
1917 phiTpetra = pLCoarse->get<decltype(phiTpetra)>("Phi Pointer");
1918
1919 MatrixPtr_Type phiMatrix = Teuchos::rcp( new Matrix_Type( phiTpetra ) );
1920 int numberOfBlocks;
1921 {
1922 ParameterListPtr_Type parameterList;
1923 if (!problem_.is_null())
1924 parameterList = problem_->getParameterList();
1925 else if(!timeProblem_.is_null())
1926 parameterList = timeProblem_->getParameterList();
1927 numberOfBlocks = parameterList->get("Number of blocks",1);
1928 }
1929
1930 std::vector<MapConstPtr_Type> blockMaps(numberOfBlocks);
1931 MultiVectorPtr_Type phiMV;
1932 phiMatrix->toMV( phiMV );
1933
1934 for (UN i = 0; i < numberOfBlocks; i++) {
1935 int dofsPerNode = pLPrec->get( "DofsPerNode" + std::to_string(i+1), 1);
1936 MapConstPtr_Type map;
1937 if (i==3) {
1938 map = timeProblem_->getDomain(i)->getInterfaceMapVecFieldUnique();
1939 }
1940 else {
1941 if (dofsPerNode>1) {
1942 if (!problem_.is_null()) {
1943 TEUCHOS_TEST_FOR_EXCEPTION( problem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1944 map = problem_->getDomain(i)->getMapVecFieldUnique();
1945 }
1946 else if(!timeProblem_.is_null()){
1947 TEUCHOS_TEST_FOR_EXCEPTION( timeProblem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1948 map = timeProblem_->getDomain(i)->getMapVecFieldUnique();
1949 }
1950 }
1951 else{
1952 if (!problem_.is_null()) {
1953 TEUCHOS_TEST_FOR_EXCEPTION( problem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1954 map = problem_->getDomain(i)->getMapUnique();
1955 }
1956 else if(!timeProblem_.is_null()){
1957 TEUCHOS_TEST_FOR_EXCEPTION( timeProblem_->getDomain(i)->getFEType() == "P0", std::logic_error, "Vector field map not implemented for P0 elements.");
1958 map = timeProblem_->getDomain(i)->getMapUnique();
1959 }
1960 }
1961 }
1962 blockMaps[i] = map;
1963 }
1964
1965
1966 BlockMultiVectorPtr_Type phiBlockMV = Teuchos::rcp( new BlockMultiVector_Type ( blockMaps, phiMV->getNumVectors() ) );
1967 phiBlockMV->setMergedVector( phiMV );
1968 phiBlockMV->split();
1969
1970 ParameterListPtr_Type pLProblem;
1971 if (!problem_.is_null())
1972 pLProblem = problem_->getParameterList();
1973 else if(!timeProblem_.is_null())
1974 pLProblem = timeProblem_->getParameterList();
1975
1976 BlockMultiVectorPtr_Type phiSumBlockMV = phiBlockMV->sumColumns();
1977
1978 for (int i=0; i<phiBlockMV->size(); i++) {
1979 bool exportThisBlock = !pLProblem->sublist("Exporter").get( "Exclude coarse functions block" + std::to_string(i+1), false);
1980
1981 if (exportThisBlock){
1982 Teuchos::RCP<ExporterParaView<SC,LO,GO,NO> > exporter(new ExporterParaView<SC,LO,GO,NO>());
1983
1984 DomainConstPtr_Type domain;
1985 if (!problem_.is_null())
1986 domain = problem_->getDomain(i);
1987 else if(!timeProblem_.is_null())
1988 domain = timeProblem_->getDomain(i);
1989
1990 int dofsPerNode = domain->getDimension();
1991 std::string fileName = pLProblem->sublist("Exporter").get( "Name coarse functions block" + std::to_string(i+1), "phi");
1992
1993 MeshPtr_Type meshNonConst = Teuchos::rcp_const_cast<Mesh_Type>( domain->getMesh() );
1994 exporter->setup(fileName, meshNonConst, domain->getFEType());
1995
1996 for (int j=0; j<phiBlockMV->getNumVectors(); j++) {
1997
1998 MultiVectorConstPtr_Type exportPhi = phiBlockMV->getBlock( i )->getVector( j );
1999
2000 std::string varName = fileName + std::to_string(j);
2001 if ( pLPrec->get( "DofsPerNode" + std::to_string(i+1), 1) > 1 )
2002 exporter->addVariable( exportPhi, varName, "Vector", dofsPerNode, domain->getMapUnique() );
2003 else
2004 exporter->addVariable( exportPhi, varName, "Scalar", 1, domain->getMapUnique() );
2005
2006 }
2007 MultiVectorConstPtr_Type exportSumPhi = phiSumBlockMV->getBlock(i);
2008 std::string varName = fileName + "Sum";
2009 if ( pLPrec->get( "DofsPerNode" + std::to_string(i+1), 1) > 1 )
2010 exporter->addVariable( exportSumPhi, varName, "Vector", dofsPerNode, domain->getMapUnique() );
2011 else
2012 exporter->addVariable( exportSumPhi, varName, "Scalar", 1, domain->getMapUnique() );
2013
2014 exporter->save(0.0);
2015 }
2016 }
2017
2018}
2019
2020
2021}
2022
2023#endif
FEDD::Preconditioner::setPressureProjection
void setPressureProjection(BlockMultiVectorPtr_Type pressureProjection) const
Setting pressure projection from specific prolem to be used in preconditioner.
Definition Preconditioner_def.hpp:175
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36