43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP
44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_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 "Teuchos_FancyOStream.hpp"
71 template<
typename TRAITS,
typename LO,
typename GO>
76 : globalIndexer_(indexer)
77 , globalDataKey_(
"Residual Scatter Container")
79 std::string scatterName = p.
get<std::string>(
"Scatter Name");
84 const std::vector<std::string>& names =
91 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
97 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
98 local_side_id_ = p.
get<
int>(
"Local Side ID");
102 scatterFields_.resize(names.size());
103 for (std::size_t eq = 0; eq < names.size(); ++eq) {
107 this->addDependentField(scatterFields_[eq]);
110 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
112 applyBC_.resize(names.size());
113 for (std::size_t eq = 0; eq < names.size(); ++eq) {
115 this->addDependentField(applyBC_[eq]);
120 this->addEvaluatedField(*scatterHolder_);
122 if (p.
isType<std::string>(
"Global Data Key"))
123 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
125 this->setName(scatterName+
" Scatter Residual");
129 template<
typename TRAITS,
typename LO,
typename GO>
134 fieldIds_.resize(scatterFields_.size());
137 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
139 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
140 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
144 num_nodes = scatterFields_[0].extent(1);
148 template<
typename TRAITS,
typename LO,
typename GO>
155 if(epetraContainer_==Teuchos::null) {
160 dirichletCounter_ = Teuchos::null;
167 dirichletCounter_ = epetraContainer->
get_f();
173 template<
typename TRAITS,
typename LO,
typename GO>
178 std::string blockId = this->wda(workset).block_id;
179 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
183 epetraContainer->
get_f() :
184 epetraContainer->
get_x();
191 auto LIDs = globalIndexer_->getLIDs();
192 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
193 Kokkos::deep_copy(LIDs_h, LIDs);
197 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
198 int fieldNum = fieldIds_[fieldIndex];
199 auto field = PHX::as_view(scatterFields_[fieldIndex]);
200 auto field_h = Kokkos::create_mirror_view(
field);
201 Kokkos::deep_copy(field_h,
field);
203 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
204 std::size_t cellLocalId = localCellIds[worksetCellIndex];
208 const std::pair<std::vector<int>,std::vector<int> > & indicePair
209 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
210 const std::vector<int> & elmtOffset = indicePair.first;
211 const std::vector<int> & basisIdMap = indicePair.second;
214 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
215 int offset = elmtOffset[basis];
216 int lid = LIDs_h(cellLocalId, offset);
220 int basisId = basisIdMap[basis];
223 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
226 (*r)[lid] = field_h(worksetCellIndex,basisId);
229 if(dirichletCounter_!=Teuchos::null)
230 (*dirichletCounter_)[lid] = 1.0;
234 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
237 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
238 int offset = elmtOffset[basis];
239 int lid = LIDs_h(cellLocalId, offset);
243 (*r)[lid] = field_h(worksetCellIndex,basis);
246 if(dirichletCounter_!=Teuchos::null)
247 (*dirichletCounter_)[lid] = 1.0;
259 template<
typename TRAITS,
typename LO,
typename GO>
264 : globalIndexer_(indexer)
265 , globalDataKey_(
"Residual Scatter Container")
267 std::string scatterName = p.
get<std::string>(
"Scatter Name");
272 const std::vector<std::string>& names =
279 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
285 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
286 local_side_id_ = p.
get<
int>(
"Local Side ID");
290 scatterFields_.resize(names.size());
291 for (std::size_t eq = 0; eq < names.size(); ++eq) {
295 this->addDependentField(scatterFields_[eq]);
298 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
300 applyBC_.resize(names.size());
301 for (std::size_t eq = 0; eq < names.size(); ++eq) {
303 this->addDependentField(applyBC_[eq]);
308 this->addEvaluatedField(*scatterHolder_);
310 if (p.
isType<std::string>(
"Global Data Key"))
311 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
313 this->setName(scatterName+
" Scatter Tangent");
317 template<
typename TRAITS,
typename LO,
typename GO>
322 fieldIds_.resize(scatterFields_.size());
325 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
327 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
328 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
332 num_nodes = scatterFields_[0].extent(1);
336 template<
typename TRAITS,
typename LO,
typename GO>
343 if(epetraContainer_==Teuchos::null) {
348 dirichletCounter_ = Teuchos::null;
355 dirichletCounter_ = epetraContainer->
get_f();
360 using Teuchos::rcp_dynamic_cast;
363 std::vector<std::string> activeParameters =
368 dfdp_vectors_.clear();
369 for(std::size_t i=0;i<activeParameters.size();i++) {
371 dfdp_vectors_.push_back(vec);
376 template<
typename TRAITS,
typename LO,
typename GO>
381 std::string blockId = this->wda(workset).block_id;
382 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
386 epetraContainer->
get_f() :
387 epetraContainer->
get_x();
394 auto LIDs = globalIndexer_->getLIDs();
395 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
396 Kokkos::deep_copy(LIDs_h, LIDs);
397 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
398 int fieldNum = fieldIds_[fieldIndex];
399 auto scatterField_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
400 Kokkos::deep_copy(scatterField_h, scatterFields_[fieldIndex].get_static_view());
403 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
404 std::size_t cellLocalId = localCellIds[worksetCellIndex];
408 const std::pair<std::vector<int>,std::vector<int> > & indicePair
409 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
410 const std::vector<int> & elmtOffset = indicePair.first;
411 const std::vector<int> & basisIdMap = indicePair.second;
414 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
415 int offset = elmtOffset[basis];
416 int lid = LIDs_h(cellLocalId, offset);
420 int basisId = basisIdMap[basis];
423 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
426 ScalarT value = scatterField_h(worksetCellIndex,basisId);
431 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
432 (*dfdp_vectors_[d])[lid] = 0.0;
434 for(
int d=0;d<value.size();d++) {
435 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
439 if(dirichletCounter_!=Teuchos::null)
440 (*dirichletCounter_)[lid] = 1.0;
444 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
447 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
448 int offset = elmtOffset[basis];
449 int lid = LIDs_h(cellLocalId, offset);
453 ScalarT value = scatterField_h(worksetCellIndex,basis);
458 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
459 (*dfdp_vectors_[d])[lid] = 0.0;
461 for(
int d=0;d<value.size();d++) {
462 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
466 if(dirichletCounter_!=Teuchos::null)
467 (*dirichletCounter_)[lid] = 1.0;
478 template<
typename TRAITS,
typename LO,
typename GO>
483 : globalIndexer_(indexer)
484 , colGlobalIndexer_(cIndexer)
485 , globalDataKey_(
"Residual Scatter Container")
486 , preserveDiagonal_(false)
488 std::string scatterName = p.
get<std::string>(
"Scatter Name");
493 const std::vector<std::string>& names =
502 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
503 local_side_id_ = p.
get<
int>(
"Local Side ID");
506 scatterFields_.resize(names.size());
507 for (std::size_t eq = 0; eq < names.size(); ++eq) {
511 this->addDependentField(scatterFields_[eq]);
514 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
516 applyBC_.resize(names.size());
517 for (std::size_t eq = 0; eq < names.size(); ++eq) {
519 this->addDependentField(applyBC_[eq]);
524 this->addEvaluatedField(*scatterHolder_);
526 if (p.
isType<std::string>(
"Global Data Key"))
527 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
529 if (p.
isType<
bool>(
"Preserve Diagonal"))
530 preserveDiagonal_ = p.
get<
bool>(
"Preserve Diagonal");
532 this->setName(scatterName+
" Scatter Residual (Jacobian)");
536 template<
typename TRAITS,
typename LO,
typename GO>
541 fieldIds_.resize(scatterFields_.size());
544 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
546 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
547 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
551 num_nodes = scatterFields_[0].extent(1);
552 num_eq = scatterFields_.size();
556 template<
typename TRAITS,
typename LO,
typename GO>
563 if(epetraContainer_==Teuchos::null) {
568 dirichletCounter_ = Teuchos::null;
575 dirichletCounter_ = epetraContainer->
get_f();
581 template<
typename TRAITS,
typename LO,
typename GO>
587 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
589 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
592 std::string blockId = this->wda(workset).block_id;
593 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
606 auto LIDs = globalIndexer_->getLIDs();
607 auto colLIDs = colGlobalIndexer->
getLIDs();
608 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
609 auto colLIDs_h = Kokkos::create_mirror_view(colLIDs);
610 Kokkos::deep_copy(LIDs_h, LIDs);
611 Kokkos::deep_copy(colLIDs_h, colLIDs);
613 std::vector<
typename decltype(scatterFields_[0].get_static_view())::HostMirror> scatterFields_h;
614 for ( std::size_t i=0; i< scatterFields_.size(); ++i) {
615 scatterFields_h.push_back(Kokkos::create_mirror_view(scatterFields_[i].get_static_view()));
616 Kokkos::deep_copy(scatterFields_h[i], scatterFields_[i].get_static_view());
619 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
620 std::size_t cellLocalId = localCellIds[worksetCellIndex];
625 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
626 int fieldNum = fieldIds_[fieldIndex];
629 const std::pair<std::vector<int>,std::vector<int> > & indicePair
630 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
631 const std::vector<int> & elmtOffset = indicePair.first;
632 const std::vector<int> & basisIdMap = indicePair.second;
635 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
636 int offset = elmtOffset[basis];
637 int row = LIDs_h(cellLocalId, offset);
641 int basisId = basisIdMap[basis];
644 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
650 int * rowIndices = 0;
651 double * rowValues = 0;
655 for(
int i=0;i<numEntries;i++) {
656 if(preserveDiagonal_) {
657 if(row!=rowIndices[i])
666 const ScalarT scatterField = (scatterFields_h[fieldIndex])(worksetCellIndex,basisId);
669 (*r)[row] = scatterField.val();
670 if(dirichletCounter_!=Teuchos::null) {
672 (*dirichletCounter_)[row] = 1.0;
676 std::vector<double> jacRow(scatterField.size(),0.0);
678 if(!preserveDiagonal_) {
680 &colLIDs_h(cellLocalId,0));
panzer::Traits::Tangent::ScalarT ScalarT
const Teuchos::RCP< Epetra_Vector > get_x() const
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
virtual int getElementBlockGIDCount(const std::size_t &blockIndex) const =0
How any GIDs are associate with each element in a particular element block.
const Kokkos::View< const panzer::LocalOrdinal **, Kokkos::LayoutRight, PHX::Device > getLIDs() const
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
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)
panzer::Traits::Jacobian::ScalarT ScalarT
const Teuchos::RCP< Epetra_Vector > get_f() const
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)
Pushes residual values into the residual vector for a Newton-based solve.