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