43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP
44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP
47 #include "Teuchos_Assert.hpp"
49 #include "Phalanx_DataLayout.hpp"
58 #include "Phalanx_DataLayout_MDALayout.hpp"
60 #include "Teuchos_FancyOStream.hpp"
67 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
71 : globalIndexer_(indexer)
72 , globalDataKey_(
"Residual Scatter Container")
74 std::string scatterName = p.
get<std::string>(
"Scatter Name");
79 const std::vector<std::string>& names =
86 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
92 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
93 local_side_id_ = p.
get<
int>(
"Local Side ID");
97 scatterFields_.resize(names.size());
98 for (std::size_t eq = 0; eq < names.size(); ++eq) {
102 this->addDependentField(scatterFields_[eq]);
105 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
107 applyBC_.resize(names.size());
108 for (std::size_t eq = 0; eq < names.size(); ++eq) {
110 this->addDependentField(applyBC_[eq]);
115 this->addEvaluatedField(*scatterHolder_);
117 if (p.
isType<std::string>(
"Global Data Key"))
118 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
120 this->setName(scatterName+
" Scatter Residual");
124 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
129 fieldIds_.resize(scatterFields_.size());
132 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
134 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
135 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
139 num_nodes = scatterFields_[0].extent(1);
143 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
148 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(globalDataKey_));
150 if(tpetraContainer_==Teuchos::null) {
153 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
155 dirichletCounter_ = Teuchos::null;
160 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
162 dirichletCounter_ = tpetraContainer->get_f();
168 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
172 std::vector<GO> GIDs;
173 std::vector<LO> LIDs;
176 std::string blockId = this->wda(workset).block_id;
177 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
180 tpetraContainer_->get_f() :
181 tpetraContainer_->get_x();
194 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
195 int fieldNum = fieldIds_[fieldIndex];
196 auto scatterFields_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
197 Kokkos::deep_copy(scatterFields_h, scatterFields_[fieldIndex].get_static_view());
200 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
201 std::size_t cellLocalId = localCellIds[worksetCellIndex];
203 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
206 LIDs.resize(GIDs.size());
207 for(std::size_t i=0;i<GIDs.size();i++)
208 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
212 const std::pair<std::vector<int>,std::vector<int> > & indicePair
213 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
214 const std::vector<int> & elmtOffset = indicePair.first;
215 const std::vector<int> & basisIdMap = indicePair.second;
218 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
219 int offset = elmtOffset[basis];
220 LO lid = LIDs[offset];
224 int basisId = basisIdMap[basis];
227 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
230 r_array[lid] = scatterFields_h(worksetCellIndex,basisId);
237 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
240 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
241 int offset = elmtOffset[basis];
242 LO lid = LIDs[offset];
246 r_array[lid] = scatterFields_h(worksetCellIndex,basis);
261 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
265 : globalIndexer_(indexer)
266 , globalDataKey_(
"Residual Scatter Container")
268 std::string scatterName = p.
get<std::string>(
"Scatter Name");
273 const std::vector<std::string>& names =
280 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
286 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
287 local_side_id_ = p.
get<
int>(
"Local Side ID");
291 scatterFields_.resize(names.size());
292 for (std::size_t eq = 0; eq < names.size(); ++eq) {
296 this->addDependentField(scatterFields_[eq]);
299 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
301 applyBC_.resize(names.size());
302 for (std::size_t eq = 0; eq < names.size(); ++eq) {
304 this->addDependentField(applyBC_[eq]);
309 this->addEvaluatedField(*scatterHolder_);
311 if (p.
isType<std::string>(
"Global Data Key"))
312 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
314 this->setName(scatterName+
" Scatter Tangent");
318 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
323 fieldIds_.resize(scatterFields_.size());
326 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
328 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
329 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
333 num_nodes = scatterFields_[0].extent(1);
337 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
342 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(globalDataKey_));
344 if(tpetraContainer_==Teuchos::null) {
347 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
349 dirichletCounter_ = Teuchos::null;
354 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
356 dirichletCounter_ = tpetraContainer->get_f();
361 using Teuchos::rcp_dynamic_cast;
364 std::vector<std::string> activeParameters =
369 dfdp_vectors_.clear();
370 for(std::size_t i=0;i<activeParameters.size();i++) {
372 rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
374 dfdp_vectors_.push_back(vec_array);
379 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
383 std::vector<GO> GIDs;
384 std::vector<LO> LIDs;
387 std::string blockId = this->wda(workset).block_id;
388 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
391 tpetraContainer_->get_f() :
392 tpetraContainer_->get_x();
404 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
405 std::size_t cellLocalId = localCellIds[worksetCellIndex];
407 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
411 for(std::size_t i=0;i<GIDs.size();i++)
412 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
415 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
416 int fieldNum = fieldIds_[fieldIndex];
420 const std::pair<std::vector<int>,std::vector<int> > & indicePair
421 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
422 const std::vector<int> & elmtOffset = indicePair.first;
423 const std::vector<int> & basisIdMap = indicePair.second;
426 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
427 int offset = elmtOffset[basis];
428 LO lid = LIDs[offset];
432 int basisId = basisIdMap[basis];
435 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
438 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
443 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
444 dfdp_vectors_[d][lid] = 0.0;
446 for(
int d=0;d<value.size();d++) {
447 dfdp_vectors_[d][lid] = value.fastAccessDx(d);
455 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
458 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
459 int offset = elmtOffset[basis];
460 LO lid = LIDs[offset];
464 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
469 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
470 dfdp_vectors_[d][lid] = 0.0;
472 for(
int d=0;d<value.size();d++) {
473 dfdp_vectors_[d][lid] = value.fastAccessDx(d);
488 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
492 : globalIndexer_(indexer)
493 , globalDataKey_(
"Residual Scatter Container")
495 std::string scatterName = p.
get<std::string>(
"Scatter Name");
500 const std::vector<std::string>& names =
509 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
510 local_side_id_ = p.
get<
int>(
"Local Side ID");
513 scatterFields_.resize(names.size());
514 for (std::size_t eq = 0; eq < names.size(); ++eq) {
518 this->addDependentField(scatterFields_[eq]);
521 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
523 applyBC_.resize(names.size());
524 for (std::size_t eq = 0; eq < names.size(); ++eq) {
526 this->addDependentField(applyBC_[eq]);
531 this->addEvaluatedField(*scatterHolder_);
533 if (p.
isType<std::string>(
"Global Data Key"))
534 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
536 this->setName(scatterName+
" Scatter Residual (Jacobian)");
540 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
545 fieldIds_.resize(scatterFields_.size());
547 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
549 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
550 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
554 num_nodes = scatterFields_[0].extent(1);
555 num_eq = scatterFields_.size();
559 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
564 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(globalDataKey_));
566 if(tpetraContainer_==Teuchos::null) {
569 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
571 dirichletCounter_ = Teuchos::null;
576 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
578 dirichletCounter_ = tpetraContainer->get_f();
584 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
588 std::vector<GO> GIDs;
591 std::string blockId = this->wda(workset).block_id;
592 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
606 auto LIDs = globalIndexer_->getLIDs();
607 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
608 Kokkos::deep_copy(LIDs_h, LIDs);
610 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
611 int fieldNum = fieldIds_[fieldIndex];
612 auto scatterFields_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
613 Kokkos::deep_copy(scatterFields_h, scatterFields_[fieldIndex].get_static_view());
614 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
615 std::size_t cellLocalId = localCellIds[worksetCellIndex];
617 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
620 const std::pair<std::vector<int>,std::vector<int> > & indicePair
621 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
622 const std::vector<int> & elmtOffset = indicePair.first;
623 const std::vector<int> & basisIdMap = indicePair.second;
626 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
627 int offset = elmtOffset[basis];
628 int lid = LIDs_h(cellLocalId, offset);
632 int basisId = basisIdMap[basis];
635 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
640 std::size_t sz = Jac->getNumEntriesInLocalRow(lid);
641 std::size_t numEntries = 0;
642 typename LOC::CrsMatrixType::nonconst_local_inds_host_view_type rowIndices(
"indices", sz);
643 typename LOC::CrsMatrixType::nonconst_values_host_view_type rowValues(
"values", sz);
646 Jac->getLocalRowCopy(lid,rowIndices,rowValues,numEntries);
648 for(std::size_t i=0;i<numEntries;i++)
651 Jac->replaceLocalValues(lid,rowIndices,rowValues);
654 GO gid = GIDs[offset];
655 const ScalarT scatterField = scatterFields_h(worksetCellIndex,basisId);
657 r_array[lid] = scatterField.val();
661 std::vector<double> jacRow(scatterField.size(),0.0);
663 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
664 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
667 Jac->replaceGlobalValues(gid, GIDs, jacRow);
panzer::Traits::Tangent::ScalarT ScalarT
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void resize(const size_type n, const T &val=T())
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pushes residual values into the residual vector for a Newton-based solve.
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
panzer::Traits::Jacobian::ScalarT ScalarT