43#ifndef PANZER_SCATTER_RESIDUAL_TPETRA_IMPL_HPP
44#define PANZER_SCATTER_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"
62#include "Tpetra_Vector.hpp"
63#include "Tpetra_Map.hpp"
64#include "Tpetra_CrsMatrix.hpp"
70template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
73 const Teuchos::ParameterList& p)
74 : globalIndexer_(indexer)
75 , globalDataKey_(
"Residual Scatter Container")
77 std::string scatterName = p.get<std::string>(
"Scatter Name");
79 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
82 const std::vector<std::string>& names =
83 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
86 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
88 Teuchos::RCP<PHX::DataLayout> dl =
89 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
92 scatterFields_.resize(names.size());
93 scratch_offsets_.resize(names.size());
94 for (std::size_t eq = 0; eq < names.size(); ++eq) {
95 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
98 this->addDependentField(scatterFields_[eq]);
102 this->addEvaluatedField(*scatterHolder_);
104 if (p.isType<std::string>(
"Global Data Key"))
105 globalDataKey_ = p.get<std::string>(
"Global Data Key");
107 this->setName(scatterName+
" Scatter Residual");
111template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
116 fieldIds_.resize(scatterFields_.size());
117 const Workset & workset_0 = (*d.worksets_)[0];
118 std::string blockId = this->wda(workset_0).
block_id;
122 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
124 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
125 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
127 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldIds_[fd]);
128 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",
offsets.size());
129 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(
offsets.data(),
offsets.size()));
131 scratch_lids_ = PHX::View<LO**>(
"lids",scatterFields_[0].extent(0),
132 globalIndexer_->getElementBlockGIDCount(blockId));
137template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
144 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
146 if(tpetraContainer_==Teuchos::null) {
148 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
149 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
158template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
161 const Teuchos::ParameterList& p)
162 : globalIndexer_(indexer)
163 , globalDataKey_(
"Residual Scatter Container")
165 std::string scatterName = p.get<std::string>(
"Scatter Name");
167 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
170 const std::vector<std::string>& names =
171 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
174 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
176 Teuchos::RCP<PHX::DataLayout> dl =
177 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
180 scatterFields_.resize(names.size());
181 for (std::size_t eq = 0; eq < names.size(); ++eq) {
182 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
185 this->addDependentField(scatterFields_[eq]);
189 this->addEvaluatedField(*scatterHolder_);
191 if (p.isType<std::string>(
"Global Data Key"))
192 globalDataKey_ = p.get<std::string>(
"Global Data Key");
194 this->setName(scatterName+
" Scatter Tangent");
198template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
203 fieldIds_.resize(scatterFields_.size());
205 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
207 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
208 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
213template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
218 using Teuchos::rcp_dynamic_cast;
223 std::vector<std::string> activeParameters =
224 rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject(
"PARAMETER_NAMES"))->getActiveParameters();
226 dfdp_vectors_.
clear();
227 for(std::size_t i=0;i<activeParameters.size();i++) {
228 RCP<typename LOC::VectorType> vec =
229 rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
230 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
231 dfdp_vectors_.push_back(vec_array);
236template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
241 std::string blockId = this->wda(workset).block_id;
242 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
250 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
251 std::size_t cellLocalId = localCellIds[worksetCellIndex];
253 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
256 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
257 int fieldNum = fieldIds_[fieldIndex];
258 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
261 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
262 int offset = elmtOffset[basis];
263 LO lid = LIDs[offset];
264 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
265 for(
int d=0;d<value.size();d++)
266 dfdp_vectors_[d][lid] += value.fastAccessDx(d);
276template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
279 const Teuchos::ParameterList& p)
280 : globalIndexer_(indexer)
281 , globalDataKey_(
"Residual Scatter Container")
282 , my_derivative_size_(0)
283 , other_derivative_size_(0)
285 std::string scatterName = p.get<std::string>(
"Scatter Name");
287 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
290 const std::vector<std::string>& names =
291 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
294 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
296 Teuchos::RCP<PHX::DataLayout> dl =
297 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
300 scatterFields_.resize(names.size());
301 scratch_offsets_.resize(names.size());
302 for (std::size_t eq = 0; eq < names.size(); ++eq) {
303 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
306 this->addDependentField(scatterFields_[eq]);
310 this->addEvaluatedField(*scatterHolder_);
312 if (p.isType<std::string>(
"Global Data Key"))
313 globalDataKey_ = p.get<std::string>(
"Global Data Key");
315 this->setName(scatterName+
" Scatter Residual (Jacobian)");
319template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
324 fieldIds_.resize(scatterFields_.size());
326 const Workset & workset_0 = (*d.worksets_)[0];
327 std::string blockId = this->wda(workset_0).
block_id;
330 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
332 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
333 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
335 int fieldNum = fieldIds_[fd];
336 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
337 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",
offsets.size());
338 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(
offsets.data(),
offsets.size()));
341 my_derivative_size_ = globalIndexer_->getElementBlockGIDCount(blockId);
342 if (Teuchos::nonnull(workset_0.
other)) {
343 auto otherBlockId = workset_0.
other->block_id;
344 other_derivative_size_ = globalIndexer_->getElementBlockGIDCount(otherBlockId);
346 scratch_lids_ = Kokkos::View<LO**, Kokkos::LayoutRight, PHX::Device>(
347 "lids", scatterFields_[0].extent(0), my_derivative_size_ + other_derivative_size_ );
348 scratch_vals_ = Kokkos::View<typename Sacado::ScalarType<ScalarT>::type**, Kokkos::LayoutRight, PHX::Device>(
349 "vals", scatterFields_[0].extent(0), my_derivative_size_ + other_derivative_size_ );
353template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
360 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
362 if(tpetraContainer_==Teuchos::null) {
364 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
365 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
374template <
typename ScalarT,
typename LO,
typename GO,
typename NodeT,
typename LocalMatrixT>
375class ScatterResidual_Jacobian_Functor {
377 typedef typename PHX::Device execution_space;
378 typedef PHX::MDField<const ScalarT,Cell,NODE>
FieldType;
381 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device>
r_data;
384 Kokkos::View<const LO**, Kokkos::LayoutRight, PHX::Device>
lids;
385 Kokkos::View<typename Sacado::ScalarType<ScalarT>::type**, Kokkos::LayoutRight, PHX::Device>
vals;
390 KOKKOS_INLINE_FUNCTION
391 void operator()(
const unsigned int cell)
const
393 int numIds =
lids.extent(1);
396 for(std::size_t basis=0; basis <
offsets.extent(0); basis++) {
397 typename FieldType::array_type::reference_type scatterField =
field(cell,basis);
399 LO lid =
lids(cell,offset);
403 Kokkos::atomic_add(&
r_data(lid,0), scatterField.val());
406 for(
int sensIndex=0;sensIndex<numIds;++sensIndex)
407 vals(cell,sensIndex) = scatterField.fastAccessDx(sensIndex);
410 jac.sumIntoValues(lid, &
lids(cell,0), numIds, &
vals(cell,0),
true,
true);
415template <
typename ScalarT,
typename LO,
typename GO,
typename NodeT>
416class ScatterResidual_Residual_Functor {
418 typedef typename PHX::Device execution_space;
419 typedef PHX::MDField<const ScalarT,Cell,NODE>
FieldType;
421 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device>
r_data;
423 PHX::View<const LO**>
lids;
428 KOKKOS_INLINE_FUNCTION
429 void operator()(
const unsigned int cell)
const
433 for(std::size_t basis=0; basis <
offsets.extent(0); basis++) {
435 LO lid =
lids(cell,offset);
436 Kokkos::atomic_add(&
r_data(lid,0),
field(cell,basis));
446template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
453 std::string blockId = this->wda(workset).block_id;
455 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->
get_f();
457 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
459 ScatterResidual_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
460 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
461 functor.lids = scratch_lids_;
464 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
465 functor.offsets = scratch_offsets_[fieldIndex];
466 functor.field = scatterFields_[fieldIndex];
468 Kokkos::parallel_for(workset.num_cells,functor);
473template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
479 typedef typename LOC::CrsMatrixType::local_matrix_device_type LocalMatrixT;
482 std::string blockId = this->wda(workset).block_id;
484 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->
get_f();
485 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
490 if (Teuchos::nonnull(workset.other)) {
491 auto my_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(0,my_derivative_size_));
492 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,my_scratch_lids);
493 auto other_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(my_derivative_size_,my_derivative_size_ + other_derivative_size_));
494 globalIndexer_->getElementLIDs(workset.other->cell_local_ids_k,other_scratch_lids);
497 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
500 ScatterResidual_Jacobian_Functor<ScalarT,LO,GO,NodeT,LocalMatrixT> functor;
501 functor.fillResidual = (r!=Teuchos::null);
502 if(functor.fillResidual)
503 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
504 functor.jac = Jac->getLocalMatrixDevice();
505 functor.lids = scratch_lids_;
506 functor.vals = scratch_vals_;
509 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
510 functor.offsets = scratch_offsets_[fieldIndex];
511 functor.field = scatterFields_[fieldIndex];
513 Kokkos::parallel_for(workset.num_cells,functor);
PHX::View< const int * > offsets
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.
Kokkos::View< const LO **, Kokkos::LayoutRight, PHX::Device > lids
Kokkos::View< typename Sacado::ScalarType< ScalarT >::type **, Kokkos::LayoutRight, PHX::Device > vals
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > r_data
panzer::Traits::Tangent::ScalarT ScalarT
Pushes residual values into the residual vector for a Newton-based solve.
void clear()
Wipe out stored data.
const Teuchos::RCP< VectorType > get_f() const
std::string block_id
DEPRECATED - use: getElementBlock()
Teuchos::RCP< WorksetDetails > other
FieldType
The type of discretization to use for a field pattern.