Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_BlockedTpetraLinearObjFactory_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Panzer: A partial differential equation assembly
4// engine for strongly coupled complex multiphysics systems
5//
6// Copyright 2011 NTESS and the Panzer contributors.
7// SPDX-License-Identifier: BSD-3-Clause
8// *****************************************************************************
9// @HEADER
10
11#ifndef __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
12#define __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
13
14// Panzer
16#include <utility>
17#ifdef PANZER_HAVE_EPETRA_STACK
18#include "Panzer_EpetraVector_Write_GlobalEvaluationData.hpp" // JMG: Remove this eventually.
19#endif
22
23#include "KokkosSparse_SortCrs.hpp"
24
25// Thyra
26#include "Thyra_DefaultBlockedLinearOp.hpp"
27#include "Thyra_DefaultProductVectorSpace.hpp"
28#include "Thyra_SpmdVectorBase.hpp"
29#include "Thyra_TpetraLinearOp.hpp"
30#include "Thyra_TpetraThyraWrappers.hpp"
31
32// Tpetra
33#include "Tpetra_CrsMatrix.hpp"
34#include "Tpetra_MultiVector.hpp"
35#include "Tpetra_Vector.hpp"
36
37namespace panzer {
38
39using Teuchos::RCP;
40
41// ************************************************************
42// class BlockedTpetraLinearObjFactory
43// ************************************************************
44
45template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
47BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
48 const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
49 : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
50{
51 for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
52 gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
53
55
56 // build and register the gather/scatter evaluators with
57 // the base class.
59}
60
61template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
63BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
64 const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
65 : gidProviders_(gidProviders), comm_(comm)
66{
68}
69
70template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
74
75// LinearObjectFactory functions
77
78template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
79Teuchos::RCP<LinearObjContainer>
82{
83 std::vector<Teuchos::RCP<const MapType> > blockMaps;
84 std::size_t blockDim = gidProviders_.size();
85 for(std::size_t i=0;i<blockDim;i++)
86 blockMaps.push_back(getMap(i));
87
88 Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
89 container->setMapsForBlocks(blockMaps);
90
91 return container;
92}
93
94template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
95Teuchos::RCP<LinearObjContainer>
98{
99 std::vector<Teuchos::RCP<const MapType> > blockMaps;
100 std::size_t blockDim = gidProviders_.size();
101 for(std::size_t i=0;i<blockDim;i++)
102 blockMaps.push_back(getGhostedMap(i));
103
104 Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
105 container->setMapsForBlocks(blockMaps);
106
107 return container;
108}
109
110template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
113{
114 using Teuchos::is_null;
115
116 typedef LinearObjContainer LOC;
117
118 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
119 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
120
121 // Operations occur if the GLOBAL container has the correct targets!
122 // Users set the GLOBAL continer arguments
123 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
124 globalToGhostThyraVector(b_in.get_x(),b_out.get_x());
125
126 if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
127 globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt());
128
129 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
130 globalToGhostThyraVector(b_in.get_f(),b_out.get_f());
131}
132
133template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
136{
137 using Teuchos::is_null;
138
139 typedef LinearObjContainer LOC;
140
141 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
142 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
143
144 // Operations occur if the GLOBAL container has the correct targets!
145 // Users set the GLOBAL continer arguments
146 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
147 ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x());
148
149 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
150 ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f());
151
152 if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
153 ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
154}
155
156template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
159 const LinearObjContainer & globalBCRows,
160 LinearObjContainer & ghostedObjs,
161 bool zeroVectorRows, bool adjustX) const
162{
163 using Teuchos::RCP;
164 using Teuchos::rcp_dynamic_cast;
166 using Thyra::PhysicallyBlockedLinearOpBase;
167 using Thyra::VectorBase;
169
170 std::size_t blockDim = gidProviders_.size();
171
172 // first cast to block LOCs
173 const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
174 const BTLOC & b_globalBCRows = Teuchos::dyn_cast<const BTLOC>(globalBCRows);
175 BTLOC & b_ghosted = Teuchos::dyn_cast<BTLOC>(ghostedObjs);
176
177 TEUCHOS_ASSERT(b_localBCRows.get_f()!=Teuchos::null);
178 TEUCHOS_ASSERT(b_globalBCRows.get_f()!=Teuchos::null);
179
180 // cast each component as needed to their product form
181 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.get_A());
182 RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_f());
183 RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.get_f(),true);
184 RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.get_f(),true);
185
186 if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
187
188 // sanity check!
189 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
190 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
191 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
192 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
193 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
194
195 for(std::size_t i=0;i<blockDim;i++) {
196 // grab epetra vector
197 RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),true)->getConstTpetraVector();
198 RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),true)->getConstTpetraVector();
199
200 // pull out epetra values
201 RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
202 RCP<VectorType> t_f;
203 if(th_f==Teuchos::null)
204 t_f = Teuchos::null;
205 else
206 t_f = rcp_dynamic_cast<ThyraVector>(th_f,true)->getTpetraVector();
207
208 for(std::size_t j=0;j<blockDim;j++) {
209 RCP<const MapType> map_i = getGhostedMap(i);
210 RCP<const MapType> map_j = getGhostedMap(j);
211
212 // pull out epetra values
213 RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
214
215 // don't do anyting if opertor is null
216 RCP<CrsMatrixType> t_A;
217 if(th_A==Teuchos::null)
218 t_A = Teuchos::null;
219 else {
220 RCP<OperatorType> t_A_op = rcp_dynamic_cast<ThyraLinearOp>(th_A,true)->getTpetraOperator();
221 t_A = rcp_dynamic_cast<CrsMatrixType>(t_A_op,true);
222 }
223
224 // adjust Block operator
225 if(t_A!=Teuchos::null) {
226 t_A->resumeFill();
227 }
228
229 adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
230
231 if(t_A!=Teuchos::null) {
232 //t_A->fillComplete(map_j,map_i);
233 }
234
235 t_f = Teuchos::null; // this is so we only adjust it once on the first pass
236 }
237 }
238}
239
240template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
243 const VectorType & global_bcs,
244 const Teuchos::Ptr<VectorType> & f,
245 const Teuchos::Ptr<CrsMatrixType> & A,
246 bool zeroVectorRows) const
247{
248 if(f==Teuchos::null && A==Teuchos::null)
249 return;
250
251 Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
252
253 Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
254 Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
255
256 TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
257 for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
258 if(global_bcs_array[i]==0.0)
259 continue;
260
261 if(local_bcs_array[i]==0.0 || zeroVectorRows) {
262 // this boundary condition was NOT set by this processor
263
264 // if they exist put 0.0 in each entry
265 if(!Teuchos::is_null(f))
266 f_array[i] = 0.0;
267 if(!Teuchos::is_null(A)) {
268 std::size_t numEntries = 0;
269 std::size_t sz = A->getNumEntriesInLocalRow(i);
270 typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
271 typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
272
273 A->getLocalRowCopy(i,indices,values,numEntries);
274
275 for(std::size_t c=0;c<numEntries;c++)
276 values(c) = 0.0;
277
278 A->replaceLocalValues(i,indices,values);
279 }
280 }
281 else {
282 // this boundary condition was set by this processor
283
284 ScalarT scaleFactor = global_bcs_array[i];
285
286 // if they exist scale linear objects by scale factor
287 if(!Teuchos::is_null(f))
288 f_array[i] /= scaleFactor;
289 if(!Teuchos::is_null(A)) {
290 std::size_t numEntries = 0;
291 std::size_t sz = A->getNumEntriesInLocalRow(i);
292 typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
293 typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
294
295 A->getLocalRowCopy(i,indices,values,numEntries);
296
297 for(std::size_t c=0;c<numEntries;c++)
298 values(c) /= scaleFactor;
299
300 A->replaceLocalValues(i,indices,values);
301 }
302 }
303 }
304}
305
306template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
308applyDirichletBCs(const LinearObjContainer & /* counter */,
309 LinearObjContainer & /* result */) const
310{
311 TEUCHOS_ASSERT(false); // not yet implemented
312}
313
315//
316// buildReadOnlyDomainContainer()
317//
319template<typename Traits, typename ScalarT, typename LocalOrdinalT,
320 typename GlobalOrdinalT, typename NodeT>
321Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
322BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
323 NodeT>::
324buildReadOnlyDomainContainer() const
325{
326 using std::vector;
327 using Teuchos::RCP;
328 using Teuchos::rcp;
331 LocalOrdinalT, GlobalOrdinalT, NodeT>;
332 vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
333 for (int i(0); i < getBlockColCount(); ++i)
334 {
335 auto tvroged = rcp(new TVROGED);
336 tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
337 gedBlocks.push_back(tvroged);
338 }
339 auto ged = rcp(new BVROGED);
340 ged->initialize(getGhostedThyraDomainSpace(), getThyraDomainSpace(),
341 gedBlocks);
342 return ged;
343} // end of buildReadOnlyDomainContainer()
344
345#ifdef PANZER_HAVE_EPETRA_STACK
347//
348// buildWriteDomainContainer()
349//
351template<typename Traits, typename ScalarT, typename LocalOrdinalT,
352 typename GlobalOrdinalT, typename NodeT>
353Teuchos::RCP<WriteVector_GlobalEvaluationData>
354BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
355 NodeT>::
356buildWriteDomainContainer() const
357{
358 using std::logic_error;
359 using Teuchos::rcp;
361 auto ged = rcp(new EVWGED);
362 TEUCHOS_TEST_FOR_EXCEPTION(true, logic_error, "NOT YET IMPLEMENTED")
363 return ged;
364} // end of buildWriteDomainContainer()
365#endif // PANZER_HAVE_EPETRA_STACK
366
367template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
373
374template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
376initializeContainer(int mem,LinearObjContainer & loc) const
377{
378 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
379 initializeContainer(mem,bloc);
380}
381
382template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
385{
386 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
387 initializeGhostedContainer(mem,bloc);
388}
389
390// Generic methods
392
393template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
395initializeContainer(int mem,BTLOC & loc) const
396{
397 typedef LinearObjContainer LOC;
398
399 loc.clear();
400
401 if((mem & LOC::X) == LOC::X)
402 loc.set_x(getThyraDomainVector());
403
404 if((mem & LOC::DxDt) == LOC::DxDt)
405 loc.set_dxdt(getThyraDomainVector());
406
407 if((mem & LOC::F) == LOC::F)
408 loc.set_f(getThyraRangeVector());
409
410 if((mem & LOC::Mat) == LOC::Mat)
411 loc.set_A(getThyraMatrix());
412}
413
414template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
416initializeGhostedContainer(int mem,BTLOC & loc) const
417{
418 typedef LinearObjContainer LOC;
419
420 loc.clear();
421
422 if((mem & LOC::X) == LOC::X)
423 loc.set_x(getGhostedThyraDomainVector());
424
425 if((mem & LOC::DxDt) == LOC::DxDt)
426 loc.set_dxdt(getGhostedThyraDomainVector());
427
428 if((mem & LOC::F) == LOC::F) {
429 loc.set_f(getGhostedThyraRangeVector());
431 }
432
433 if((mem & LOC::Mat) == LOC::Mat) {
434 loc.set_A(getGhostedThyraMatrix());
436 }
437}
438
439template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
441addExcludedPair(int rowBlock,int colBlock)
442{
443 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
444}
445
446template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
448addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs)
449{
450 for(std::size_t i=0;i<exPairs.size();i++)
451 excludedPairs_.insert(exPairs[i]);
452}
453
454template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
455Teuchos::RCP<const GlobalIndexer>
461
462template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
464makeRoomForBlocks(std::size_t blockCnt)
465{
466 maps_.resize(blockCnt);
467 ghostedMaps_.resize(blockCnt);
468 importers_.resize(blockCnt);
469 exporters_.resize(blockCnt);
470}
471
472// Thyra methods
474
475template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
476Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
479{
480 if(domainSpace_==Teuchos::null) {
481 // loop over all vectors and build the vector space
482 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
483 for(std::size_t i=0;i<gidProviders_.size();i++)
484 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
485
486 domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
487 }
488
489 return domainSpace_;
490}
491
492template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
493Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
495getThyraRangeSpace() const
496{
497 if(rangeSpace_==Teuchos::null) {
498 // loop over all vectors and build the vector space
499 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
500 for(std::size_t i=0;i<gidProviders_.size();i++)
501 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
502
503 rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
504 }
505
506 return rangeSpace_;
507}
508
509template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
510Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
512getThyraDomainSpace(int blk) const
513{
514 if(domainSpace_==Teuchos::null) {
515 getThyraDomainSpace();
516 }
517
518 return domainSpace_->getBlock(blk);
519}
520
521template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
522Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
524getThyraRangeSpace(int blk) const
525{
526 if(rangeSpace_==Teuchos::null) {
527 getThyraRangeSpace();
528 }
529
530 return rangeSpace_->getBlock(blk);
531}
532
533template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
534Teuchos::RCP<Thyra::VectorBase<ScalarT> >
537{
538 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
539 Thyra::createMember<ScalarT>(*getThyraDomainSpace());
540 Thyra::assign(vec.ptr(),0.0);
541
542 Teuchos::RCP<Thyra::ProductVectorBase<ScalarT> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<ScalarT> >(vec);
543 for(std::size_t i=0;i<gidProviders_.size();i++) {
544 TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
545 Teuchos::as<int>(getMap(i)->getLocalNumElements()));
546 }
547
548 return vec;
549}
550
551template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
552Teuchos::RCP<Thyra::VectorBase<ScalarT> >
555{
556 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
557 Thyra::createMember<ScalarT>(*getThyraRangeSpace());
558 Thyra::assign(vec.ptr(),0.0);
559
560 return vec;
561}
562
563template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
564Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
566getThyraMatrix() const
567{
568 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
569
570 // get the block dimension
571 std::size_t blockDim = gidProviders_.size();
572
573 // this operator will be square
574 blockedOp->beginBlockFill(blockDim,blockDim);
575
576 // loop over each block
577 for(std::size_t i=0;i<blockDim;i++) {
578 for(std::size_t j=0;j<blockDim;j++) {
579 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
580 // build (i,j) block matrix and add it to blocked operator
581 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
582 blockedOp->setNonconstBlock(i,j,block);
583 }
584 }
585 }
586
587 // all done
588 blockedOp->endBlockFill();
589
590 return blockedOp;
591}
592
593template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
594Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
597{
598 if(ghostedDomainSpace_==Teuchos::null) {
599 // loop over all vectors and build the vector space
600 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
601 for(std::size_t i=0;i<gidProviders_.size();i++)
602 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
603
604 ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
605 }
606
607 return ghostedDomainSpace_;
608}
609
610template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
611Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
614{
615 if(ghostedRangeSpace_==Teuchos::null) {
616 // loop over all vectors and build the vector space
617 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
618 for(std::size_t i=0;i<gidProviders_.size();i++)
619 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
620
621 ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
622 }
623
624 return ghostedRangeSpace_;
625}
626
627template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
628Teuchos::RCP<Thyra::VectorBase<ScalarT> >
631{
632 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
633 Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
634 Thyra::assign(vec.ptr(),0.0);
635
636 return vec;
637}
638
639template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
640Teuchos::RCP<Thyra::VectorBase<ScalarT> >
643{
644 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
645 Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
646 Thyra::assign(vec.ptr(),0.0);
647
648 return vec;
649}
650
651template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
652Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
655{
656 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
657
658 // get the block dimension
659 std::size_t blockDim = gidProviders_.size();
660
661 // this operator will be square
662 blockedOp->beginBlockFill(blockDim,blockDim);
663
664 // loop over each block
665 for(std::size_t i=0;i<blockDim;i++) {
666 for(std::size_t j=0;j<blockDim;j++) {
667 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
668 // build (i,j) block matrix and add it to blocked operator
669 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
670 blockedOp->setNonconstBlock(i,j,block);
671 }
672 }
673 }
674
675 // all done
676 blockedOp->endBlockFill();
677
678 return blockedOp;
679}
680
681template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
683ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
684 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
685{
686 using Teuchos::RCP;
687 using Teuchos::rcp_dynamic_cast;
689
690 std::size_t blockDim = gidProviders_.size();
691
692 // get product vectors
693 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
694 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
695
696 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
697 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
698
699 for(std::size_t i=0;i<blockDim;i++) {
700 // first get each Tpetra vector
701 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
702 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
703
704 // use Tpetra to do global communication
705 ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
706 }
707}
708
709template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
712{
713 using Teuchos::RCP;
714 using Teuchos::rcp_dynamic_cast;
715 using Teuchos::dyn_cast;
717 using Thyra::PhysicallyBlockedLinearOpBase;
718
719 std::size_t blockDim = gidProviders_.size();
720
721 // get product vectors
722 const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
723 PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
724
725 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
726 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
727 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
728 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
729
730 for(std::size_t i=0;i<blockDim;i++) {
731 for(std::size_t j=0;j<blockDim;j++) {
732 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
733 // extract the blocks
734 RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
735 RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
736
737 // sanity check
738 TEUCHOS_ASSERT(th_in!=Teuchos::null);
739 TEUCHOS_ASSERT(th_out!=Teuchos::null);
740
741 // get the epetra version of the blocks
742 RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<const ThyraLinearOp>(th_in,true)->getConstTpetraOperator();
743 RCP<OperatorType> tp_op_out = rcp_dynamic_cast<ThyraLinearOp>(th_out,true)->getTpetraOperator();
744
745 RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<const CrsMatrixType>(tp_op_in,true);
746 RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<CrsMatrixType>(tp_op_out,true);
747
748 // use Tpetra to do global communication
749 ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
750 }
751 }
752 }
753}
754
755template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
757globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
758 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
759{
760 using Teuchos::RCP;
761 using Teuchos::rcp_dynamic_cast;
763
764 std::size_t blockDim = gidProviders_.size();
765
766 // get product vectors
767 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
768 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
769
770 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
771 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
772
773 for(std::size_t i=0;i<blockDim;i++) {
774 // first get each Tpetra vector
775 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
776 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
777
778 // use Tpetra to do global communication
779 globalToGhostTpetraVector(i,*tp_in,*tp_out);
780 }
781}
782
783// Tpetra methods
785
786template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
788ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out) const
789{
790 using Teuchos::RCP;
791
792 // do the global distribution
793 RCP<const ExportType> exporter = getGhostedExport(i);
794 out.putScalar(0.0);
795 out.doExport(in,*exporter,Tpetra::ADD);
796}
797
798template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
800ghostToGlobalTpetraMatrix(int blockRow,const CrsMatrixType & in,CrsMatrixType & out) const
801{
802 using Teuchos::RCP;
803
804 RCP<const MapType> map_i = out.getRangeMap();
805 RCP<const MapType> map_j = out.getDomainMap();
806
807 // do the global distribution
808 RCP<const ExportType> exporter = getGhostedExport(blockRow);
809
810 out.resumeFill();
811 out.setAllToScalar(0.0);
812 out.doExport(in,*exporter,Tpetra::ADD);
813 out.fillComplete(map_j,map_i);
814}
815
816template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
818globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out) const
819{
820 using Teuchos::RCP;
821
822 // do the global distribution
823 RCP<const ImportType> importer = getGhostedImport(i);
824 out.putScalar(0.0);
825 out.doImport(in,*importer,Tpetra::INSERT);
826}
827
828// get the map from the matrix
829template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
830Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
832getMap(int i) const
833{
834 if(maps_[i]==Teuchos::null)
835 maps_[i] = buildTpetraMap(i);
836
837 return maps_[i];
838}
839
840template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
841Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
843getGhostedMap(int i) const
844{
845 if(ghostedMaps_[i]==Teuchos::null)
846 ghostedMaps_[i] = buildTpetraGhostedMap(i);
847
848 return ghostedMaps_[i];
849}
850
851// get the graph of the crs matrix
852template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
853Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
855getGraph(int i,int j) const
856{
857 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
858
859 typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
860 Teuchos::RCP<const CrsGraphType> graph;
861 if(itr==graphs_.end()) {
862 graph = buildTpetraGraph(i,j);
863 graphs_[std::make_pair(i,j)] = graph;
864 }
865 else
866 graph = itr->second;
867
868 TEUCHOS_ASSERT(graph!=Teuchos::null);
869 return graph;
870}
871
872template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
873Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
875getGhostedGraph(int i,int j) const
876{
877 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
878
879 typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
880 Teuchos::RCP<const CrsGraphType> ghostedGraph;
881 if(itr==ghostedGraphs_.end()) {
882 ghostedGraph = buildTpetraGhostedGraph(i,j);
883 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
884 }
885 else
886 ghostedGraph = itr->second;
887
888 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
889 return ghostedGraph;
890}
891
892template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
893Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
895getGhostedImport(int i) const
896{
897 if(importers_[i]==Teuchos::null)
898 importers_[i] = Teuchos::rcp(new ImportType(getMap(i),getGhostedMap(i)));
899
900 return importers_[i];
901}
902
903template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
904Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
906getGhostedExport(int i) const
907{
908 if(exporters_[i]==Teuchos::null)
909 exporters_[i] = Teuchos::rcp(new ExportType(getGhostedMap(i),getMap(i)));
910
911 return exporters_[i];
912}
913
914template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
915Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
917buildTpetraMap(int i) const
918{
919 std::vector<GlobalOrdinalT> indices;
920
921 // get the global indices
922 getGlobalIndexer(i)->getOwnedIndices(indices);
923
924 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
925}
926
927// build the ghosted map
928template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
929Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
931buildTpetraGhostedMap(int i) const
932{
933 std::vector<GlobalOrdinalT> indices;
934
935 // get the global indices
936 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
937
938 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
939}
940
941// get the graph of the crs matrix
942template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
943Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
945buildTpetraGraph(int i,int j) const
946{
947 using Teuchos::RCP;
948 using Teuchos::rcp;
949
950 // build the map and allocate the space for the graph and
951 // grab the ghosted graph
952 RCP<const MapType> map_i = getMap(i);
953 RCP<const MapType> map_j = getMap(j);
954
955 RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,0));
956 RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
957
958 // perform the communication to finish building graph
959 RCP<const ExportType> exporter = getGhostedExport(i);
960 graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
961 graph->fillComplete(map_j,map_i);
962
963 return graph;
964}
965
966template <class LocalOrdinalT>
968 LocalOrdinalT row;
969 LocalOrdinalT col;
970};
971
972template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
973Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
975buildTpetraGhostedGraph(int i,int j) const
976{
977 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::BlockedTpetraLinearObjFactory::buildTpetraGhostedGraph",BTLOF);
978
979 using Teuchos::RCP;
980 using Teuchos::rcp;
981 using exec_space = typename CrsGraphType::execution_space;
982 using memory_space = typename NodeT::memory_space;
983
984 // build the map and allocate the space for the graph and
985 // grab the ghosted graph
986 RCP<const MapType> map_i = getGhostedMap(i);
987 RCP<const MapType> map_j = getGhostedMap(j);
988
989 std::vector<std::string> elementBlockIds;
990
991 Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
992
993 rowProvider = getGlobalIndexer(i);
994 colProvider = getGlobalIndexer(j);
995
996 gidProviders_[0]->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
997 // same element blocks
998
999 RCP<CrsGraphType> graph;
1000 if constexpr (NodeT::is_gpu) {
1001
1002 // Gather elements from mesh blocks.
1003 size_t numElements;
1004 Kokkos::View<LocalOrdinalT *, memory_space> elementsFromBlocks;
1005 {
1006 auto numElementBlocks = elementBlockIds.size();
1007
1008 std::vector<size_t> elementBlockOffsets(numElementBlocks + 1);
1009 elementBlockOffsets[0] = 0;
1010
1011 numElements = 0;
1012 size_t blockNo = 0;
1013 std::vector<std::string>::const_iterator blockItr;
1014 for (blockItr = elementBlockIds.begin();
1015 blockItr != elementBlockIds.end(); ++blockItr) {
1016 std::string blockId = *blockItr;
1017 const std::vector<LocalOrdinalT> &elements =
1018 gidProviders_[0]->getElementBlock(
1019 blockId); // each sub provider "should" have the
1020 // same elements in each element block
1021 numElements += elements.size();
1022 ++blockNo;
1023 elementBlockOffsets[blockNo] = numElements;
1024 }
1025 elementsFromBlocks = Kokkos::View<LocalOrdinalT *, memory_space>(
1026 "elementsFromBlocks", numElements);
1027 blockNo = 0;
1028 for (blockItr = elementBlockIds.begin();
1029 blockItr != elementBlockIds.end(); ++blockItr) {
1030 std::string blockId = *blockItr;
1031 const std::vector<LocalOrdinalT> &elements =
1032 gidProviders_[0]->getElementBlock(
1033 blockId); // each sub provider "should" have the
1034 // same elements in each element block
1035 Kokkos::View<const LocalOrdinalT *, Kokkos::HostSpace,
1036 Kokkos::MemoryTraits<Kokkos::Unmanaged>>
1037 elements_h(elements.data(), elements.size());
1038 Kokkos::deep_copy(
1039 Kokkos::subview(
1040 elementsFromBlocks,
1041 Kokkos::make_pair(elementBlockOffsets[blockNo],
1042 elementBlockOffsets[blockNo + 1])),
1043 elements_h);
1044 ++blockNo;
1045 }
1046 }
1047
1048 {
1049
1050 using local_graph_type = typename CrsGraphType::local_graph_device_type;
1051 using rowptr_type =
1052 typename local_graph_type::row_map_type::non_const_type;
1053 using colidx_type =
1054 typename local_graph_type::entries_type::non_const_type;
1055
1056 using entries_map_type =
1057 Kokkos::UnorderedMap<entry_type<LocalOrdinalT>, void, exec_space>;
1058
1059 auto numRows = map_i->getLocalNumElements();
1060
1061 // We are overallocating by 1 here. This simplifies the logic below. But
1062 // we have to remember to take a subview in the end.
1063 rowptr_type rowptr("ghostedGraph_rowptr", numRows + 2);
1064
1065 auto rowLIDs = rowProvider->getLIDs();
1066 auto colLIDs = colProvider->getLIDs();
1067
1068 auto numDoFsPerElementRow = rowLIDs.extent(1);
1069 auto numDoFsPerElementCol = colLIDs.extent(1);
1070
1071 auto capacity =
1072 numElements * numDoFsPerElementRow * numDoFsPerElementCol;
1073 entries_map_type entries(capacity);
1074
1075 while (true) {
1076
1077 // Loop over all elements and record the entries that we need in the
1078 // graph. Also start building the rowptr.
1079 Kokkos::parallel_for(
1080 "collect_entries", Kokkos::RangePolicy<exec_space>(0, numElements),
1081 KOKKOS_LAMBDA(const LocalOrdinalT k) {
1082 auto elementId = elementsFromBlocks(k);
1084 for (size_t dofNoRow = 0; dofNoRow < numDoFsPerElementRow;
1085 ++dofNoRow) {
1086 entry.row = rowLIDs(elementId, dofNoRow);
1087 for (size_t dofNoCol = 0; dofNoCol < numDoFsPerElementCol;
1088 ++dofNoCol) {
1089 entry.col = colLIDs(elementId, dofNoCol);
1090 auto result = entries.insert(entry);
1091 if (result.success()) {
1092 // New entry. We offset by 2 here.
1093 Kokkos::atomic_inc(&rowptr(entry.row + 2));
1094 }
1095 }
1096 }
1097 });
1098
1099 if (!entries.failed_insert()) {
1100 auto numEntries = entries.size();
1101
1102 // Prefix sum to get offsets.
1103 // This is not the correct rowptr yet.
1104 // We have essentially shifted everything by one position.
1105 // This is useful for when we fill.
1106 typename rowptr_type::value_type numEntries2;
1107 Kokkos::parallel_scan(
1108 "prefix_sum", Kokkos::RangePolicy<exec_space>(0, numRows + 2),
1109 KOKKOS_LAMBDA(const size_t rlid,
1110 typename rowptr_type::value_type &nnz,
1111 const bool is_final) {
1112 nnz += rowptr(rlid);
1113 if (is_final)
1114 rowptr(rlid) = nnz;
1115 },
1116 numEntries2);
1117 TEUCHOS_ASSERT_EQUALITY(numEntries, numEntries2);
1118
1119 // The column indices.
1120 colidx_type colidx(
1121 Kokkos::ViewAllocateWithoutInitializing("ghostedGraph_colidx"),
1122 numEntries);
1123
1124 // Fill the column indices.
1125 // We are using the rowptr to figure out offsets.
1126 // After this step the rowptr is correct.
1127 Kokkos::parallel_for(
1128 "fill", Kokkos::RangePolicy<exec_space>(0, entries.capacity()),
1129 KOKKOS_LAMBDA(const uint32_t c) {
1130 if (entries.valid_at(c)) {
1131 auto entry = entries.key_at(c);
1132 auto offset =
1133 Kokkos::atomic_fetch_inc(&rowptr(entry.row + 1));
1134 colidx(offset) = entry.col;
1135 }
1136 });
1137
1138 // Sort the rows.
1139 KokkosSparse::sort_crs_graph(rowptr, colidx);
1140
1141 // Create the graph
1142 graph = rcp(new CrsGraphType(
1143 map_i, map_j,
1144 Kokkos::subview(rowptr, Kokkos::make_pair((decltype(numRows))0,
1145 numRows + 1)),
1146 colidx));
1147 graph->fillComplete(getMap(j), getMap(i));
1148
1149 break;
1150 } else {
1151 // We ended up not having enough capacity in the UnorderedMap.
1152 // Bump it up and try again.
1153 std::cout << "Insufficient capacity: " << capacity << std::endl;
1154 capacity *= 2;
1155 Kokkos::deep_copy(rowptr, 0);
1156 entries = entries_map_type(capacity);
1157 }
1158 }
1159 }
1160 } else {
1161 // Count number of entries in each row of graph; needed for graph
1162 // constructor
1163 std::vector<size_t> nEntriesPerRow(map_i->getLocalNumElements(), 0);
1164 std::vector<std::string>::const_iterator blockItr;
1165 for (blockItr = elementBlockIds.begin(); blockItr != elementBlockIds.end();
1166 ++blockItr) {
1167 std::string blockId = *blockItr;
1168 // grab elements for this block
1169 const std::vector<LocalOrdinalT> &elements =
1170 gidProviders_[0]->getElementBlock(
1171 blockId); // each sub provider "should" have the
1172 // same elements in each element block
1173
1174 // get information about number of indicies
1175 std::vector<GlobalOrdinalT> row_gids;
1176 std::vector<GlobalOrdinalT> col_gids;
1177
1178 // loop over the elemnts
1179 for (std::size_t elmt = 0; elmt < elements.size(); elmt++) {
1180
1181 rowProvider->getElementGIDs(elements[elmt], row_gids);
1182 colProvider->getElementGIDs(elements[elmt], col_gids);
1183 for (std::size_t row = 0; row < row_gids.size(); row++) {
1184 LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
1185 nEntriesPerRow[lid] += col_gids.size();
1186 }
1187 }
1188 }
1189 Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
1190 graph = rcp(new CrsGraphType(map_i, map_j, nEntriesPerRowView));
1191
1192 // graph information about the mesh
1193 for (blockItr = elementBlockIds.begin(); blockItr != elementBlockIds.end();
1194 ++blockItr) {
1195 std::string blockId = *blockItr;
1196
1197 // grab elements for this block
1198 const std::vector<LocalOrdinalT> &elements =
1199 gidProviders_[0]->getElementBlock(
1200 blockId); // each sub provider "should" have the
1201 // same elements in each element block
1202
1203 // get information about number of indicies
1204 std::vector<GlobalOrdinalT> row_gids;
1205 std::vector<GlobalOrdinalT> col_gids;
1206
1207 // loop over the elemnts
1208 for (std::size_t elmt = 0; elmt < elements.size(); elmt++) {
1209
1210 rowProvider->getElementGIDs(elements[elmt], row_gids);
1211 colProvider->getElementGIDs(elements[elmt], col_gids);
1212 for (std::size_t row = 0; row < row_gids.size(); row++)
1213 graph->insertGlobalIndices(row_gids[row], col_gids);
1214 }
1215 }
1216
1217 // finish filling the graph: Make sure the colmap and row maps coincide to
1218 // minimize calls to LID lookups
1219 graph->fillComplete(getMap(j), getMap(i));
1220 }
1221
1222 return graph;
1223}
1224
1225template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1226Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1228getTpetraMatrix(int i,int j) const
1229{
1230 Teuchos::RCP<const MapType> map_i = getMap(i);
1231 Teuchos::RCP<const MapType> map_j = getMap(j);
1232
1233 Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
1234 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
1235 mat->fillComplete(map_j,map_i);
1236
1237 return mat;
1238}
1239
1240template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1241Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1243getGhostedTpetraMatrix(int i,int j) const
1244{
1245 Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1246 Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1247
1248 Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1249 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
1250 mat->fillComplete(getMap(j),getMap(i));
1251
1252 return mat;
1253}
1254
1255template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1256Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1258getTpetraDomainVector(int i) const
1259{
1260 Teuchos::RCP<const MapType> tMap = getMap(i);
1261 return Teuchos::rcp(new VectorType(tMap));
1262}
1263
1264template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1265Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1268{
1269 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1270 return Teuchos::rcp(new VectorType(tMap));
1271}
1272
1273template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1274Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1276getTpetraRangeVector(int i) const
1277{
1278 Teuchos::RCP<const MapType> tMap = getMap(i);
1279 return Teuchos::rcp(new VectorType(tMap));
1280}
1281
1282template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1283Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1285getGhostedTpetraRangeVector(int i) const
1286{
1287 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1288 return Teuchos::rcp(new VectorType(tMap));
1289}
1290
1291template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1292int
1298
1299template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1300int
1306
1307template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1309beginFill(LinearObjContainer & loc) const
1310{
1311 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1312 if(tloc.get_A()!=Teuchos::null)
1313 tloc.beginFill();
1314}
1315
1316template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1318endFill(LinearObjContainer & loc) const
1319{
1320 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1321 if(tloc.get_A()!=Teuchos::null)
1322 tloc.endFill();
1323}
1324
1325}
1326
1327#endif // __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we're performing.
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraRangeVector() const
Get a range vector.
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
std::vector< Teuchos::RCP< const GlobalIndexer > > gidProviders_
void initializeGhostedContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraDomainVector() const
Get a domain vector.
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
virtual Teuchos::RCP< const CrsGraphType > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
virtual Teuchos::RCP< const MapType > buildTpetraGhostedMap(int i) const
virtual Teuchos::RCP< const MapType > buildTpetraMap(int i) const
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGraph(int i, int j) const
Teuchos::RCP< VectorType > getGhostedTpetraRangeVector(int i) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void makeRoomForBlocks(std::size_t blockCnt)
Allocate the space in the std::vector objects so we can fill with appropriate Tpetra data.
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraRangeVector() const
Get a range vector.
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGhostedGraph(int i, int j) const
virtual Teuchos::RCP< const MapType > getGhostedMap(int i) const
get the ghosted map from the matrix
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
Teuchos::RCP< VectorType > getGhostedTpetraDomainVector(int i) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
void globalToGhostTpetraVector(int i, const VectorType &in, VectorType &out) const
BlockedTpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const BlockedDOFManager > &gidProvider)
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
void ghostToGlobalTpetraVector(int i, const VectorType &in, VectorType &out) const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< ScalarT > &in, Thyra::LinearOpBase< ScalarT > &out) const
void ghostToGlobalTpetraMatrix(int blockRow, const CrsMatrixType &in, CrsMatrixType &out) const
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
virtual Teuchos::RCP< const ImportType > getGhostedImport(int i) const
get importer for converting an overalapped object to a "normal" object
Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > ImportType
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a Thyra operator.
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraDomainVector() const
Get a domain vector.
virtual Teuchos::RCP< const CrsGraphType > getGraph(int i, int j) const
get the graph of the crs matrix
Teuchos::RCP< const panzer::BlockedDOFManager > getGlobalIndexer() const
Teuchos::RCP< Thyra::BlockedLinearOpBase< ScalarT > > getGhostedThyraMatrix() const
Get a Thyra operator.
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors.
This class provides a boundary exchange communication mechanism for vectors.
void buildGatherScatterEvaluators(const BuilderT &builder)