74 const Teuchos::RCP<const panzer::GlobalIndexer> & ,
75 const Teuchos::ParameterList& p,
76 bool useDiscreteAdjoint)
77 : globalIndexer_(indexer)
78 , globalDataKey_(
"Residual Scatter Container")
79 , useDiscreteAdjoint_(useDiscreteAdjoint)
81 std::string scatterName = p.get<std::string>(
"Scatter Name");
83 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
86 const std::vector<std::string>& names =
87 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
90 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
92 Teuchos::RCP<PHX::DataLayout> dl =
93 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
96 scatterFields_.resize(names.size());
97 for (std::size_t eq = 0; eq < names.size(); ++eq) {
98 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
101 this->addDependentField(scatterFields_[eq]);
105 this->addEvaluatedField(*scatterHolder_);
107 if (p.isType<std::string>(
"Global Data Key"))
108 globalDataKey_ = p.get<std::string>(
"Global Data Key");
110 this->setName(scatterName+
" Scatter Residual");
149 std::string blockId = this->wda(workset).block_id;
150 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
152 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
160 auto LIDs = globalIndexer_->getLIDs();
161 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
162 Kokkos::deep_copy(LIDs_h, LIDs);
164 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
165 int fieldNum = fieldIds_[fieldIndex];
166 auto field = PHX::as_view(scatterFields_[fieldIndex]);
167 auto field_h = Kokkos::create_mirror_view(
field);
168 Kokkos::deep_copy(field_h,
field);
170 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
171 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
172 std::size_t cellLocalId = localCellIds[worksetCellIndex];
175 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
176 int offset = elmtOffset[basis];
177 int lid = LIDs_h(cellLocalId, offset);
178 (*r)[lid] += field_h(worksetCellIndex,basis);
191 const Teuchos::RCP<const panzer::GlobalIndexer> & ,
192 const Teuchos::ParameterList& p,
194 : globalIndexer_(indexer)
196 std::string scatterName = p.get<std::string>(
"Scatter Name");
198 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
201 const std::vector<std::string>& names =
202 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
205 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
207 Teuchos::RCP<PHX::DataLayout> dl =
208 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
211 scatterFields_.resize(names.size());
212 for (std::size_t eq = 0; eq < names.size(); ++eq) {
213 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
216 this->addDependentField(scatterFields_[eq]);
220 this->addEvaluatedField(*scatterHolder_);
222 this->setName(scatterName+
" Scatter Tangent");
265 std::string blockId = this->wda(workset).block_id;
266 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
274 auto LIDs = globalIndexer_->getLIDs();
275 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
276 Kokkos::deep_copy(LIDs_h, LIDs);
277 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
278 int fieldNum = fieldIds_[fieldIndex];
279 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
280 auto scatterField_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
281 Kokkos::deep_copy(scatterField_h, scatterFields_[fieldIndex].get_static_view());
283 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
284 std::size_t cellLocalId = localCellIds[worksetCellIndex];
286 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
288 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
289 int offset = elmtOffset[basis];
290 int lid = LIDs_h(cellLocalId, offset);
293 ScalarT value = scatterField_h(worksetCellIndex,basis);
294 for(
int d=0;d<value.size();d++)
295 (*dfdp_vectors_[d])[lid] += value.fastAccessDx(d);
391 std::vector<double> jacRow;
393 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
396 std::string blockId = this->wda(workset).block_id;
397 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
399 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
400 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
402 const Teuchos::RCP<const panzer::GlobalIndexer>&
403 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
411 auto LIDs = globalIndexer_->getLIDs();
412 auto colLIDs = colGlobalIndexer->getLIDs();
413 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
414 auto colLIDs_h = Kokkos::create_mirror_view(colLIDs);
415 Kokkos::deep_copy(LIDs_h, LIDs);
416 Kokkos::deep_copy(colLIDs_h, colLIDs);
417 std::vector<
typename decltype(scatterFields_[0].get_static_view())::HostMirror> scatterFields_h;
418 for ( std::size_t i=0; i< scatterFields_.size(); ++i) {
419 scatterFields_h.push_back(Kokkos::create_mirror_view(scatterFields_[i].get_static_view()));
420 Kokkos::deep_copy(scatterFields_h[i], scatterFields_[i].get_static_view());
423 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
424 std::size_t cellLocalId = localCellIds[worksetCellIndex];
426 std::vector<int> cLIDs;
427 for (
int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
428 cLIDs.push_back(colLIDs_h(cellLocalId, i));
429 if (Teuchos::nonnull(workset.other)) {
430 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
431 for (
int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
432 cLIDs.push_back(colLIDs_h(other_cellLocalId, i));
436 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
437 int fieldNum = fieldIds_[fieldIndex];
438 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
441 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
442 const ScalarT scatterField = (scatterFields_h[fieldIndex])(worksetCellIndex,rowBasisNum);
443 int rowOffset = elmtOffset[rowBasisNum];
444 int row = LIDs_h(cellLocalId,rowOffset);
448 r->SumIntoMyValue(row,0,scatterField.val());
451 if(useDiscreteAdjoint_) {
453 jacRow.resize(scatterField.size());
455 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
456 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
458 for(std::size_t c=0;c<cLIDs.size();c++) {
459 int err = Jac->SumIntoMyValues(cLIDs[c], 1, &jacRow[c],&row);
460 TEUCHOS_ASSERT_EQUALITY(err,0);
464 int err = Jac->SumIntoMyValues(
466 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
470 TEUCHOS_ASSERT_EQUALITY(err,0);