43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_IMPL_HPP
44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_IMPL_HPP
47 #include "Teuchos_Assert.hpp"
49 #include "Phalanx_DataLayout.hpp"
51 #include "Epetra_Map.h"
52 #include "Epetra_Vector.h"
53 #include "Epetra_CrsMatrix.h"
62 #include "Phalanx_DataLayout_MDALayout.hpp"
64 #include "Thyra_SpmdVectorBase.hpp"
65 #include "Thyra_ProductVectorBase.hpp"
66 #include "Thyra_DefaultProductVector.hpp"
67 #include "Thyra_BlockedLinearOpBase.hpp"
68 #include "Thyra_get_Epetra_Operator.hpp"
70 #include "Teuchos_FancyOStream.hpp"
72 #include <unordered_map>
79 template<
typename TRAITS,
typename LO,
typename GO>
85 : rowIndexers_(rIndexers)
86 , colIndexers_(cIndexers)
87 , globalDataKey_(
"Residual Scatter Container")
89 std::string scatterName = p.
get<std::string>(
"Scatter Name");
94 const std::vector<std::string>& names =
101 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
107 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
108 local_side_id_ = p.
get<
int>(
"Local Side ID");
112 scatterFields_.resize(names.size());
113 for (std::size_t eq = 0; eq < names.size(); ++eq) {
117 this->addDependentField(scatterFields_[eq]);
120 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
122 applyBC_.resize(names.size());
123 for (std::size_t eq = 0; eq < names.size(); ++eq) {
125 this->addDependentField(applyBC_[eq]);
130 this->addEvaluatedField(*scatterHolder_);
132 if (p.
isType<std::string>(
"Global Data Key"))
133 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
135 this->setName(scatterName+
" Scatter Residual");
139 template<
typename TRAITS,
typename LO,
typename GO>
144 indexerIds_.resize(scatterFields_.size());
145 subFieldIds_.resize(scatterFields_.size());
148 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
150 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
153 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
157 num_nodes = scatterFields_[0].extent(1);
161 template<
typename TRAITS,
typename LO,
typename GO>
168 using Teuchos::rcp_dynamic_cast;
173 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
179 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_));
180 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc->getDataObject(globalDataKey_));
183 if(blockedContainer!=Teuchos::null)
185 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
186 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
187 else if(epetraContainer!=Teuchos::null)
189 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
190 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
196 template<
typename TRAITS,
typename LO,
typename GO>
202 using Teuchos::ptrFromRef;
203 using Teuchos::rcp_dynamic_cast;
206 using Thyra::SpmdVectorBase;
210 std::string blockId = this->wda(workset).block_id;
211 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
221 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
222 int rowIndexer = indexerIds_[fieldIndex];
223 int subFieldNum = subFieldIds_[fieldIndex];
225 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
226 ->getNonconstLocalData(ptrFromRef(local_dc));
229 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
230 ->getNonconstLocalData(ptrFromRef(local_r));
232 auto subRowIndexer = rowIndexers_[rowIndexer];
233 auto LIDs = subRowIndexer->getLIDs();
234 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
235 Kokkos::deep_copy(LIDs_h, LIDs);
237 auto field = scatterFields_[fieldIndex].get_view();
238 auto field_h = Kokkos::create_mirror_view(
field);
239 Kokkos::deep_copy(field_h,
field);
241 BCFieldType::array_type::HostMirror applyBC_h;
243 auto applyBC = applyBC_[fieldIndex].get_static_view();
244 applyBC_h = Kokkos::create_mirror_view(applyBC);
245 Kokkos::deep_copy(applyBC_h, applyBC);
249 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
250 std::size_t cellLocalId = localCellIds[worksetCellIndex];
254 const std::pair<std::vector<int>,std::vector<int> > & indicePair
255 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
256 const std::vector<int> & elmtOffset = indicePair.first;
257 const std::vector<int> & basisIdMap = indicePair.second;
260 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
261 int offset = elmtOffset[basis];
262 int lid = LIDs_h(cellLocalId, offset);
266 int basisId = basisIdMap[basis];
268 if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
271 local_r[lid] = field_h(worksetCellIndex,basisId);
278 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
281 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
282 int offset = elmtOffset[basis];
283 int lid = LIDs_h(cellLocalId, offset);
287 local_r[lid] = field_h(worksetCellIndex,basis);
302 template<
typename TRAITS,
typename LO,
typename GO>
308 : rowIndexers_(rIndexers)
309 , colIndexers_(cIndexers)
310 , globalDataKey_(
"Residual Scatter Container")
312 std::string scatterName = p.
get<std::string>(
"Scatter Name");
317 const std::vector<std::string>& names =
324 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
330 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
331 local_side_id_ = p.
get<
int>(
"Local Side ID");
335 scatterFields_.resize(names.size());
336 for (std::size_t eq = 0; eq < names.size(); ++eq) {
340 this->addDependentField(scatterFields_[eq]);
343 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
345 applyBC_.resize(names.size());
346 for (std::size_t eq = 0; eq < names.size(); ++eq) {
348 this->addDependentField(applyBC_[eq]);
353 this->addEvaluatedField(*scatterHolder_);
355 if (p.
isType<std::string>(
"Global Data Key"))
356 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
358 this->setName(scatterName+
" Scatter Tangent");
362 template<
typename TRAITS,
typename LO,
typename GO>
367 indexerIds_.resize(scatterFields_.size());
368 subFieldIds_.resize(scatterFields_.size());
371 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
373 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
376 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
380 num_nodes = scatterFields_[0].extent(1);
384 template<
typename TRAITS,
typename LO,
typename GO>
391 using Teuchos::rcp_dynamic_cast;
396 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
402 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_));
403 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc->getDataObject(globalDataKey_));
406 if(blockedContainer!=Teuchos::null)
408 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
409 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
410 else if(epetraContainer!=Teuchos::null)
412 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
413 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
419 template<
typename TRAITS,
typename LO,
typename GO>
427 using Teuchos::ptrFromRef;
428 using Teuchos::rcp_dynamic_cast;
431 using Thyra::SpmdVectorBase;
435 std::string blockId = this->wda(workset).block_id;
436 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
446 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
447 int rowIndexer = indexerIds_[fieldIndex];
448 int subFieldNum = subFieldIds_[fieldIndex];
450 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
451 ->getNonconstLocalData(ptrFromRef(local_dc));
454 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
455 ->getNonconstLocalData(ptrFromRef(local_r));
457 auto subRowIndexer = rowIndexers_[rowIndexer];
459 auto LIDs = subRowIndexer->getLIDs();
460 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
461 Kokkos::deep_copy(LIDs_h, LIDs);
463 auto field = scatterFields_[fieldIndex].get_view();
464 auto field_h = Kokkos::create_mirror_view(
field);
465 Kokkos::deep_copy(field_h,
field);
467 BCFieldType::array_type::HostMirror applyBC_h;
469 auto applyBC = applyBC_[fieldIndex].get_static_view();
470 applyBC_h = Kokkos::create_mirror_view(applyBC);
471 Kokkos::deep_copy(applyBC_h, applyBC);
475 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
476 std::size_t cellLocalId = localCellIds[worksetCellIndex];
480 const std::pair<std::vector<int>,std::vector<int> > & indicePair
481 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
482 const std::vector<int> & elmtOffset = indicePair.first;
483 const std::vector<int> & basisIdMap = indicePair.second;
486 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
487 int offset = elmtOffset[basis];
488 int lid = LIDs_h(cellLocalId, offset);
492 int basisId = basisIdMap[basis];
494 if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
497 local_r[lid] = field_h(worksetCellIndex,basisId).val();
504 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
507 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
508 int offset = elmtOffset[basis];
509 int lid = LIDs_h(cellLocalId, offset);
513 local_r[lid] = field_h(worksetCellIndex,basis).val();
527 template<
typename TRAITS,
typename LO,
typename GO>
533 : rowIndexers_(rIndexers)
534 , colIndexers_(cIndexers)
535 , globalDataKey_(
"Residual Scatter Container")
537 std::string scatterName = p.
get<std::string>(
"Scatter Name");
542 const std::vector<std::string>& names =
551 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
552 local_side_id_ = p.
get<
int>(
"Local Side ID");
555 scatterFields_.resize(names.size());
556 for (std::size_t eq = 0; eq < names.size(); ++eq) {
560 this->addDependentField(scatterFields_[eq]);
563 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
565 applyBC_.resize(names.size());
566 for (std::size_t eq = 0; eq < names.size(); ++eq) {
568 this->addDependentField(applyBC_[eq]);
573 this->addEvaluatedField(*scatterHolder_);
575 if (p.
isType<std::string>(
"Global Data Key"))
576 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
578 if(colIndexers_.size()==0)
579 colIndexers_ = rowIndexers_;
581 this->setName(scatterName+
" Scatter Residual (Jacobian)");
585 template<
typename TRAITS,
typename LO,
typename GO>
590 indexerIds_.resize(scatterFields_.size());
591 subFieldIds_.resize(scatterFields_.size());
594 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
596 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
599 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
603 num_nodes = scatterFields_[0].extent(1);
604 num_eq = scatterFields_.size();
608 template<
typename TRAITS,
typename LO,
typename GO>
614 using Teuchos::rcp_dynamic_cast;
618 = rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
624 blockContainer = rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_),
true);
632 template<
typename TRAITS,
typename LO,
typename GO>
638 using Teuchos::ptrFromRef;
639 using Teuchos::rcp_dynamic_cast;
641 using Thyra::SpmdVectorBase;
644 std::string blockId = this->wda(workset).block_id;
645 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
647 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
649 std::vector<int> blockOffsets;
659 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
660 int rowIndexer = indexerIds_[fieldIndex];
661 int subFieldNum = subFieldIds_[fieldIndex];
665 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
666 ->getNonconstLocalData(ptrFromRef(local_dc));
670 if(r_!=Teuchos::null)
671 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
672 ->getNonconstLocalData(ptrFromRef(local_r));
674 auto subRowIndexer = rowIndexers_[rowIndexer];
675 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
676 const std::vector<int> & subElmtOffset = subIndicePair.first;
677 const std::vector<int> & subBasisIdMap = subIndicePair.second;
679 auto rLIDs = subRowIndexer->getLIDs();
680 auto rLIDs_h = Kokkos::create_mirror_view(rLIDs);
681 Kokkos::deep_copy(rLIDs_h, rLIDs);
683 auto field = scatterFields_[fieldIndex].get_view();
684 auto field_h = Kokkos::create_mirror_view(
field);
685 Kokkos::deep_copy(field_h,
field);
687 BCFieldType::array_type::HostMirror applyBC_h;
689 auto applyBC = applyBC_[fieldIndex].get_static_view();
690 applyBC_h = Kokkos::create_mirror_view(applyBC);
691 Kokkos::deep_copy(applyBC_h, applyBC);
695 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
696 std::size_t cellLocalId = localCellIds[worksetCellIndex];
699 for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
700 int offset = subElmtOffset[basis];
701 int lid = rLIDs_h(cellLocalId, offset);
705 int basisId = subBasisIdMap[basis];
707 if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
711 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
712 int start = blockOffsets[colIndexer];
713 int end = blockOffsets[colIndexer+1];
719 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
720 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
722 if(subJac==Teuchos::null) {
731 jacEpetraBlocks[blockIndex] = subJac;
735 int * rowIndices = 0;
736 double * rowValues = 0;
740 for(
int i=0;i<numEntries;i++)
744 const ScalarT scatterField = field_h(worksetCellIndex,basisId);
746 if(r_!=Teuchos::null)
747 local_r[lid] = scatterField.val();
751 std::vector<double> jacRow(scatterField.size(),0.0);
753 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
754 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
756 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
757 int start = blockOffsets[colIndexer];
758 int end = blockOffsets[colIndexer+1];
763 auto subColIndexer = colIndexers_[colIndexer];
765 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
766 auto cLIDs_h = Kokkos::create_mirror_view(cLIDs);
767 Kokkos::deep_copy(cLIDs_h, cLIDs);
772 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
773 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
776 if(subJac==Teuchos::null) {
785 jacEpetraBlocks[blockIndex] = subJac;
789 int err = subJac->
ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs_h[0]);
791 std::stringstream ss;
792 ss <<
"Failed inserting row: " <<
" (" << lid <<
"): ";
793 for(
int i=0;i<end-start;i++)
794 ss << cLIDs_h[i] <<
" ";
796 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
798 ss <<
"scatter field = ";
799 scatterFields_[fieldIndex].print(ss);
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer >> &ugis, std::vector< int > &blockOffsets)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evaluateFields(typename TRAITS::EvalData d)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer >> &ugis)
ScatterDirichletResidual_BlockedEpetra(const Teuchos::ParameterList &p)
virtual void preEvaluate(typename Traits::PreEvalData d)=0
panzer::Traits::Jacobian::ScalarT ScalarT
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral...
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)