Panzer  Version of the Day
Panzer_ScatterResidual_BlockedTpetra_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP
44 #define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP
45 
46 #include "Teuchos_RCP.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 #include "Phalanx_DataLayout.hpp"
50 
51 #include "Panzer_GlobalIndexer.hpp"
53 #include "Panzer_PureBasis.hpp"
56 #include "Panzer_HashUtils.hpp"
58 
59 #include "Thyra_ProductVectorBase.hpp"
60 #include "Thyra_BlockedLinearOpBase.hpp"
61 #include "Thyra_TpetraVector.hpp"
62 #include "Thyra_TpetraLinearOp.hpp"
63 #include "Tpetra_CrsMatrix.hpp"
64 #include "KokkosSparse_CrsMatrix.hpp"
65 
66 #include "Phalanx_DataLayout_MDALayout.hpp"
67 
68 #include "Teuchos_FancyOStream.hpp"
69 
70 template <typename EvalT,typename TRAITS,typename LO,typename GO,typename NodeT>
72 ScatterResidual_BlockedTpetra(const Teuchos::RCP<const BlockedDOFManager> & /* indexer */,
73  const Teuchos::ParameterList& p)
74 {
75  std::string scatterName = p.get<std::string>("Scatter Name");
76  Teuchos::RCP<PHX::FieldTag> scatterHolder =
77  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
78 
79  // get names to be evaluated
80  const std::vector<std::string>& names =
81  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
82 
83  Teuchos::RCP<PHX::DataLayout> dl =
84  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
85 
86  // build the vector of fields that this is dependent on
87  for (std::size_t eq = 0; eq < names.size(); ++eq) {
88  PHX::MDField<const ScalarT,Cell,NODE> field = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
89 
90  // tell the field manager that we depend on this field
91  this->addDependentField(field.fieldTag());
92  }
93 
94  // this is what this evaluator provides
95  this->addEvaluatedField(*scatterHolder);
96 
97  this->setName(scatterName+" Scatter Residual");
98 }
99 
100 // **********************************************************************
101 // Specialization: Residual
102 // **********************************************************************
103 
104 template <typename TRAITS,typename LO,typename GO,typename NodeT>
106 ScatterResidual_BlockedTpetra(const Teuchos::RCP<const BlockedDOFManager> & indexer,
107  const Teuchos::ParameterList& p)
108  : globalIndexer_(indexer)
109  , globalDataKey_("Residual Scatter Container")
110 {
111  std::string scatterName = p.get<std::string>("Scatter Name");
112  scatterHolder_ =
113  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
114 
115  // get names to be evaluated
116  const std::vector<std::string>& names =
117  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
118 
119  // grab map from evaluated names to field names
120  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
121 
122  Teuchos::RCP<PHX::DataLayout> dl =
123  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
124 
125  // build the vector of fields that this is dependent on
126  scatterFields_.resize(names.size());
127  for (std::size_t eq = 0; eq < names.size(); ++eq) {
128  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
129 
130  // tell the field manager that we depend on this field
131  this->addDependentField(scatterFields_[eq]);
132  }
133 
134  // this is what this evaluator provides
135  this->addEvaluatedField(*scatterHolder_);
136 
137  if (p.isType<std::string>("Global Data Key"))
138  globalDataKey_ = p.get<std::string>("Global Data Key");
139 
140  this->setName(scatterName+" Scatter Residual");
141 }
142 
143 // **********************************************************************
144 template <typename TRAITS,typename LO,typename GO,typename NodeT>
146 postRegistrationSetup(typename TRAITS::SetupData d,
147  PHX::FieldManager<TRAITS>& /* fm */)
148 {
149  const Workset & workset_0 = (*d.worksets_)[0];
150  const std::string blockId = this->wda(workset_0).block_id;
151 
152  fieldIds_.resize(scatterFields_.size());
153  fieldOffsets_.resize(scatterFields_.size());
154  fieldGlobalIndexers_.resize(scatterFields_.size());
155  productVectorBlockIndex_.resize(scatterFields_.size());
156  int maxElementBlockGIDCount = -1;
157  for(std::size_t fd=0; fd < scatterFields_.size(); ++fd) {
158  const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
159  const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
160  productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
161  fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
162  fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName); // Field number in the sub-global-indexer
163 
164  const std::vector<int>& offsets = fieldGlobalIndexers_[fd]->getGIDFieldOffsets(blockId,fieldIds_[fd]);
165  fieldOffsets_[fd] = PHX::View<int*>("ScatterResidual_BlockedTpetra(Residual):fieldOffsets",offsets.size());
166  auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
167  for (std::size_t i=0; i < offsets.size(); ++i)
168  hostOffsets(i) = offsets[i];
169  Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
170 
171  maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
172  }
173 
174  // We will use one workset lid view for all fields, but has to be
175  // sized big enough to hold the largest elementBlockGIDCount in the
176  // ProductVector.
177  worksetLIDs_ = PHX::View<LO**>("ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
178  scatterFields_[0].extent(0),
179  maxElementBlockGIDCount);
180 }
181 
182 // **********************************************************************
183 template <typename TRAITS,typename LO,typename GO,typename NodeT>
185 preEvaluate(typename TRAITS::PreEvalData d)
186 {
187  // extract linear object container
188  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
189 }
190 
191 // **********************************************************************
192 template <typename TRAITS,typename LO,typename GO,typename NodeT>
194 evaluateFields(typename TRAITS::EvalData workset)
195 {
196  using Teuchos::RCP;
197  using Teuchos::rcp_dynamic_cast;
198  using Thyra::VectorBase;
200 
201  const auto& localCellIds = this->wda(workset).cell_local_ids_k;
202  const RCP<ProductVectorBase<double>> thyraBlockResidual = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f(),true);
203 
204  // Loop over scattered fields
205  int currentWorksetLIDSubBlock = -1;
206  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
207  // workset LIDs only change for different sub blocks
208  if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
209  fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
210  currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
211  }
212 
213  auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),true))->getTpetraVector());
214  const auto& kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::OverwriteAll);
215 
216  // Class data fields for lambda capture
217  const auto& fieldOffsets = fieldOffsets_[fieldIndex];
218  const auto& worksetLIDs = worksetLIDs_;
219  const auto& fieldValues = scatterFields_[fieldIndex].get_static_view();
220 
221  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
222  for(int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
223  const int lid = worksetLIDs(cell,fieldOffsets(basis));
224  Kokkos::atomic_add(&kokkosResidual(lid,0), fieldValues(cell,basis));
225  }
226  });
227  }
228 }
229 
230 // **********************************************************************
231 // Specialization: Jacobian
232 // **********************************************************************
233 
234 template <typename TRAITS,typename LO,typename GO,typename NodeT>
236 ScatterResidual_BlockedTpetra(const Teuchos::RCP<const BlockedDOFManager> & indexer,
237  const Teuchos::ParameterList& p)
238  : globalIndexer_(indexer)
239  , globalDataKey_("Residual Scatter Container")
240 {
241  std::string scatterName = p.get<std::string>("Scatter Name");
242  scatterHolder_ =
243  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
244 
245  // get names to be evaluated
246  const std::vector<std::string>& names =
247  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
248 
249  // grab map from evaluated names to field names
250  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
251 
252  Teuchos::RCP<PHX::DataLayout> dl =
253  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
254 
255  // build the vector of fields that this is dependent on
256  scatterFields_.resize(names.size());
257  for (std::size_t eq = 0; eq < names.size(); ++eq) {
258  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
259 
260  // tell the field manager that we depend on this field
261  this->addDependentField(scatterFields_[eq]);
262  }
263 
264  // this is what this evaluator provides
265  this->addEvaluatedField(*scatterHolder_);
266 
267  if (p.isType<std::string>("Global Data Key"))
268  globalDataKey_ = p.get<std::string>("Global Data Key");
269 
270  this->setName(scatterName+" Scatter Residual (Jacobian)");
271 }
272 
273 // **********************************************************************
274 template <typename TRAITS,typename LO,typename GO,typename NodeT>
276 postRegistrationSetup(typename TRAITS::SetupData d,
277  PHX::FieldManager<TRAITS>& /* fm */)
278 {
279  const Workset & workset_0 = (*d.worksets_)[0];
280  const std::string blockId = this->wda(workset_0).block_id;
281 
282  fieldIds_.resize(scatterFields_.size());
283  fieldOffsets_.resize(scatterFields_.size());
284  productVectorBlockIndex_.resize(scatterFields_.size());
285  for (std::size_t fd=0; fd < scatterFields_.size(); ++fd) {
286  const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
287  const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
288  productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
289  const auto& fieldGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
290  fieldIds_[fd] = fieldGlobalIndexer->getFieldNum(fieldName); // Field number in the sub-global-indexer
291 
292  const std::vector<int>& offsets = globalIndexer_->getGIDFieldOffsets(blockId,globalFieldNum);
293  fieldOffsets_[fd] = PHX::View<int*>("ScatterResidual_BlockedTpetra(Jacobian):fieldOffsets",offsets.size());
294  auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
295  for (std::size_t i=0; i < offsets.size(); ++i)
296  hostOffsets(i) = offsets[i];
297  Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
298  }
299 
300  // This is sized differently than the Residual implementation since
301  // we need the LIDs for all sub-blocks, not just the single
302  // sub-block for the field residual scatter.
303  int elementBlockGIDCount = 0;
304  for (const auto blockDOFMgr : globalIndexer_->getFieldDOFManagers())
305  elementBlockGIDCount += blockDOFMgr->getElementBlockGIDCount(blockId);
306 
307  worksetLIDs_ = PHX::View<LO**>("ScatterResidual_BlockedTpetra(Jacobian):worksetLIDs",
308  scatterFields_[0].extent(0),
309  elementBlockGIDCount);
310 
311  // Compute the block offsets
312  const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
313  const int numBlocks = static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
314  blockOffsets_ = PHX::View<LO*>("ScatterResidual_BlockedTpetra(Jacobian):blockOffsets_",
315  numBlocks+1); // Number of fields, plus a sentinel
316  const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
317  for (int blk=0;blk<numBlocks;blk++) {
318  int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
319  hostBlockOffsets(blk) = blockOffset;
320  }
321  hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
322  Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
323 
324  // Make sure the that hard coded derivative dimension in the
325  // evaluate call is large enough to hold all derivatives for each
326  // sub block load
327  for (int blk=0;blk<numBlocks;blk++) {
328  const int blockDerivativeSize = hostBlockOffsets(blk+1) - hostBlockOffsets(blk);
329  TEUCHOS_TEST_FOR_EXCEPTION(blockDerivativeSize > 256, std::runtime_error,
330  "ERROR: the derivative dimension for sub block "
331  << blk << "with a value of " << blockDerivativeSize
332  << "is larger than the size allocated for cLIDs and vals "
333  << "in the evaluate call! You must manually increase the "
334  << "size and recompile!");
335  }
336 }
337 
338 // **********************************************************************
339 template <typename TRAITS,typename LO,typename GO,typename NodeT>
341 preEvaluate(typename TRAITS::PreEvalData d)
342 {
343  using Teuchos::RCP;
344  using Teuchos::rcp_dynamic_cast;
345 
346  // extract linear object container
347  blockedContainer_ = rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_));
348 
349  if(blockedContainer_==Teuchos::null) {
350  RCP<const LOCPair_GlobalEvaluationData> gdata = rcp_dynamic_cast<const LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true);
351  blockedContainer_ = rcp_dynamic_cast<const ContainerType>(gdata->getGhostedLOC());
352  }
353 }
354 
355 // **********************************************************************
356 template <typename TRAITS,typename LO,typename GO,typename NodeT>
358 evaluateFields(typename TRAITS::EvalData workset)
359 {
360  using Teuchos::RCP;
361  using Teuchos::rcp_dynamic_cast;
362  using Thyra::VectorBase;
365 
366  const auto& localCellIds = this->wda(workset).cell_local_ids_k;
367 
368  const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
369  const RCP<const ContainerType> blockedContainer = blockedContainer_;
370  const RCP<ProductVectorBase<double>> thyraBlockResidual = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f());
371  const bool haveResidual = Teuchos::nonnull(thyraBlockResidual);
372  const RCP<BlockedLinearOpBase<double>> Jac = rcp_dynamic_cast<BlockedLinearOpBase<double> >(blockedContainer_->get_A(),true);
373 
374  // Get the local data for the sub-block crs matrices. First allocate
375  // on host and then deep_copy to device. The sub-blocks are
376  // unmanaged since they are allocated and ref counted separately on
377  // host.
378  using LocalMatrixType = KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>, size_t>;
379  typename PHX::View<LocalMatrixType**>::HostMirror
380  hostJacTpetraBlocks("panzer::ScatterResidual_BlockTpetra<Jacobian>::hostJacTpetraBlocks", numFieldBlocks,numFieldBlocks);
381 
382  PHX::View<int**> blockExistsInJac = PHX::View<int**>("blockExistsInJac_",numFieldBlocks,numFieldBlocks);
383  auto hostBlockExistsInJac = Kokkos::create_mirror_view(blockExistsInJac);
384 
385  for (int row=0; row < numFieldBlocks; ++row) {
386  for (int col=0; col < numFieldBlocks; ++col) {
387  const auto thyraTpetraOperator = rcp_dynamic_cast<Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(Jac->getNonconstBlock(row,col),false);
388  if (nonnull(thyraTpetraOperator)) {
389 
390  // HACK to enforce views in the CrsGraph to be
391  // Unmanaged. Passing in the MemoryTrait<Unmanaged> doesn't
392  // work as the CrsGraph in the CrsMatrix ignores the
393  // MemoryTrait. Need to use the runtime constructor by passing
394  // in points to ensure Unmanaged. See:
395  // https://github.com/kokkos/kokkos/issues/1581
396 
397  // These two lines are the original code we can revert to when #1581 is fixed.
398  // const auto crsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),true);
399  // new (&hostJacTpetraBlocks(row,col)) KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>> (crsMatrix->getLocalMatrix());
400 
401  // Instead do this
402  {
403  // Grab the local managed matrix and graph
404  const auto tpetraCrsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),true);
405  const auto managedMatrix = tpetraCrsMatrix->getLocalMatrixDevice();
406  const auto managedGraph = managedMatrix.graph;
407 
408  // Create runtime unmanaged versions
409  using StaticCrsGraphType = typename LocalMatrixType::StaticCrsGraphType;
410  StaticCrsGraphType unmanagedGraph;
411  unmanagedGraph.entries = typename StaticCrsGraphType::entries_type(managedGraph.entries.data(),managedGraph.entries.extent(0));
412  unmanagedGraph.row_map = typename StaticCrsGraphType::row_map_type(managedGraph.row_map.data(),managedGraph.row_map.extent(0));
413  unmanagedGraph.row_block_offsets = typename StaticCrsGraphType::row_block_type(managedGraph.row_block_offsets.data(),managedGraph.row_block_offsets.extent(0));
414 
415  typename LocalMatrixType::values_type unmanagedValues(managedMatrix.values.data(),managedMatrix.values.extent(0));
416  LocalMatrixType unmanagedMatrix(managedMatrix.values.label(), managedMatrix.numCols(), unmanagedValues, unmanagedGraph);
417  new (&hostJacTpetraBlocks(row,col)) LocalMatrixType(unmanagedMatrix);
418  }
419 
420  hostBlockExistsInJac(row,col) = 1;
421  }
422  else {
423  hostBlockExistsInJac(row,col) = 0;
424  }
425  }
426  }
427  typename PHX::View<LocalMatrixType**>
428  jacTpetraBlocks("panzer::ScatterResidual_BlockedTpetra<Jacobian>::jacTpetraBlocks",numFieldBlocks,numFieldBlocks);
429  Kokkos::deep_copy(jacTpetraBlocks,hostJacTpetraBlocks);
430  Kokkos::deep_copy(blockExistsInJac,hostBlockExistsInJac);
431 
432  // worksetLIDs is larger for Jacobian than Residual fill. Need the
433  // entire set of field offsets for derivative indexing no matter
434  // which block row you are scattering. The residual only needs the
435  // lids for the sub-block that it is scattering to. The subviews
436  // below are to offset the LID blocks correctly.
437  const auto& globalIndexers = globalIndexer_->getFieldDOFManagers();
438  auto blockOffsets_h = Kokkos::create_mirror_view(blockOffsets_);
439  Kokkos::deep_copy(blockOffsets_h, blockOffsets_);
440  for (size_t block=0; block < globalIndexers.size(); ++block) {
441  const auto subviewOfBlockLIDs = Kokkos::subview(worksetLIDs_,Kokkos::ALL(), std::make_pair(blockOffsets_h(block),blockOffsets_h(block+1)));
442  globalIndexers[block]->getElementLIDs(localCellIds,subviewOfBlockLIDs);
443  }
444 
445  // Loop over scattered fields
446  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
447 
448  const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
449  typename Tpetra::Vector<double,LO,GO,PHX::Device>::dual_view_type::t_dev kokkosResidual;
450  if (haveResidual) {
451  auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(blockRowIndex),true))->getTpetraVector());
452  kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::OverwriteAll);
453  }
454 
455  // Class data fields for lambda capture
456  const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
457  const PHX::View<const LO**> worksetLIDs = worksetLIDs_;
458  const PHX::View<const ScalarT**> fieldValues = scatterFields_[fieldIndex].get_static_view();
459  const PHX::View<const LO*> blockOffsets = blockOffsets_;
460 
461  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
462  LO cLIDs[256];
463  typename Sacado::ScalarType<ScalarT>::type vals[256];
464 
465  for(int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
466  typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
467  typename FieldType::array_type::reference_type tmpFieldVal = fieldValues(cell,basis);
468  const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
469 
470  if (haveResidual)
471  Kokkos::atomic_add(&kokkosResidual(rowLID,0), tmpFieldVal.val());
472 
473  for (int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
474  if (blockExistsInJac(blockRowIndex,blockColIndex)) {
475  const int start = blockOffsets(blockColIndex);
476  const int stop = blockOffsets(blockColIndex+1);
477  const int sensSize = stop-start;
478  // Views may be padded. Use contiguous arrays here
479  for (int i=0; i < sensSize; ++i) {
480  cLIDs[i] = worksetLIDs(cell,start+i);
481  vals[i] = tmpFieldVal.fastAccessDx(start+i);
482  }
483  jacTpetraBlocks(blockRowIndex,blockColIndex).sumIntoValues(rowLID,cLIDs,sensSize,vals,true,true);
484  }
485  }
486  }
487  });
488 
489  }
490 
491  // Placement delete on view of matrices
492  for (int row=0; row < numFieldBlocks; ++row) {
493  for (int col=0; col < numFieldBlocks; ++col) {
494  if (hostBlockExistsInJac(row,col)) {
495  hostJacTpetraBlocks(row,col).~CrsMatrix();
496  }
497  }
498  }
499 
500 }
501 
502 // **********************************************************************
503 
504 #endif
Pushes residual values into the residual vector for a Newton-based solve.
void postRegistrationSetup(typename TRAITS::SetupData, PHX::FieldManager< TRAITS > &)
ScatterResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager > &)
FieldType
The type of discretization to use for a field pattern.
PHX::View< const int * > offsets
std::string block_id
DEPRECATED - use: getElementBlock()
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we&#39;ll contribute, or in which we&#39;ll store, the result of computing this integral...