43 #ifndef PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP 44 #define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP 46 #include "Teuchos_RCP.hpp" 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 58 #include "Thyra_SpmdVectorBase.hpp" 59 #include "Thyra_ProductVectorBase.hpp" 60 #include "Thyra_BlockedLinearOpBase.hpp" 62 #include "Phalanx_DataLayout_MDALayout.hpp" 64 #include "Teuchos_FancyOStream.hpp" 66 template <
typename EvalT,
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
69 const Teuchos::ParameterList& p)
71 std::string scatterName = p.get<std::string>(
"Scatter Name");
72 Teuchos::RCP<PHX::FieldTag> scatterHolder =
73 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
76 const std::vector<std::string>& names =
77 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
79 Teuchos::RCP<PHX::DataLayout> dl =
80 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
83 for (std::size_t eq = 0; eq < names.size(); ++eq) {
84 PHX::MDField<const ScalarT,Cell,NODE>
field = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
87 this->addDependentField(
field.fieldTag());
91 this->addEvaluatedField(*scatterHolder);
93 this->setName(scatterName+
" Scatter Residual");
100 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
103 const Teuchos::ParameterList& p)
104 : globalIndexer_(indexer)
105 , globalDataKey_(
"Residual Scatter Container")
107 std::string scatterName = p.get<std::string>(
"Scatter Name");
109 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
112 const std::vector<std::string>& names =
113 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
116 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
118 Teuchos::RCP<PHX::DataLayout> dl =
119 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
123 for (std::size_t eq = 0; eq < names.size(); ++eq) {
124 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
131 this->addEvaluatedField(*scatterHolder_);
133 if (p.isType<std::string>(
"Global Data Key"))
134 globalDataKey_ = p.get<std::string>(
"Global Data Key");
136 this->setName(scatterName+
" Scatter Residual");
140 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
150 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
151 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
159 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
164 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc.getDataObject(globalDataKey_),
true);
168 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
173 using Teuchos::ArrayRCP;
174 using Teuchos::ptrFromRef;
175 using Teuchos::rcp_dynamic_cast;
178 using Thyra::SpmdVectorBase;
181 std::vector<std::pair<int,GO> > GIDs;
182 std::vector<LO> LIDs;
185 std::string blockId = this->wda(workset).block_id;
186 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
188 RCP<const ContainerType> blockedContainer = blockedContainer_;
190 RCP<ProductVectorBase<double> > r = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true);
198 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
199 std::size_t cellLocalId = localCellIds[worksetCellIndex];
201 globalIndexer_->getElementGIDs(cellLocalId,GIDs,blockId);
204 LIDs.resize(GIDs.size());
205 for(std::size_t i=0;i<GIDs.size();i++) {
207 RCP<const MapType> r_map = blockedContainer->getMapForBlock(GIDs[i].first);
209 LIDs[i] = r_map->getLocalElement(GIDs[i].second);
213 Teuchos::ArrayRCP<double> local_r;
214 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
215 int fieldNum = fieldIds_[fieldIndex];
216 int indexerId = globalIndexer_->getFieldBlock(fieldNum);
219 RCP<SpmdVectorBase<double> > block_r = rcp_dynamic_cast<SpmdVectorBase<double> >(r->getNonconstVectorBlock(indexerId));
220 block_r->getNonconstLocalData(ptrFromRef(local_r));
222 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
238 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
241 const Teuchos::ParameterList& p)
242 : globalIndexer_(indexer)
243 , globalDataKey_(
"Residual Scatter Container")
245 std::string scatterName = p.get<std::string>(
"Scatter Name");
247 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
250 const std::vector<std::string>& names =
251 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
254 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
256 Teuchos::RCP<PHX::DataLayout> dl =
257 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
261 for (std::size_t eq = 0; eq < names.size(); ++eq) {
262 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
269 this->addEvaluatedField(*scatterHolder_);
271 if (p.isType<std::string>(
"Global Data Key"))
272 globalDataKey_ = p.get<std::string>(
"Global Data Key");
274 this->setName(scatterName+
" Scatter Residual (Jacobian)");
278 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
288 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
289 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
297 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
302 using Teuchos::rcp_dynamic_cast;
305 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(d.gedc.getDataObject(globalDataKey_));
307 if(blockedContainer_==Teuchos::null) {
308 RCP<const LOCPair_GlobalEvaluationData> gdata = rcp_dynamic_cast<
const LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true);
309 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(gdata->getGhostedLOC());
314 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
319 using Teuchos::ArrayRCP;
320 using Teuchos::ptrFromRef;
321 using Teuchos::rcp_dynamic_cast;
324 using Thyra::SpmdVectorBase;
328 std::vector<std::pair<int,GO> > GIDs;
329 std::vector<LO> LIDs;
330 std::vector<double> jacRow;
333 std::string blockId = this->wda(workset).block_id;
334 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
336 RCP<const ContainerType> blockedContainer = blockedContainer_;
338 RCP<ProductVectorBase<double> > r = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f());
339 Teuchos::RCP<BlockedLinearOpBase<double> > Jac = rcp_dynamic_cast<BlockedLinearOpBase<double> >(blockedContainer->get_A());
341 int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
342 std::vector<int> blockOffsets(numFieldBlocks+1);
343 for(
int blk=0;blk<numFieldBlocks;blk++) {
344 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
345 blockOffsets[blk] = blockOffset;
348 std::unordered_map<std::pair<int,int>,Teuchos::RCP<CrsMatrixType>,
panzer::pair_hash> jacTpetraBlocks;
356 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
357 std::size_t cellLocalId = localCellIds[worksetCellIndex];
359 globalIndexer_->getElementGIDs(cellLocalId,GIDs,blockId);
362 LIDs.resize(GIDs.size());
363 for(std::size_t i=0;i<GIDs.size();i++) {
365 RCP<const MapType> r_map = blockedContainer->getMapForBlock(GIDs[i].first);
367 LIDs[i] = r_map->getLocalElement(GIDs[i].second);
371 Teuchos::ArrayRCP<double> local_r;
372 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
373 int fieldNum = fieldIds_[fieldIndex];
374 int blockRowIndex = globalIndexer_->getFieldBlock(fieldNum);
377 if(r!=Teuchos::null) {
378 RCP<SpmdVectorBase<double> > block_r = rcp_dynamic_cast<SpmdVectorBase<double> >(r->getNonconstVectorBlock(blockRowIndex));
379 block_r->getNonconstLocalData(ptrFromRef(local_r));
382 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
385 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
387 int rowOffset = elmtOffset[rowBasisNum];
388 int r_lid = LIDs[rowOffset];
391 if(local_r!=Teuchos::null)
392 local_r[r_lid] += (scatterField.val());
394 blockOffsets[numFieldBlocks] = scatterField.size();
397 jacRow.resize(scatterField.size());
399 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex) {
400 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
403 for(
int blockColIndex=0;blockColIndex<numFieldBlocks;blockColIndex++) {
404 int start = blockOffsets[blockColIndex];
405 int end = blockOffsets[blockColIndex+1];
411 std::pair<int,int> blockIndex = std::make_pair(blockRowIndex,blockColIndex);
412 Teuchos::RCP<CrsMatrixType> subJac = jacTpetraBlocks[blockIndex];
415 if(subJac==Teuchos::null) {
416 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac->getNonconstBlock(blockIndex.first,blockIndex.second);
419 if(Teuchos::is_null(tOp))
422 Teuchos::RCP<OperatorType> tpetra_Op = rcp_dynamic_cast<
ThyraLinearOp>(tOp)->getTpetraOperator();
424 jacTpetraBlocks[blockIndex] = subJac;
428 subJac->sumIntoLocalValues(r_lid, Teuchos::arrayViewFromVector(LIDs).view(start,end-start),
429 Teuchos::arrayViewFromVector(jacRow).view(start,end-start));
Pushes residual values into the residual vector for a Newton-based solve.
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
ScatterResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &indexer)
PHX::MDField< const ScalarT, Cell, IP > field
void evaluateFields(typename TRAITS::EvalData d)
Thyra::TpetraLinearOp< RealType, LO, GO, NodeT > ThyraLinearOp
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
Tpetra::CrsMatrix< RealType, LO, GO, NodeT > CrsMatrixType
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
panzer::Traits::Jacobian::ScalarT ScalarT