43 #ifndef PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP
44 #define PANZER_SCATTER_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"
63 #include "Phalanx_DataLayout_MDALayout.hpp"
65 #include "Teuchos_FancyOStream.hpp"
71 template<
typename TRAITS,
typename LO,
typename GO>
76 bool useDiscreteAdjoint)
77 : globalIndexer_(indexer)
78 , globalDataKey_(
"Residual Scatter Container")
79 , useDiscreteAdjoint_(useDiscreteAdjoint)
81 std::string scatterName = p.
get<std::string>(
"Scatter Name");
86 const std::vector<std::string>& names =
96 scatterFields_.resize(names.size());
97 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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) {
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;
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);
188 template<
typename TRAITS,
typename LO,
typename GO>
194 : globalIndexer_(indexer)
196 std::string scatterName = p.
get<std::string>(
"Scatter Name");
201 const std::vector<std::string>& names =
211 scatterFields_.resize(names.size());
212 for (std::size_t eq = 0; eq < names.size(); ++eq) {
216 this->addDependentField(scatterFields_[eq]);
220 this->addEvaluatedField(*scatterHolder_);
222 this->setName(scatterName+
" Scatter Tangent");
226 template<
typename TRAITS,
typename LO,
typename GO>
231 fieldIds_.resize(scatterFields_.size());
233 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
235 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
236 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
241 template<
typename TRAITS,
typename LO,
typename GO>
246 using Teuchos::rcp_dynamic_cast;
249 std::vector<std::string> activeParameters =
252 dfdp_vectors_.clear();
253 for(std::size_t i=0;i<activeParameters.size();i++) {
255 dfdp_vectors_.push_back(vec);
260 template<
typename TRAITS,
typename LO,
typename GO>
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);
305 template<
typename TRAITS,
typename LO,
typename GO>
310 bool useDiscreteAdjoint)
311 : globalIndexer_(indexer)
312 , colGlobalIndexer_(cIndexer)
313 , globalDataKey_(
"Residual Scatter Container")
314 , useDiscreteAdjoint_(useDiscreteAdjoint)
316 std::string scatterName = p.
get<std::string>(
"Scatter Name");
321 const std::vector<std::string>& names =
331 scatterFields_.resize(names.size());
332 for (std::size_t eq = 0; eq < names.size(); ++eq) {
336 this->addDependentField(scatterFields_[eq]);
340 this->addEvaluatedField(*scatterHolder_);
342 if (p.
isType<std::string>(
"Global Data Key"))
343 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
344 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
345 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
348 if(useDiscreteAdjoint)
351 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
355 template<
typename TRAITS,
typename LO,
typename GO>
362 fieldIds_.resize(scatterFields_.size());
364 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
366 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
367 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
372 template<
typename TRAITS,
typename LO,
typename GO>
379 if(epetraContainer_==Teuchos::null) {
387 template<
typename TRAITS,
typename LO,
typename GO>
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;
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));
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++) {
466 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
const Kokkos::View< const panzer::LocalOrdinal **, Kokkos::LayoutRight, PHX::Device > getLIDs() const
panzer::Traits::Tangent::ScalarT ScalarT
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T * ptrFromStlVector(std::vector< T > &v)
Pushes residual values into the residual vector for a Newton-based solve.
panzer::Traits::Jacobian::ScalarT ScalarT
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
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)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)