11 #ifndef PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP
12 #define PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP
15 #include "Teuchos_Assert.hpp"
17 #include "Phalanx_DataLayout.hpp"
19 #include "Epetra_Map.h"
20 #include "Epetra_Vector.h"
21 #include "Epetra_CrsMatrix.h"
31 #include "Phalanx_DataLayout_MDALayout.hpp"
33 #include "Teuchos_FancyOStream.hpp"
39 template<
typename TRAITS,
typename LO,
typename GO>
44 bool useDiscreteAdjoint)
45 : globalIndexer_(indexer)
46 , globalDataKey_(
"Residual Scatter Container")
47 , useDiscreteAdjoint_(useDiscreteAdjoint)
49 std::string scatterName = p.
get<std::string>(
"Scatter Name");
54 const std::vector<std::string>& names =
64 scatterFields_.resize(names.size());
65 for (std::size_t eq = 0; eq < names.size(); ++eq) {
69 this->addDependentField(scatterFields_[eq]);
73 this->addEvaluatedField(*scatterHolder_);
75 if (p.
isType<std::string>(
"Global Data Key"))
76 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
78 this->setName(scatterName+
" Scatter Residual");
82 template<
typename TRAITS,
typename LO,
typename GO>
87 fieldIds_.resize(scatterFields_.size());
89 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
91 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
92 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
97 template<
typename TRAITS,
typename LO,
typename GO>
104 if(epetraContainer_==Teuchos::null) {
112 template<
typename TRAITS,
typename LO,
typename GO>
117 std::string blockId = this->wda(workset).block_id;
118 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
128 auto LIDs = globalIndexer_->getLIDs();
129 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
130 Kokkos::deep_copy(LIDs_h, LIDs);
132 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
133 int fieldNum = fieldIds_[fieldIndex];
134 auto field = PHX::as_view(scatterFields_[fieldIndex]);
135 auto field_h = Kokkos::create_mirror_view(
field);
136 Kokkos::deep_copy(field_h,
field);
138 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
139 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
140 std::size_t cellLocalId = localCellIds[worksetCellIndex];
143 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
144 int offset = elmtOffset[basis];
145 int lid = LIDs_h(cellLocalId, offset);
146 (*r)[lid] += field_h(worksetCellIndex,basis);
156 template<
typename TRAITS,
typename LO,
typename GO>
162 : globalIndexer_(indexer)
164 std::string scatterName = p.
get<std::string>(
"Scatter Name");
169 const std::vector<std::string>& names =
179 scatterFields_.resize(names.size());
180 for (std::size_t eq = 0; eq < names.size(); ++eq) {
184 this->addDependentField(scatterFields_[eq]);
188 this->addEvaluatedField(*scatterHolder_);
190 this->setName(scatterName+
" Scatter Tangent");
194 template<
typename TRAITS,
typename LO,
typename GO>
199 fieldIds_.resize(scatterFields_.size());
201 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
203 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
204 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
209 template<
typename TRAITS,
typename LO,
typename GO>
214 using Teuchos::rcp_dynamic_cast;
217 std::vector<std::string> activeParameters =
220 dfdp_vectors_.clear();
221 for(std::size_t i=0;i<activeParameters.size();i++) {
223 dfdp_vectors_.push_back(vec);
228 template<
typename TRAITS,
typename LO,
typename GO>
233 std::string blockId = this->wda(workset).block_id;
234 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
242 auto LIDs = globalIndexer_->getLIDs();
243 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
244 Kokkos::deep_copy(LIDs_h, LIDs);
245 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
246 int fieldNum = fieldIds_[fieldIndex];
247 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
248 auto scatterField_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
249 Kokkos::deep_copy(scatterField_h, scatterFields_[fieldIndex].get_static_view());
251 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
252 std::size_t cellLocalId = localCellIds[worksetCellIndex];
255 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
256 int offset = elmtOffset[basis];
257 int lid = LIDs_h(cellLocalId, offset);
260 ScalarT value = scatterField_h(worksetCellIndex,basis);
261 for(
int d=0;d<value.size();d++)
262 (*dfdp_vectors_[d])[lid] += value.fastAccessDx(d);
272 template<
typename TRAITS,
typename LO,
typename GO>
277 bool useDiscreteAdjoint)
278 : globalIndexer_(indexer)
279 , colGlobalIndexer_(cIndexer)
280 , globalDataKey_(
"Residual Scatter Container")
281 , useDiscreteAdjoint_(useDiscreteAdjoint)
283 std::string scatterName = p.
get<std::string>(
"Scatter Name");
288 const std::vector<std::string>& names =
298 scatterFields_.resize(names.size());
299 for (std::size_t eq = 0; eq < names.size(); ++eq) {
303 this->addDependentField(scatterFields_[eq]);
307 this->addEvaluatedField(*scatterHolder_);
309 if (p.
isType<std::string>(
"Global Data Key"))
310 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
311 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
312 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
315 if(useDiscreteAdjoint)
318 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
322 template<
typename TRAITS,
typename LO,
typename GO>
329 fieldIds_.resize(scatterFields_.size());
331 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
333 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
334 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
339 template<
typename TRAITS,
typename LO,
typename GO>
346 if(epetraContainer_==Teuchos::null) {
354 template<
typename TRAITS,
typename LO,
typename GO>
358 std::vector<double> jacRow;
360 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
363 std::string blockId = this->wda(workset).block_id;
364 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
370 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
378 auto LIDs = globalIndexer_->
getLIDs();
379 auto colLIDs = colGlobalIndexer->
getLIDs();
380 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
381 auto colLIDs_h = Kokkos::create_mirror_view(colLIDs);
382 Kokkos::deep_copy(LIDs_h, LIDs);
383 Kokkos::deep_copy(colLIDs_h, colLIDs);
384 std::vector<
typename decltype(scatterFields_[0].get_static_view())::host_mirror_type> scatterFields_h;
385 for ( std::size_t i=0; i< scatterFields_.size(); ++i) {
386 scatterFields_h.push_back(Kokkos::create_mirror_view(scatterFields_[i].get_static_view()));
387 Kokkos::deep_copy(scatterFields_h[i], scatterFields_[i].get_static_view());
390 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
391 std::size_t cellLocalId = localCellIds[worksetCellIndex];
393 std::vector<int> cLIDs;
394 for (
int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
395 cLIDs.push_back(colLIDs_h(cellLocalId, i));
397 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
398 for (
int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
399 cLIDs.push_back(colLIDs_h(other_cellLocalId, i));
403 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
404 int fieldNum = fieldIds_[fieldIndex];
405 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
408 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
409 const ScalarT scatterField = (scatterFields_h[fieldIndex])(worksetCellIndex,rowBasisNum);
410 int rowOffset = elmtOffset[rowBasisNum];
411 int row = LIDs_h(cellLocalId,rowOffset);
415 r->SumIntoMyValue(row,0,scatterField.val());
418 if(useDiscreteAdjoint_) {
420 jacRow.resize(scatterField.size());
422 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
423 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
425 for(std::size_t c=0;c<cLIDs.size();c++) {
433 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)