43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 46 #include "Teuchos_RCP.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>
70 const Teuchos::ParameterList& p)
71 : globalIndexer_(indexer)
72 , globalDataKey_(
"Residual Scatter Container")
74 std::string scatterName = p.get<std::string>(
"Scatter Name");
76 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
79 const std::vector<std::string>& names =
80 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
83 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
86 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
88 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
89 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
90 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional ;
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) {
99 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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) {
109 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
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) {
152 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
153 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
155 dirichletCounter_ = Teuchos::null;
159 Teuchos::RCP<LOC> tpetraContainer
160 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
162 dirichletCounter_ = tpetraContainer->
get_f();
163 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
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;
179 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
180 tpetraContainer_->get_f() :
181 tpetraContainer_->get_x();
183 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
184 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
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>
264 const Teuchos::ParameterList& p)
265 : globalIndexer_(indexer)
266 , globalDataKey_(
"Residual Scatter Container")
268 std::string scatterName = p.get<std::string>(
"Scatter Name");
270 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
273 const std::vector<std::string>& names =
274 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
277 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
280 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
282 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
283 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
284 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional ;
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) {
293 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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) {
303 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
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) {
346 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
347 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
349 dirichletCounter_ = Teuchos::null;
353 Teuchos::RCP<LOC> tpetraContainer
354 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
356 dirichletCounter_ = tpetraContainer->
get_f();
357 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
361 using Teuchos::rcp_dynamic_cast;
364 std::vector<std::string> activeParameters =
368 TEUCHOS_ASSERT(!scatterIC_);
369 dfdp_vectors_.clear();
370 for(std::size_t i=0;i<activeParameters.size();i++) {
371 RCP<typename LOC::VectorType> vec =
372 rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
373 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
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;
390 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
391 tpetraContainer_->get_f() :
392 tpetraContainer_->get_x();
394 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
395 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
404 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
405 std::size_t cellLocalId = localCellIds[worksetCellIndex];
407 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
410 LIDs.resize(GIDs.size());
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>
491 const Teuchos::ParameterList& p)
492 : globalIndexer_(indexer)
493 , globalDataKey_(
"Residual Scatter Container")
495 std::string scatterName = p.get<std::string>(
"Scatter Name");
497 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
500 const std::vector<std::string>& names =
501 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
504 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
506 Teuchos::RCP<PHX::DataLayout> dl =
507 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
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) {
515 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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) {
525 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
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) {
568 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
569 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
571 dirichletCounter_ = Teuchos::null;
575 Teuchos::RCP<LOC> tpetraContainer
576 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
578 dirichletCounter_ = tpetraContainer->
get_f();
579 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
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;
594 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
595 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
597 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
598 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
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);
665 TEUCHOS_ASSERT(jacRow.size()==GIDs.size());
667 Jac->replaceGlobalValues(gid, GIDs, jacRow);
panzer::Traits::Tangent::ScalarT ScalarT
Pushes residual values into the residual vector for a Newton-based solve.
const Teuchos::RCP< VectorType > get_f() const
panzer::Traits::Jacobian::ScalarT ScalarT