Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
PrecOpFaCSI_def.hpp
1#ifndef PrecOpFaCSI_DEF_hpp
2#define PrecOpFaCSI_DEF_hpp
3#include "PrecOpFaCSI_decl.hpp"
4#include <Thyra_TpetraMultiVector_decl.hpp>
5#include <Teuchos_VerboseObject.hpp>
14
15namespace FEDD {
16using namespace Thyra;
17
18// Constructors
19
20template<class SC, class LO, class GO, class NO>
21PrecOpFaCSI<SC,LO,GO,NO>::PrecOpFaCSI()
22:PreconditionerOperator<SC,LO,GO,NO>(),
23fluidPrecMonolithic_(false),
24useSolidPreconditioner_(true)
25{
26// TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Still a problem for FaCSI, cant be used yet.");
27}
28
29template<class SC, class LO, class GO, class NO>
30PrecOpFaCSI<SC,LO,GO,NO>::PrecOpFaCSI(CommConstPtr_Type comm, bool fluidPrecMonolithic, bool useFluidPreconditioner, bool useSolidPreconditioner, bool onlyDiagonal)
31:PreconditionerOperator<SC,LO,GO,NO>(),
32fluidPrecMonolithic_(fluidPrecMonolithic),
33useFluidPreconditioner_(useFluidPreconditioner),
34useSolidPreconditioner_(useSolidPreconditioner),
35onlyDiagonal_(onlyDiagonal)
36{
37 comm_=comm;
38// TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Still a problem for FaCSI, cant be used yet.");
39}
40
41template<class SC, class LO, class GO, class NO>
42void PrecOpFaCSI<SC,LO,GO,NO>::setGE(ThyraLinOpPtr_Type C1,
43 ThyraLinOpPtr_Type C1T,
44 ThyraLinOpPtr_Type C2,
45 ThyraLinOpPtr_Type sInv,
46 ThyraLinOpPtr_Type fInv,
47 ThyraLinOpPtr_Type fF,
48 ThyraLinOpPtr_Type fBT){
49
50 setC1(C1);
51 setC1T(C1T);
52 setC2(C2);
53 setStructInv(sInv);
54 setFluidInv(fInv);
55 setFluidF(fF);
56 setFluidBT(fBT);
57
58 initialize();
59}
60
61template<class SC, class LO, class GO, class NO>
62void PrecOpFaCSI<SC,LO,GO,NO>::setGI(ThyraLinOpPtr_Type C1,
63 ThyraLinOpPtr_Type C1T,
64 ThyraLinOpPtr_Type C2,
65 ThyraLinOpPtr_Type C4,
66 ThyraLinOpPtr_Type sInv,
67 ThyraLinOpPtr_Type fInv,
68 ThyraLinOpPtr_Type fF,
69 ThyraLinOpPtr_Type fBT,
70 ThyraLinOpPtr_Type gInv){
71
72 setC1(C1);
73 setC1T(C1T);
74 setC2(C2);
75 setC4(C4);
76 setStructInv(sInv);
77 setFluidInv(fInv);
78 setFluidF(fF);
79 setFluidBT(fBT);
80 setGeoInv(gInv);
81
82 initialize();
83}
84
85template<class SC, class LO, class GO, class NO>
86void PrecOpFaCSI<SC,LO,GO,NO>::setGIShape(ThyraLinOpPtr_Type C1,
87 ThyraLinOpPtr_Type C1T,
88 ThyraLinOpPtr_Type C2,
89 ThyraLinOpPtr_Type C4,
90 ThyraLinOpPtr_Type sInv,
91 ThyraLinOpPtr_Type fInv,
92 ThyraLinOpPtr_Type fF,
93 ThyraLinOpPtr_Type fBT,
94 ThyraLinOpPtr_Type gInv,
95 ThyraLinOpPtr_Type shape_v,
96 ThyraLinOpPtr_Type shape_p){
97
98 setC1(C1);
99 setC1T(C1T);
100 setC2(C2);
101 setC4(C4);
102 setStructInv(sInv);
103 setFluidInv(fInv);
104 setFluidF(fF);
105 setFluidBT(fBT);
106 setGeoInv(gInv);
107 setShapeDeriv(shape_v, shape_p);
108
109 initialize();
110}
111
112
113template<class SC, class LO, class GO, class NO>
114void PrecOpFaCSI<SC,LO,GO,NO>::setC1(ThyraLinOpPtr_Type C1){
115 C1_ = C1;
116}
117template<class SC, class LO, class GO, class NO>
118void PrecOpFaCSI<SC,LO,GO,NO>::setC1T(ThyraLinOpPtr_Type C1T){
119 C1T_ = C1T;
120}
121template<class SC, class LO, class GO, class NO>
122void PrecOpFaCSI<SC,LO,GO,NO>::setC2(ThyraLinOpPtr_Type C2){
123 C2_ = C2;
124}
125template<class SC, class LO, class GO, class NO>
126void PrecOpFaCSI<SC,LO,GO,NO>::setC4(ThyraLinOpPtr_Type C4){
127 C4_ = C4;
128}
129template<class SC, class LO, class GO, class NO>
130void PrecOpFaCSI<SC,LO,GO,NO>::setShapeDeriv(ThyraLinOpPtr_Type shape_v, ThyraLinOpPtr_Type shape_p){
131 shape_v_ = shape_v;
132 shape_p_ = shape_p;
133}
134template<class SC, class LO, class GO, class NO>
135void PrecOpFaCSI<SC,LO,GO,NO>::setStructInv(ThyraLinOpPtr_Type sInv){
136 sInv_ = sInv;
137}
138template<class SC, class LO, class GO, class NO>
139void PrecOpFaCSI<SC,LO,GO,NO>::setFluidInv(ThyraLinOpPtr_Type fInv){
140 fInv_ = fInv;
141}
142template<class SC, class LO, class GO, class NO>
143void PrecOpFaCSI<SC,LO,GO,NO>::setGeoInv(ThyraLinOpPtr_Type gInv){
144 gInv_ = gInv;
145}
146template<class SC, class LO, class GO, class NO>
147void PrecOpFaCSI<SC,LO,GO,NO>::setFluidF(ThyraLinOpPtr_Type fF){
148 fF_ = fF;
149}
150template<class SC, class LO, class GO, class NO>
151void PrecOpFaCSI<SC,LO,GO,NO>::setFluidBT(ThyraLinOpPtr_Type fBT){
152 fBT_ = fBT;
153}
154
155template<class SC, class LO, class GO, class NO>
156void PrecOpFaCSI<SC,LO,GO,NO>::initialize(){
157 TEUCHOS_TEST_FOR_EXCEPTION(fInv_.is_null(), std::runtime_error,"Can not initialize FaCSI preconditioner: Fluid preconditioner not set.");
158 TEUCHOS_TEST_FOR_EXCEPTION(sInv_.is_null(), std::runtime_error,"Can not initialize FaCSI preconditioner: Structure preconditioner not set.");
159 TEUCHOS_TEST_FOR_EXCEPTION(C1_.is_null(), std::runtime_error,"Can not initialize FaCSI preconditioner: C1 not set.");
160 Teuchos::Array< Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesRange( 4 );
161 Teuchos::Array< Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesDomain( 4 );
162 vectorSpacesRange[0] = fF_->range();
163 vectorSpacesRange[1] = fBT_->domain();
164 vectorSpacesRange[2] = sInv_->range();
165 vectorSpacesRange[3] = C1_->range();
166
167 vectorSpacesDomain[0] = fF_->domain();
168 vectorSpacesDomain[1] = fBT_->domain();
169 vectorSpacesDomain[2] = sInv_->domain();
170 vectorSpacesDomain[3] = C1T_->domain();
171 if ( !gInv_.is_null() ) {
172 vectorSpacesRange.push_back( gInv_->range() );
173 vectorSpacesDomain.push_back( gInv_->domain() );
174 }
175 // defaultProductRange_;
176 // Teuchos::RCP<const Thyra::DefaultProductVectorSpace<SC> > defaultProductDomain_;
177
178 Teuchos::RCP<const Thyra::DefaultProductVectorSpace<SC> > pR = Thyra::productVectorSpace<SC>( vectorSpacesRange );
179 Teuchos::RCP<const Thyra::DefaultProductVectorSpace<SC> > pD = Thyra::productVectorSpace<SC>( vectorSpacesDomain );
180// Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > pVSR = Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<SC> >(pR);
181// Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > pVSD = Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<SC> >(pD);
182//
183// this->defaultProductRange_ = Thyra::multiVectorProductVectorSpace<SC>( pVSR , 1);
184// this->defaultProductDomain_ = Thyra::multiVectorProductVectorSpace<SC>( pVSD, 1);
185
186 this->defaultProductRange_ = pR;
187 this->defaultProductDomain_ = pD;
188}
189
190template<class SC, class LO, class GO, class NO>
191void PrecOpFaCSI<SC,LO,GO,NO>::applyIt(
192 const EOpTransp M_trans,
193 const MultiVectorBase<SC> &X_in,
194 const Ptr<MultiVectorBase<SC> > &Y_inout,
195 const SC alpha,
196 const SC beta
197 ) const
198{
199 applyImpl(M_trans, X_in, Y_inout, alpha, beta);
200
201}
202
203template<class SC, class LO, class GO, class NO>
204void PrecOpFaCSI<SC,LO,GO,NO>::applyImpl(
205 const EOpTransp M_trans,
206 const MultiVectorBase<SC> &X_in,
207 const Ptr<MultiVectorBase<SC> > &Y_inout,
208 const SC alpha,
209 const SC beta
210 ) const
211{
212 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
213
214 using Teuchos::rcpFromRef;
215 typedef Teuchos::ScalarTraits<SC> ST;
216 typedef RCP<MultiVectorBase<SC> > MultiVectorPtr;
217 typedef RCP<const MultiVectorBase<SC> > ConstMultiVectorPtr;
218 typedef RCP<const LinearOpBase<SC> > ConstLinearOpPtr;
219
220 int rank = comm_->getRank();
221 if (onlyDiagonal_) {
222
223 Teuchos::RCP<const Thyra::ProductMultiVectorBase<SC> > X
224 = Teuchos::rcp_dynamic_cast<const Thyra::ProductMultiVectorBase<SC> > ( rcpFromRef(X_in) );
225
226 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > Y
227 = Teuchos::rcp_dynamic_cast< Thyra::ProductMultiVectorBase<SC> > ( rcpFromPtr(Y_inout) );
228
229 Y_inout->assign(0.);
230
231 Teuchos::RCP< const MultiVectorBase< SC > > X_s = X->getMultiVectorBlock(2);
232 Teuchos::RCP< MultiVectorBase< SC > > Y_s = Y->getNonconstMultiVectorBlock(2);
233
234
235 Teuchos::RCP< const MultiVectorBase< SC > > X_fv = X->getMultiVectorBlock(0);
236 Teuchos::RCP< MultiVectorBase< SC > > Y_fv = Y->getNonconstMultiVectorBlock(0);
237 Teuchos::RCP< const MultiVectorBase< SC > > X_fp = X->getMultiVectorBlock(1);
238 Teuchos::RCP< MultiVectorBase< SC > > Y_fp = Y->getNonconstMultiVectorBlock(1);
239
240
241 Teuchos::RCP< MultiVectorBase< SC > > Y_l = Y->getNonconstMultiVectorBlock(3);
242 Teuchos::RCP<const MultiVectorBase< SC > > X_l = X->getMultiVectorBlock(3);
243
244 assign(Y_fv.ptr(), *X_fv);
245 assign(Y_fp.ptr(), *X_fp);
246 assign(Y_s.ptr(), *X_s);
247 assign(Y_l.ptr(), *X_l);
248
249// Teuchos::RCP< Thyra::TpetraMultiVector< SC, LO, GO, NO > > YsTpetra =
250// Teuchos::rcp_dynamic_cast< Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( Y_s );
251//
252//
257// Teuchos::RCP< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > XsTpetra =
258// Teuchos::rcp_dynamic_cast< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( X_s );
259// Teuchos::RCP< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > XfvT =
260// Teuchos::rcp_dynamic_cast< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( X_fv );
261// Teuchos::RCP< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > XfpT =
262// Teuchos::rcp_dynamic_cast< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( X_fp );
263// Teuchos::RCP< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > XlT =
264// Teuchos::rcp_dynamic_cast< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( X_l );
265// std::cout << "Xs:" << std::endl;
266// XsTpetra->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
267// comm_->barrier(); comm_->barrier(); comm_->barrier();
268// std::cout << "Xfv:" << std::endl;
269// XfvT->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
270// comm_->barrier(); comm_->barrier(); comm_->barrier();
271// std::cout << "Xfp:" << std::endl;
272// XfpT->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
273// comm_->barrier(); comm_->barrier(); comm_->barrier();
274// std::cout << "Xl:" << std::endl;
275// XlT->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
276// comm_->barrier(); comm_->barrier(); comm_->barrier();
277//
278//
279//
280// if (useSolidPreconditioner_)
281// sInv_->apply(NOTRANS, *X_s, Y_s.ptr(), 1., 0.);
282// else
283// assign(Y_s.ptr(), *X_s);
284//
285//
289//
290//
291// if (productRangeFluid_.is_null()) {
292// Teuchos::Array< Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesRangeFluid( 2 );
293//
294// vectorSpacesRangeFluid[0] = X_fv->range();
295// vectorSpacesRangeFluid[1] = X_fp->range();
296//
297// productRangeFluid_ = Thyra::productVectorSpace<SC>( vectorSpacesRangeFluid );
298// }
299//
300// Teuchos::Array< Teuchos::RCP< Thyra::MultiVectorBase< SC > > > X_fluid( 2 );
301// Teuchos::Array< Teuchos::RCP< Thyra::MultiVectorBase< SC > > > Y_fluid( 2 );
302//
303// X_fluid[0] = Y_fv;
304// X_fluid[1] = Y_fp;
305//
306// Y_fluid[0] = Y_fv;
307// Y_fluid[1] = Y_fp;
308//
309//
310//
311// // std::cout << "Y_fp before mono" << std::endl;
312// // Y_fp->describe(*out,Teuchos::VERB_EXTREME);
313// // comm_->barrier(); comm_->barrier(); comm_->barrier();
314//
315// if (useFluidPreconditioner_){
316//
317// if (fluidPrecMonolithic_) {
318// Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > fMonoVS = fInv_->domain();
319//
320// if ( X_fmono_.is_null() ){
321// X_fmono_ = createMembers( fMonoVS, X_fluid[0]->domain()->dim() );
322// Y_fmono_ = createMembers( fMonoVS, Y_fluid[0]->domain()->dim() );
323// }
324// //We can/should speedup this process
325// copyToMono(X_fluid);
326// fInv_->apply(NOTRANS, *X_fmono_, Y_fmono_.ptr(), 1., 0.);
327// copyFromMono(Y_fluid);
328// }
329// else{
330// Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > prodX_f = Thyra::defaultProductMultiVector<SC>( productRangeFluid_, X_fluid );
331// Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > prodY_f = Thyra::defaultProductMultiVector<SC>( productRangeFluid_, Y_fluid );
332//
333// Teuchos::RCP< MultiVectorBase<SC> > X_f = Teuchos::rcp_dynamic_cast<MultiVectorBase< SC > >(prodX_f);
334// Teuchos::RCP< MultiVectorBase<SC> > Y_f = Teuchos::rcp_dynamic_cast<MultiVectorBase< SC > >(prodY_f);
335//
336// fInv_->apply(NOTRANS, *X_f, Y_f.ptr(), 1., 0.);
337// }
338// }
339// else{
340// assign(Y_fv.ptr(), *X_fv);
341// assign(Y_fp.ptr(), *X_fp);
342// }
343//
344// Teuchos::RCP< Thyra::TpetraMultiVector< SC, LO, GO, NO > > YfvTpetra =
345// Teuchos::rcp_dynamic_cast< Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( Y_fv );
346// Teuchos::RCP< Thyra::TpetraMultiVector< SC, LO, GO, NO > > YfpTpetra =
347// Teuchos::rcp_dynamic_cast< Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( Y_fp );
348//
356//
357//
358// assign(Y_l.ptr(), *X_l);
359//
362//
363//
365
366 }
367 else{
368 Teuchos::RCP<const Thyra::ProductMultiVectorBase<SC> > X
369 = Teuchos::rcp_dynamic_cast<const Thyra::ProductMultiVectorBase<SC> > ( rcpFromRef(X_in) );
370
371 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > Y
372 = Teuchos::rcp_dynamic_cast< Thyra::ProductMultiVectorBase<SC> > ( rcpFromPtr(Y_inout) );
373
374
375// Teuchos::RCP<const Thyra::ProductMultiVectorBase<SC> >
376// X = Thyra::castOrCreateSingleBlockProductMultiVector<SC>( this->defaultProductDomain_, rcpFromRef(X_in) );
377// Teuchos::RCP<Thyra::ProductMultiVectorBase<SC> >
378// Y = Thyra::nonconstCastOrCreateSingleBlockProductMultiVector<SC>( this->defaultProductRange_, rcpFromPtr(Y_inout) );
379 Y_inout->assign(0.);
380// std::cout << "Xin" << std::endl;
381// X_in.describe(*out,Teuchos::VERB_EXTREME);
382// comm_->barrier(); comm_->barrier(); comm_->barrier();
383//
384// std::cout << "Xrcp" << std::endl;
385// X->describe(*out,Teuchos::VERB_EXTREME);
386// comm_->barrier(); comm_->barrier(); comm_->barrier();
387// std::cout << "Y_inoutin" << std::endl;
388// Y_inout->describe(*out,Teuchos::VERB_EXTREME);
389// comm_->barrier(); comm_->barrier(); comm_->barrier();
390 Teuchos::RCP< const MultiVectorBase< SC > > X_s = X->getMultiVectorBlock(2);
391 Teuchos::RCP< MultiVectorBase< SC > > Y_s = Y->getNonconstMultiVectorBlock(2);
392
393
394 Teuchos::RCP< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > XsTpetra =
395 Teuchos::rcp_dynamic_cast< const Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( X_s );
396// XsTpetra->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
397
398// std::cout << "Xs" << std::endl;
399// X_s->describe(*out,Teuchos::VERB_EXTREME);
400// comm_->barrier(); comm_->barrier(); comm_->barrier();
401
402// typedef Thyra::TpetraVector< SC, LO, GO, NO > TpetraVector_Type;
403// typedef Teuchos::RCP<TpetraVector_Type> TpetraVectorPtr_Type;
404// typedef Teuchos::RCP<const TpetraVector_Type> TpetraVectorConstPtr_Type;
405//
406// TpetraVectorConstPtr_Type tpetraThyraXs = Teuchos::rcp_dynamic_cast<const TpetraVector_Type>( X_s );
407//
408// std::cout << "Xs tpetra:" << std::endl;
409// tpetraThyraXs->getConstTpetraVector()->describe(*out,Teuchos::VERB_EXTREME);
410//
411// comm_->barrier(); comm_->barrier(); comm_->barrier();
412
413
414 // apply solid preconditioner
415 if (useSolidPreconditioner_)
416 sInv_->apply(NOTRANS, *X_s, Y_s.ptr(), 1., 0.);
417 else
418 assign(Y_s.ptr(), *X_s);
419
420// std::cout << "Ys" << std::endl;
421// Y_s->describe(*out,Teuchos::VERB_EXTREME);
422// comm_->barrier(); comm_->barrier(); comm_->barrier();
423 Teuchos::RCP< const MultiVectorBase< SC > > X_g;
424 Teuchos::RCP< MultiVectorBase< SC > > Y_g;
425 // apply geometry preconditioner
426 if ( !gInv_.is_null() ) {
427 X_g = X->getMultiVectorBlock(4);
428 Y_g = Y->getNonconstMultiVectorBlock(4);
429
430 assign(Y_g.ptr(), *X_g);
431
432 C4_->apply(NOTRANS, *Y_s, Y_g.ptr(), -1., 1.);
433
434 gInv_->apply(NOTRANS, *Y_g, Y_g.ptr(), 1., 0.);
435 }
436
437 Teuchos::RCP< const MultiVectorBase< SC > > X_l = X->getMultiVectorBlock(3);
438 Teuchos::RCP< MultiVectorBase< SC > > Y_l = Y->getNonconstMultiVectorBlock(3);
439
440 assign(Y_l.ptr(), *X_l);
441// std::cout << "Yl pre c2 apply" << std::endl;
442// Y_l->describe(*out,Teuchos::VERB_EXTREME);
443// comm_->barrier(); comm_->barrier(); comm_->barrier();
444
445 C2_->apply( NOTRANS, *Y_s, Y_l.ptr(), -1., 1. );
446// C2_->apply( NOTRANS, *Y_s, Y_l.ptr(), 1., 0. );
447// std::cout << "Yl after c2 apply" << std::endl;
448// Y_l->describe(*out,Teuchos::VERB_EXTREME);
449// comm_->barrier(); comm_->barrier(); comm_->barrier();
450// scale(-1., Y_l.ptr());
451// update(1., *X_l, Y_l.ptr());
452// std::cout << "Yl2" << std::endl;
453// Y_l->describe(*out,Teuchos::VERB_EXTREME);
454// comm_->barrier(); comm_->barrier(); comm_->barrier();
455
456 Teuchos::RCP< const MultiVectorBase< SC > > X_fv = X->getMultiVectorBlock(0);
457 Teuchos::RCP< MultiVectorBase< SC > > Y_fv = Y->getNonconstMultiVectorBlock(0);
458 assign(Y_fv.ptr(), *X_fv);
459 // if (Z_fv_.is_null())
460// Z_fv_ = X_fv->clone_mv();
461// else
462
463
464
465
466 Teuchos::RCP< const MultiVectorBase< SC > > X_fp = X->getMultiVectorBlock(1);
467 Teuchos::RCP< MultiVectorBase< SC > > Y_fp = Y->getNonconstMultiVectorBlock(1);
468 assign(Y_fp.ptr(), *X_fp);
469// if (Z_fp_.is_null())
470// Z_fp_ = X_fv->clone_mv();
471// else
472// assign(Z_fp_.ptr(), *X_fp);
473
474// std::cout << "zp" << std::endl;
475// Z_fp->describe(*out,Teuchos::VERB_EXTREME);
476// comm_->barrier(); comm_->barrier(); comm_->barrier();
477// std::cout << "Y-fp set:" << std::endl;
478// Y_fp->describe(*out,Teuchos::VERB_EXTREME);
479// comm_->barrier(); comm_->barrier(); comm_->barrier();
480 if (!shape_v_.is_null() && !shape_p_.is_null()) {
481 shape_v_->apply(NOTRANS, *Y_g, Y_fv.ptr(), -1., 1.);
482 shape_p_->apply(NOTRANS, *Y_g, Y_fp.ptr(), -1., 1.);
483 }
484
485
486 if (Z_fv_.is_null())
487 Z_fv_ = Y_fv->clone_mv();
488 else
489 assign(Z_fv_.ptr(), *Y_fv);
490
491// std::cout << "Z_fv_ set:" << std::endl;
492// Z_fv_->describe(*out,Teuchos::VERB_EXTREME);
493// comm_->barrier(); comm_->barrier(); comm_->barrier();
494//
495// std::cout << "Z_fp_ set:" << std::endl;
496// Y_fp->describe(*out,Teuchos::VERB_EXTREME);
497// comm_->barrier(); comm_->barrier(); comm_->barrier();
498
499
500 // fluid condensation
501 if (tmp_l_.is_null())
502 tmp_l_ = Y_l->clone_mv();
503
504 C1_->apply(NOTRANS, *Y_fv, tmp_l_.ptr(), 1., 0);
505 C1T_->apply(NOTRANS, *tmp_l_, Y_fv.ptr(), -1., 1.);
506
507 C1T_->apply(NOTRANS, *Y_l, Y_fv.ptr(), 1., 1.);
508
509// std::cout << "W_fv" << std::endl;
510// Y_fv->describe(*out,Teuchos::VERB_EXTREME);
511// comm_->barrier(); comm_->barrier(); comm_->barrier();
512
513 if (productRangeFluid_.is_null()) {
514 Teuchos::Array< Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesRangeFluid( 2 );
515
516 vectorSpacesRangeFluid[0] = X_fv->range();
517 vectorSpacesRangeFluid[1] = X_fp->range();
518
519 productRangeFluid_ = Thyra::productVectorSpace<SC>( vectorSpacesRangeFluid );
520 }
521
522 Teuchos::Array< Teuchos::RCP< Thyra::MultiVectorBase< SC > > > X_fluid( 2 );
523 Teuchos::Array< Teuchos::RCP< Thyra::MultiVectorBase< SC > > > Y_fluid( 2 );
524
525 X_fluid[0] = Y_fv;
526 X_fluid[1] = Y_fp;
527
528 Y_fluid[0] = Y_fv;
529 Y_fluid[1] = Y_fp;
530
531
532// std::cout << "Y_fp before mono" << std::endl;
533// Y_fp->describe(*out,Teuchos::VERB_EXTREME);
534// comm_->barrier(); comm_->barrier(); comm_->barrier();
535
536 if (useFluidPreconditioner_){
537
538 if (fluidPrecMonolithic_) {
539 Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > fMonoVS = fInv_->domain();
540
541 if ( X_fmono_.is_null() ){
542 X_fmono_ = createMembers( fMonoVS, X_fluid[0]->domain()->dim() );
543 Y_fmono_ = createMembers( fMonoVS, Y_fluid[0]->domain()->dim() );
544 }
545 //We can/should speedup this process
546 copyToMono(X_fluid);
547 fInv_->apply(NOTRANS, *X_fmono_, Y_fmono_.ptr(), 1., 0.);
548 copyFromMono(Y_fluid);
549 }
550 else{
551 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > prodX_f = Thyra::defaultProductMultiVector<SC>( productRangeFluid_, X_fluid );
552 Teuchos::RCP< Thyra::ProductMultiVectorBase<SC> > prodY_f = Thyra::defaultProductMultiVector<SC>( productRangeFluid_, Y_fluid );
553
554 Teuchos::RCP< MultiVectorBase<SC> > X_f = Teuchos::rcp_dynamic_cast<MultiVectorBase< SC > >(prodX_f);
555 Teuchos::RCP< MultiVectorBase<SC> > Y_f = Teuchos::rcp_dynamic_cast<MultiVectorBase< SC > >(prodY_f);
556
557 fInv_->apply(NOTRANS, *X_f, Y_f.ptr(), 1., 0.);
558 }
559 }
560 else{
561 assign(Y_fv.ptr(), *X_fv);
562 assign(Y_fp.ptr(), *X_fp);
563 }
564
565// std::cout << "Y_fv after mono" << std::endl;
566// Y_fv->describe(*out,Teuchos::VERB_EXTREME);
567// comm_->barrier(); comm_->barrier(); comm_->barrier();
568// std::cout << "Y_fp after mono" << std::endl;
569// Y_fp->describe(*out,Teuchos::VERB_EXTREME);
570// comm_->barrier(); comm_->barrier(); comm_->barrier();
571
572
573 fBT_->apply(NOTRANS, *Y_fp, Z_fv_.ptr(), -1., 1.);
574
575// std::cout << "W_fv after BT" << std::endl;
576// W_fv->describe(*out,Teuchos::VERB_EXTREME);
577// comm_->barrier(); comm_->barrier(); comm_->barrier();
578
579 fF_->apply(NOTRANS, *Y_fv, Z_fv_.ptr(), -1., 1.);
580
581// std::cout << "W_fv after F" << std::endl;
582// W_fv->describe(*out,Teuchos::VERB_EXTREME);
583// comm_->barrier(); comm_->barrier(); comm_->barrier();
584
585// std::cout << "W_fv after update" << std::endl;
586// Z_fv_->describe(*out,Teuchos::VERB_EXTREME);
587// comm_->barrier(); comm_->barrier(); comm_->barrier();
588
589 C1_->apply(NOTRANS, *Z_fv_, Y_l.ptr(), 1., 0.);
590
591// std::cout << "Y_l after C1" << std::endl;
592// Y_l->describe(*out,Teuchos::VERB_EXTREME);
593// comm_->barrier(); comm_->barrier(); comm_->barrier();
594// std::cout << "Y_inoutout" << std::endl;
595// Y_inout->describe(*out,Teuchos::VERB_EXTREME);
596// comm_->barrier(); comm_->barrier(); comm_->barrier();
597// std::cout << "FaCSI done!" << std::endl;
598// TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,"End of FaCSI apply.");
599
600
601
602// Teuchos::RCP< MultiVectorBase< SC > > Y_fv = Y->getNonconstMultiVectorBlock(0);
603// Teuchos::RCP< MultiVectorBase< SC > > Y_fp = Y->getNonconstMultiVectorBlock(1);
604// Teuchos::RCP< MultiVectorBase< SC > > Y_s = Y->getNonconstMultiVectorBlock(2);
605// Teuchos::RCP< MultiVectorBase< SC > > Y_l = Y->getNonconstMultiVectorBlock(3);
606
607 Teuchos::RCP< Thyra::TpetraMultiVector< SC, LO, GO, NO > > Y_fvT =
608 Teuchos::rcp_dynamic_cast< Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( Y_fv );
609 Teuchos::RCP< Thyra::TpetraMultiVector< SC, LO, GO, NO > > Y_fpT =
610 Teuchos::rcp_dynamic_cast< Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( Y_fp );
611 Teuchos::RCP< Thyra::TpetraMultiVector< SC, LO, GO, NO > > Y_sT =
612 Teuchos::rcp_dynamic_cast< Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( Y_s );
613 Teuchos::RCP< Thyra::TpetraMultiVector< SC, LO, GO, NO > > Y_lT =
614 Teuchos::rcp_dynamic_cast< Thyra::TpetraMultiVector< SC, LO, GO, NO > > ( Y_l );
615
616// std::cout << "Yfv:" << std::endl;
617// Y_fvT->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
618// comm_->barrier(); comm_->barrier(); comm_->barrier();
619// std::cout << "Yfp:" << std::endl;
620// Y_fpT->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
621// comm_->barrier(); comm_->barrier(); comm_->barrier();
622// std::cout << "Ys:" << std::endl;
623// Y_sT->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
624// comm_->barrier(); comm_->barrier(); comm_->barrier();
625// std::cout << "Yl:" << std::endl;
626// Y_lT->getConstTpetraMultiVector()->describe(*out,Teuchos::VERB_EXTREME);
627// comm_->barrier(); comm_->barrier(); comm_->barrier();
628
629
630
631
632 }
633}
634// private
635template<class SC, class LO, class GO, class NO>
636void PrecOpFaCSI<SC,LO,GO,NO>::copyToMono( Teuchos::Array< Teuchos::RCP< Thyra::MultiVectorBase< SC > > > X_fluid ) const{
637
638 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<SC> > mpiVS_v = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<SC> >(X_fluid[0]->range());
639 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<SC> > mpiVS_p = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<SC> >(X_fluid[1]->range());
640
641 const LO localOffset_v = mpiVS_v->localOffset();
642 const LO localSubDim_v = mpiVS_v->localSubDim();
643
644 const LO localOffset_p = mpiVS_p->localOffset();
645 const LO localSubDim_p = mpiVS_p->localSubDim();
646
647 Teuchos::RCP<Thyra::DetachedMultiVectorView<SC> > thyData_v =
648 Teuchos::rcp(new Thyra::DetachedMultiVectorView<SC>( X_fluid[0] ,Range1D(localOffset_v,localOffset_v+localSubDim_v-1) ) );
649
650 Teuchos::RCP<Thyra::DetachedMultiVectorView<SC> > thyData_p =
651 Teuchos::rcp(new Thyra::DetachedMultiVectorView<SC>( X_fluid[1] ,Range1D(localOffset_p,localOffset_p+localSubDim_p-1) ) );
652
653 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<SC> > mpiVS_mono = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<SC> >(X_fmono_->range());
654
655 const LO localOffset_mono = mpiVS_mono->localOffset();
656 const LO localSubDim_mono = mpiVS_mono->localSubDim();
657
658 Teuchos::RCP<Thyra::DetachedMultiVectorView<SC> > thyData_fluid =
659 Teuchos::rcp(new Thyra::DetachedMultiVectorView<SC>( X_fmono_ ,Range1D(localOffset_mono,localOffset_mono+localSubDim_mono-1) ) );
660
661 for (int j=0; j<X_fluid[0]->domain()->dim(); j++) {
662
663 for (LO i=0; i < localSubDim_v; i++)
664 (*thyData_fluid)(i,j) = (*thyData_v)(i,j);
665
666 for (LO i=0; i<localSubDim_p; i++)
667 (*thyData_fluid)( i + localSubDim_v, j ) = (*thyData_p)(i,j);
668
669 }
670}
671
672template<class SC, class LO, class GO, class NO>
673void PrecOpFaCSI<SC,LO,GO,NO>::copyFromMono(Teuchos::Array< Teuchos::RCP< Thyra::MultiVectorBase< SC > > > Y_fluid) const{
674
675
676 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<SC> > mpiVS_v = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<SC> >(Y_fluid[0]->range());
677 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<SC> > mpiVS_p = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<SC> >(Y_fluid[1]->range());
678
679 const LO localOffset_v = mpiVS_v->localOffset();
680 const LO localSubDim_v = mpiVS_v->localSubDim();
681
682 const LO localOffset_p = mpiVS_p->localOffset();
683 const LO localSubDim_p = mpiVS_p->localSubDim();
684
685 Teuchos::RCP<Thyra::DetachedMultiVectorView<SC> > thyData_v =
686 Teuchos::rcp(new Thyra::DetachedMultiVectorView<SC>( Y_fluid[0] ,Range1D(localOffset_v,localOffset_v+localSubDim_v-1) ) );
687
688 Teuchos::RCP<Thyra::DetachedMultiVectorView<SC> > thyData_p =
689 Teuchos::rcp(new Thyra::DetachedMultiVectorView<SC>( Y_fluid[1] ,Range1D(localOffset_p,localOffset_p+localSubDim_p-1) ) );
690
691 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<SC> > mpiVS_mono = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<SC> >(Y_fmono_->range());
692
693 const LO localOffset_mono = mpiVS_mono->localOffset();
694 const LO localSubDim_mono = mpiVS_mono->localSubDim();
695
696 Teuchos::RCP<Thyra::DetachedMultiVectorView<SC> > thyData_fluid =
697 Teuchos::rcp(new Thyra::DetachedMultiVectorView<SC>( Y_fmono_ ,Range1D(localOffset_mono,localOffset_mono+localSubDim_mono-1) ) );
698
699 for (int j=0; j<Y_fluid[0]->domain()->dim(); j++) {
700
701 for (LO i=0; i < localSubDim_v; i++)
702 (*thyData_v)(i,j) = (*thyData_fluid)(i,j);
703
704 for (LO i=0; i<localSubDim_p; i++)
705 (*thyData_p)(i,j) = (*thyData_fluid)( i + localSubDim_v, j );
706
707 }
708}
709
710}
711
712#endif
FEDD::PrecOpFaCSI::applyImpl
virtual void applyImpl(const Thyra::EOpTransp M_trans, const Thyra::MultiVectorBase< SC > &X, const Teuchos::Ptr< Thyra::MultiVectorBase< SC > > &Y, const SC alpha, const SC beta) const
Definition PrecOpFaCSI_def.hpp:204
FEDD::PreconditionerOperator
Definition PreconditionerOperator_decl.hpp:23
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5