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