Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
PrecOpFaCSI_def.hpp
1#ifndef PrecOpFaCSI_DEF_hpp
2#define PrecOpFaCSI_DEF_hpp
3
4#include <Thyra_TpetraMultiVector_decl.hpp>
5#include <Teuchos_VerboseObject.hpp>
6#include <Teuchos_Range1D.hpp>
7
16
17namespace FEDD {
18
19// Constructors
20
21template<class SC, class LO, class GO, class NO>
22PrecOpFaCSI<SC,LO,GO,NO>::PrecOpFaCSI()
23:PreconditionerOperator<SC,LO,GO,NO>(),
24fluidPrecMonolithic_(false),
25useSolidPreconditioner_(true)
26{
27// TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Still a problem for FaCSI, cant be used yet.");
28}
29
30template<class SC, class LO, class GO, class NO>
31PrecOpFaCSI<SC,LO,GO,NO>::PrecOpFaCSI(CommConstPtr_Type comm, bool fluidPrecMonolithic, bool useFluidPreconditioner, bool useSolidPreconditioner, bool onlyDiagonal)
32:PreconditionerOperator<SC,LO,GO,NO>(),
33fluidPrecMonolithic_(fluidPrecMonolithic),
34useFluidPreconditioner_(useFluidPreconditioner),
35useSolidPreconditioner_(useSolidPreconditioner),
36onlyDiagonal_(onlyDiagonal)
37{
38 comm_=comm;
39// TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Still a problem for FaCSI, cant be used yet.");
40}
41
42template<class SC, class LO, class GO, class NO>
43void PrecOpFaCSI<SC,LO,GO,NO>::setGE(ThyraLinOpPtr_Type C1,
44 ThyraLinOpPtr_Type C1T,
45 ThyraLinOpPtr_Type C2,
46 ThyraLinOpPtr_Type sInv,
47 ThyraLinOpPtr_Type fInv,
48 ThyraLinOpPtr_Type fF,
49 ThyraLinOpPtr_Type fBT){
50
51 setC1(C1);
52 setC1T(C1T);
53 setC2(C2);
54 setStructInv(sInv);
55 setFluidInv(fInv);
56 setFluidF(fF);
57 setFluidBT(fBT);
58
59 initialize();
60}
61
62template<class SC, class LO, class GO, class NO>
63void PrecOpFaCSI<SC,LO,GO,NO>::setGI(ThyraLinOpPtr_Type C1,
64 ThyraLinOpPtr_Type C1T,
65 ThyraLinOpPtr_Type C2,
66 ThyraLinOpPtr_Type C4,
67 ThyraLinOpPtr_Type sInv,
68 ThyraLinOpPtr_Type fInv,
69 ThyraLinOpPtr_Type fF,
70 ThyraLinOpPtr_Type fBT,
71 ThyraLinOpPtr_Type gInv){
72
73 setC1(C1);
74 setC1T(C1T);
75 setC2(C2);
76 setC4(C4);
77 setStructInv(sInv);
78 setFluidInv(fInv);
79 setFluidF(fF);
80 setFluidBT(fBT);
81 setGeoInv(gInv);
82
83 initialize();
84}
85
86template<class SC, class LO, class GO, class NO>
87void PrecOpFaCSI<SC,LO,GO,NO>::setGIShape(ThyraLinOpPtr_Type C1,
88 ThyraLinOpPtr_Type C1T,
89 ThyraLinOpPtr_Type C2,
90 ThyraLinOpPtr_Type C4,
91 ThyraLinOpPtr_Type sInv,
92 ThyraLinOpPtr_Type fInv,
93 ThyraLinOpPtr_Type fF,
94 ThyraLinOpPtr_Type fBT,
95 ThyraLinOpPtr_Type gInv,
96 ThyraLinOpPtr_Type shape_v,
97 ThyraLinOpPtr_Type shape_p){
98
99 setC1(C1);
100 setC1T(C1T);
101 setC2(C2);
102 setC4(C4);
103 setStructInv(sInv);
104 setFluidInv(fInv);
105 setFluidF(fF);
106 setFluidBT(fBT);
107 setGeoInv(gInv);
108 setShapeDeriv(shape_v, shape_p);
109
110 initialize();
111}
112
113
114template<class SC, class LO, class GO, class NO>
115void PrecOpFaCSI<SC,LO,GO,NO>::setC1(ThyraLinOpPtr_Type C1){
116 C1_ = C1;
117}
118template<class SC, class LO, class GO, class NO>
119void PrecOpFaCSI<SC,LO,GO,NO>::setC1T(ThyraLinOpPtr_Type C1T){
120 C1T_ = C1T;
121}
122template<class SC, class LO, class GO, class NO>
123void PrecOpFaCSI<SC,LO,GO,NO>::setC2(ThyraLinOpPtr_Type C2){
124 C2_ = C2;
125}
126template<class SC, class LO, class GO, class NO>
127void PrecOpFaCSI<SC,LO,GO,NO>::setC4(ThyraLinOpPtr_Type C4){
128 C4_ = C4;
129}
130template<class SC, class LO, class GO, class NO>
131void PrecOpFaCSI<SC,LO,GO,NO>::setShapeDeriv(ThyraLinOpPtr_Type shape_v, ThyraLinOpPtr_Type shape_p){
132 shape_v_ = shape_v;
133 shape_p_ = shape_p;
134}
135template<class SC, class LO, class GO, class NO>
136void PrecOpFaCSI<SC,LO,GO,NO>::setStructInv(ThyraLinOpPtr_Type sInv){
137 sInv_ = sInv;
138}
139template<class SC, class LO, class GO, class NO>
140void PrecOpFaCSI<SC,LO,GO,NO>::setFluidInv(ThyraLinOpPtr_Type fInv){
141 fInv_ = fInv;
142}
143template<class SC, class LO, class GO, class NO>
144void PrecOpFaCSI<SC,LO,GO,NO>::setGeoInv(ThyraLinOpPtr_Type gInv){
145 gInv_ = gInv;
146}
147template<class SC, class LO, class GO, class NO>
148void PrecOpFaCSI<SC,LO,GO,NO>::setFluidF(ThyraLinOpPtr_Type fF){
149 fF_ = fF;
150}
151template<class SC, class LO, class GO, class NO>
152void PrecOpFaCSI<SC,LO,GO,NO>::setFluidBT(ThyraLinOpPtr_Type fBT){
153 fBT_ = fBT;
154}
155
156template<class SC, class LO, class GO, class NO>
157void PrecOpFaCSI<SC,LO,GO,NO>::initialize(){
158 TEUCHOS_TEST_FOR_EXCEPTION(fInv_.is_null(), std::runtime_error,"Can not initialize FaCSI preconditioner: Fluid preconditioner not set.");
159 TEUCHOS_TEST_FOR_EXCEPTION(sInv_.is_null(), std::runtime_error,"Can not initialize FaCSI preconditioner: Structure preconditioner not set.");
160 TEUCHOS_TEST_FOR_EXCEPTION(C1_.is_null(), std::runtime_error,"Can not initialize FaCSI preconditioner: C1 not set.");
161 Teuchos::Array< Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesRange( 4 );
162 Teuchos::Array< Teuchos::RCP< const Thyra::VectorSpaceBase< SC > > > vectorSpacesDomain( 4 );
163 vectorSpacesRange[0] = fF_->range();
164 vectorSpacesRange[1] = fBT_->domain();
165 vectorSpacesRange[2] = sInv_->range();
166 vectorSpacesRange[3] = C1_->range();
167
168 vectorSpacesDomain[0] = fF_->domain();
169 vectorSpacesDomain[1] = fBT_->domain();
170 vectorSpacesDomain[2] = sInv_->domain();
171 vectorSpacesDomain[3] = C1T_->domain();
172 if ( !gInv_.is_null() ) {
173 vectorSpacesRange.push_back( gInv_->range() );
174 vectorSpacesDomain.push_back( gInv_->domain() );
175 }
176 // defaultProductRange_;
177 // Teuchos::RCP<const Thyra::DefaultProductVectorSpace<SC> > defaultProductDomain_;
178
179 Teuchos::RCP<const Thyra::DefaultProductVectorSpace<SC> > pR = Thyra::productVectorSpace<SC>( vectorSpacesRange );
180 Teuchos::RCP<const Thyra::DefaultProductVectorSpace<SC> > pD = Thyra::productVectorSpace<SC>( vectorSpacesDomain );
181// Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > pVSR = Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<SC> >(pR);
182// Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > pVSD = Teuchos::rcp_dynamic_cast<const Thyra::VectorSpaceBase<SC> >(pD);
183//
184// this->defaultProductRange_ = Thyra::multiVectorProductVectorSpace<SC>( pVSR , 1);
185// this->defaultProductDomain_ = Thyra::multiVectorProductVectorSpace<SC>( pVSD, 1);
186
187 this->defaultProductRange_ = pR;
188 this->defaultProductDomain_ = pD;
189}
190
191template<class SC, class LO, class GO, class NO>
192void PrecOpFaCSI<SC,LO,GO,NO>::applyIt(
193 const Thyra::EOpTransp M_trans,
194 const Thyra::MultiVectorBase<SC> &X_in,
195 const Thyra::Ptr<Thyra::MultiVectorBase<SC> > &Y_inout,
196 const SC alpha,
197 const SC beta
198 ) const
199{
200 applyImpl(M_trans, X_in, Y_inout, alpha, beta);
201
202}
203
204template<class SC, class LO, class GO, class NO>
205void PrecOpFaCSI<SC,LO,GO,NO>::applyImpl(
206 const Thyra::EOpTransp M_trans,
207 const Thyra::MultiVectorBase<SC> &X_in,
208 const Thyra::Ptr<Thyra::MultiVectorBase<SC> > &Y_inout,
209 const SC alpha,
210 const SC beta
211 ) const
212{
213 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
214
215 typedef Teuchos::ScalarTraits<SC> ST;
216 typedef Teuchos::RCP<Thyra::MultiVectorBase<SC> > MultiVectorPtr;
217 typedef Teuchos::RCP<const Thyra::MultiVectorBase<SC> > ConstMultiVectorPtr;
218 typedef Teuchos::RCP<const Thyra::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> > ( Teuchos::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 Thyra::MultiVectorBase< SC > > X_s = X->getMultiVectorBlock(2);
232 Teuchos::RCP< Thyra::MultiVectorBase< SC > > Y_s = Y->getNonconstMultiVectorBlock(2);
233
234
235 Teuchos::RCP< const Thyra::MultiVectorBase< SC > > X_fv = X->getMultiVectorBlock(0);
236 Teuchos::RCP< Thyra::MultiVectorBase< SC > > Y_fv = Y->getNonconstMultiVectorBlock(0);
237 Teuchos::RCP< const Thyra::MultiVectorBase< SC > > X_fp = X->getMultiVectorBlock(1);
238 Teuchos::RCP< Thyra::MultiVectorBase< SC > > Y_fp = Y->getNonconstMultiVectorBlock(1);
239
240
241 Teuchos::RCP< Thyra::MultiVectorBase< SC > > Y_l = Y->getNonconstMultiVectorBlock(3);
242 Teuchos::RCP<const Thyra::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(Thyra::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(Thyra::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< Thyra::MultiVectorBase<SC> > X_f = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase< SC > >(prodX_f);
334// Teuchos::RCP< Thyra::MultiVectorBase<SC> > Y_f = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase< SC > >(prodY_f);
335//
336// fInv_->apply(Thyra::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> > ( Teuchos::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_, Teuchos::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 Thyra::MultiVectorBase< SC > > X_s = X->getMultiVectorBlock(2);
391 Teuchos::RCP< Thyra::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(Thyra::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 Thyra::MultiVectorBase< SC > > X_g;
424 Teuchos::RCP< Thyra::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(Thyra::NOTRANS, *Y_s, Y_g.ptr(), -1., 1.);
433
434 gInv_->apply(Thyra::NOTRANS, *Y_g, Y_g.ptr(), 1., 0.);
435 }
436
437 Teuchos::RCP< const Thyra::MultiVectorBase< SC > > X_l = X->getMultiVectorBlock(3);
438 Teuchos::RCP< Thyra::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( Thyra::NOTRANS, *Y_s, Y_l.ptr(), -1., 1. );
446// C2_->apply( Thyra::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 Thyra::MultiVectorBase< SC > > X_fv = X->getMultiVectorBlock(0);
457 Teuchos::RCP< Thyra::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 Thyra::MultiVectorBase< SC > > X_fp = X->getMultiVectorBlock(1);
467 Teuchos::RCP< Thyra::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(Thyra::NOTRANS, *Y_g, Y_fv.ptr(), -1., 1.);
482 shape_p_->apply(Thyra::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(Thyra::NOTRANS, *Y_fv, tmp_l_.ptr(), 1., 0);
505 C1T_->apply(Thyra::NOTRANS, *tmp_l_, Y_fv.ptr(), -1., 1.);
506
507 C1T_->apply(Thyra::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(Thyra::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< Thyra::MultiVectorBase<SC> > X_f = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase< SC > >(prodX_f);
555 Teuchos::RCP< Thyra::MultiVectorBase<SC> > Y_f = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase< SC > >(prodY_f);
556
557 fInv_->apply(Thyra::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(Thyra::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(Thyra::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(Thyra::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] ,Teuchos::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] ,Teuchos::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_ ,Teuchos::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] ,Teuchos::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] ,Teuchos::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_ ,Teuchos::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:205
FEDD::PreconditionerOperator
Definition PreconditionerOperator_decl.hpp:23
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36