43 #ifndef PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP 44 #define PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP 46 #include "Teuchos_RCP.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" 63 #include "Phalanx_DataLayout_MDALayout.hpp" 65 #include "Teuchos_FancyOStream.hpp" 71 template<
typename TRAITS,
typename LO,
typename GO>
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");
114 template<
typename TRAITS,
typename LO,
typename GO>
119 fieldIds_.resize(scatterFields_.size());
121 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
123 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
124 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
129 template<
typename TRAITS,
typename LO,
typename GO>
136 if(epetraContainer_==Teuchos::null) {
138 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
144 template<
typename TRAITS,
typename LO,
typename GO>
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 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
161 std::size_t cellLocalId = localCellIds[worksetCellIndex];
163 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
166 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
167 int fieldNum = fieldIds_[fieldIndex];
168 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
171 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
172 int offset = elmtOffset[basis];
173 int lid = LIDs[offset];
174 (*r)[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis);
184 template<
typename TRAITS,
typename LO,
typename GO>
188 const Teuchos::ParameterList& p,
190 : globalIndexer_(indexer)
192 std::string scatterName = p.get<std::string>(
"Scatter Name");
194 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
197 const std::vector<std::string>& names =
198 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
201 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
203 Teuchos::RCP<PHX::DataLayout> dl =
204 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
207 scatterFields_.resize(names.size());
208 for (std::size_t eq = 0; eq < names.size(); ++eq) {
209 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
212 this->addDependentField(scatterFields_[eq]);
216 this->addEvaluatedField(*scatterHolder_);
218 this->setName(scatterName+
" Scatter Tangent");
222 template<
typename TRAITS,
typename LO,
typename GO>
227 fieldIds_.resize(scatterFields_.size());
229 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
231 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
232 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
237 template<
typename TRAITS,
typename LO,
typename GO>
242 using Teuchos::rcp_dynamic_cast;
245 std::vector<std::string> activeParameters =
248 dfdp_vectors_.clear();
249 for(std::size_t i=0;i<activeParameters.size();i++) {
250 RCP<Epetra_Vector> vec = rcp_dynamic_cast<
EpetraLinearObjContainer>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
251 dfdp_vectors_.push_back(vec);
256 template<
typename TRAITS,
typename LO,
typename GO>
261 std::string blockId = this->wda(workset).block_id;
262 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
270 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
271 std::size_t cellLocalId = localCellIds[worksetCellIndex];
273 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
276 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
277 int fieldNum = fieldIds_[fieldIndex];
278 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
281 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
282 int offset = elmtOffset[basis];
283 int lid = LIDs[offset];
286 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
287 for(
int d=0;d<value.size();d++)
288 (*dfdp_vectors_[d])[lid] += value.fastAccessDx(d);
298 template<
typename TRAITS,
typename LO,
typename GO>
302 const Teuchos::ParameterList& p,
303 bool useDiscreteAdjoint)
304 : globalIndexer_(indexer)
305 , colGlobalIndexer_(cIndexer)
306 , globalDataKey_(
"Residual Scatter Container")
307 , useDiscreteAdjoint_(useDiscreteAdjoint)
309 std::string scatterName = p.get<std::string>(
"Scatter Name");
311 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
314 const std::vector<std::string>& names =
315 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
318 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
320 Teuchos::RCP<PHX::DataLayout> dl =
321 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
324 scatterFields_.resize(names.size());
325 for (std::size_t eq = 0; eq < names.size(); ++eq) {
326 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
329 this->addDependentField(scatterFields_[eq]);
333 this->addEvaluatedField(*scatterHolder_);
335 if (p.isType<std::string>(
"Global Data Key"))
336 globalDataKey_ = p.get<std::string>(
"Global Data Key");
337 if (p.isType<
bool>(
"Use Discrete Adjoint"))
338 useDiscreteAdjoint = p.get<
bool>(
"Use Discrete Adjoint");
341 if(useDiscreteAdjoint)
342 { TEUCHOS_ASSERT(colGlobalIndexer_==globalIndexer_); }
344 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
348 template<
typename TRAITS,
typename LO,
typename GO>
355 fieldIds_.resize(scatterFields_.size());
357 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
359 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
360 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
365 template<
typename TRAITS,
typename LO,
typename GO>
372 if(epetraContainer_==Teuchos::null) {
374 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
380 template<
typename TRAITS,
typename LO,
typename GO>
384 std::vector<double> jacRow;
386 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
389 std::string blockId = this->wda(workset).block_id;
390 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
392 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
393 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
395 const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> >&
396 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
404 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
405 std::size_t cellLocalId = localCellIds[worksetCellIndex];
407 auto rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
408 auto initial_cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
409 std::vector<int> cLIDs;
410 for (
int i(0); i < static_cast<int>(initial_cLIDs.extent(0)); ++i)
411 cLIDs.push_back(initial_cLIDs(i));
412 if (Teuchos::nonnull(workset.other)) {
413 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
414 auto other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
415 for (
int i(0); i < static_cast<int>(other_cLIDs.extent(0)); ++i)
416 cLIDs.push_back(other_cLIDs(i));
420 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
421 int fieldNum = fieldIds_[fieldIndex];
422 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
425 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
426 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
427 int rowOffset = elmtOffset[rowBasisNum];
428 int row = rLIDs[rowOffset];
432 r->SumIntoMyValue(row,0,scatterField.val());
435 if(useDiscreteAdjoint_) {
437 jacRow.resize(scatterField.size());
439 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
440 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
442 for(std::size_t c=0;c<cLIDs.size();c++) {
443 int err = Jac->SumIntoMyValues(cLIDs[c], 1, &jacRow[c],&row);
444 TEUCHOS_ASSERT_EQUALITY(err,0);
448 int err = Jac->SumIntoMyValues(
450 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
453 TEUCHOS_ASSERT_EQUALITY(err,0);
panzer::Traits::Tangent::ScalarT ScalarT
T * ptrFromStlVector(std::vector< T > &v)
Pushes residual values into the residual vector for a Newton-based solve.
panzer::Traits::Jacobian::ScalarT ScalarT