Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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>
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 
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>
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 
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] = Kokkos::View<int*,PHX::Device>("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  PHX::Device::fence();
173  }
174 
175  // We will use one workset lid view for all fields, but has to be
176  // sized big enough to hold the largest elementBlockGIDCount in the
177  // ProductVector.
178  worksetLIDs_ = Kokkos::View<LO**,PHX::Device>("ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
179  scatterFields_[0].extent(0),
180  maxElementBlockGIDCount);
181 }
182 
183 // **********************************************************************
184 template <typename TRAITS,typename LO,typename GO,typename NodeT>
186 preEvaluate(typename TRAITS::PreEvalData d)
187 {
188  // extract linear object container
189  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
190 }
191 
192 // **********************************************************************
193 template <typename TRAITS,typename LO,typename GO,typename NodeT>
195 evaluateFields(typename TRAITS::EvalData workset)
196 {
197  using Teuchos::RCP;
198  using Teuchos::rcp_dynamic_cast;
199  using Thyra::VectorBase;
201 
202  const auto& localCellIds = this->wda(workset).cell_local_ids_k;
203  const RCP<ProductVectorBase<double>> thyraBlockResidual = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f(),true);
204 
205  // Loop over scattered fields
206  int currentWorksetLIDSubBlock = -1;
207  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
208  // workset LIDs only change for different sub blocks
209  if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
210  fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
211  currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
212  }
213 
214  const auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),true))->getTpetraVector());
215  const auto& kokkosResidual = tpetraResidual.template getLocalView<PHX::mem_space>();
216 
217  // Class data fields for lambda capture
218  const auto& fieldOffsets = fieldOffsets_[fieldIndex];
219  const auto& worksetLIDs = worksetLIDs_;
220  const auto& fieldValues = scatterFields_[fieldIndex];
221 
222  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
223  for(int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
224  const int lid = worksetLIDs(cell,fieldOffsets(basis));
225  Kokkos::atomic_add(&kokkosResidual(lid,0), fieldValues(cell,basis));
226  }
227  });
228  }
229 }
230 
231 // **********************************************************************
232 // Specialization: Jacobian
233 // **********************************************************************
234 
235 template <typename TRAITS,typename LO,typename GO,typename NodeT>
238  const Teuchos::ParameterList& p)
239  : globalIndexer_(indexer)
240  , globalDataKey_("Residual Scatter Container")
241 {
242  std::string scatterName = p.get<std::string>("Scatter Name");
243  scatterHolder_ =
244  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
245 
246  // get names to be evaluated
247  const std::vector<std::string>& names =
248  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
249 
250  // grab map from evaluated names to field names
251  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
252 
254  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
255 
256  // build the vector of fields that this is dependent on
257  scatterFields_.resize(names.size());
258  for (std::size_t eq = 0; eq < names.size(); ++eq) {
259  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
260 
261  // tell the field manager that we depend on this field
262  this->addDependentField(scatterFields_[eq]);
263  }
264 
265  // this is what this evaluator provides
266  this->addEvaluatedField(*scatterHolder_);
267 
268  if (p.isType<std::string>("Global Data Key"))
269  globalDataKey_ = p.get<std::string>("Global Data Key");
270 
271  this->setName(scatterName+" Scatter Residual (Jacobian)");
272 }
273 
274 // **********************************************************************
275 template <typename TRAITS,typename LO,typename GO,typename NodeT>
277 postRegistrationSetup(typename TRAITS::SetupData d,
278  PHX::FieldManager<TRAITS>& /* fm */)
279 {
280  const Workset & workset_0 = (*d.worksets_)[0];
281  const std::string blockId = this->wda(workset_0).block_id;
282 
283  fieldIds_.resize(scatterFields_.size());
284  fieldOffsets_.resize(scatterFields_.size());
285  productVectorBlockIndex_.resize(scatterFields_.size());
286  for (std::size_t fd=0; fd < scatterFields_.size(); ++fd) {
287  const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
288  const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
289  productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
290  const auto& fieldGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
291  fieldIds_[fd] = fieldGlobalIndexer->getFieldNum(fieldName); // Field number in the sub-global-indexer
292 
293  const std::vector<int>& offsets = globalIndexer_->getGIDFieldOffsets(blockId,globalFieldNum);
294  fieldOffsets_[fd] = Kokkos::View<int*,PHX::Device>("ScatterResidual_BlockedTpetra(Jacobian):fieldOffsets",offsets.size());
295  auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
296  for (std::size_t i=0; i < offsets.size(); ++i)
297  hostOffsets(i) = offsets[i];
298  Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
299  PHX::Device::fence();
300  }
301 
302  // This is sized differently than the Residual implementation since
303  // we need the LIDs for all sub-blocks, not just the single
304  // sub-block for the field residual scatter.
305  int elementBlockGIDCount = 0;
306  for (const auto blockDOFMgr : globalIndexer_->getFieldDOFManagers())
307  elementBlockGIDCount += blockDOFMgr->getElementBlockGIDCount(blockId);
308 
309  worksetLIDs_ = Kokkos::View<LO**,PHX::Device>("ScatterResidual_BlockedTpetra(Jacobian):worksetLIDs",
310  scatterFields_[0].extent(0),
311  elementBlockGIDCount);
312 
313  // Compute the block offsets
314  const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
315  const int numBlocks = static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
316  blockOffsets_ = Kokkos::View<LO*,PHX::Device>("ScatterResidual_BlockedTpetra(Jacobian):blockOffsets_",
317  numBlocks+1); // Number of fields, plus a sentinel
318  const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
319  for (int blk=0;blk<numBlocks;blk++) {
320  int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
321  hostBlockOffsets(blk) = blockOffset;
322  }
323  blockOffsets_(numBlocks) = blockOffsets_(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
324  Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
325 
326  // Make sure the that hard coded derivative dimension in the
327  // evaluate call is large enough to hold all derivatives for each
328  // sub block load
329  for (int blk=0;blk<numBlocks;blk++) {
330  const int blockDerivativeSize = hostBlockOffsets(blk+1) - hostBlockOffsets(blk);
331  TEUCHOS_TEST_FOR_EXCEPTION(blockDerivativeSize > 256, std::runtime_error,
332  "ERROR: the derivative dimension for sub block "
333  << blk << "with a value of " << blockDerivativeSize
334  << "is larger than the size allocated for cLIDs and vals "
335  << "in the evaluate call! You must manually increase the "
336  << "size and recompile!");
337  }
338 
339  PHX::Device::fence();
340 }
341 
342 // **********************************************************************
343 template <typename TRAITS,typename LO,typename GO,typename NodeT>
345 preEvaluate(typename TRAITS::PreEvalData d)
346 {
347  using Teuchos::RCP;
348  using Teuchos::rcp_dynamic_cast;
349 
350  // extract linear object container
351  blockedContainer_ = rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_));
352 
353  if(blockedContainer_==Teuchos::null) {
354  RCP<const LOCPair_GlobalEvaluationData> gdata = rcp_dynamic_cast<const LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true);
355  blockedContainer_ = rcp_dynamic_cast<const ContainerType>(gdata->getGhostedLOC());
356  }
357 }
358 
359 // **********************************************************************
360 template <typename TRAITS,typename LO,typename GO,typename NodeT>
362 evaluateFields(typename TRAITS::EvalData workset)
363 {
364  using Teuchos::RCP;
365  using Teuchos::rcp_dynamic_cast;
366  using Thyra::VectorBase;
369 
370  const auto& localCellIds = this->wda(workset).cell_local_ids_k;
371 
372  const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
373  const RCP<const ContainerType> blockedContainer = blockedContainer_;
374  const RCP<ProductVectorBase<double>> thyraBlockResidual = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f());
375  const bool haveResidual = Teuchos::nonnull(thyraBlockResidual);
376  const RCP<BlockedLinearOpBase<double>> Jac = rcp_dynamic_cast<BlockedLinearOpBase<double> >(blockedContainer_->get_A(),true);
377 
378  // Get the local data for the sub-block crs matrices. First allocate
379  // on host and then deep_copy to device. The sub-blocks are
380  // unmanaged since they are allocated and ref counted separately on
381  // host.
382  using LocalMatrixType = KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
383  typename Kokkos::View<LocalMatrixType**,PHX::Device>::HostMirror
384  hostJacTpetraBlocks("panzer::ScatterResidual_BlockTpetra<Jacobian>::hostJacTpetraBlocks", numFieldBlocks,numFieldBlocks);
385 
386  Kokkos::View<int**,PHX::Device> blockExistsInJac = Kokkos::View<int**,PHX::Device>("blockExistsInJac_",numFieldBlocks,numFieldBlocks);
387  auto hostBlockExistsInJac = Kokkos::create_mirror_view(blockExistsInJac);
388 
389  for (int row=0; row < numFieldBlocks; ++row) {
390  for (int col=0; col < numFieldBlocks; ++col) {
391  const auto thyraTpetraOperator = rcp_dynamic_cast<Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(Jac->getNonconstBlock(row,col),false);
392  if (nonnull(thyraTpetraOperator)) {
393 
394  // HACK to enforce views in the CrsGraph to be
395  // Unmanaged. Passing in the MemoryTrait<Unmanaged> doesn't
396  // work as the CrsGraph in the CrsMatrix ignores the
397  // MemoryTrait. Need to use the runtime constructor by passing
398  // in points to ensure Unmanaged. See:
399  // https://github.com/kokkos/kokkos/issues/1581
400 
401  // These two lines are the original code we can revert to when #1581 is fixed.
402  // const auto crsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),true);
403  // new (&hostJacTpetraBlocks(row,col)) KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>> (crsMatrix->getLocalMatrix());
404 
405  // Instead do this
406  {
407  // Grab the local managed matrix and graph
408  const auto tpetraCrsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),true);
409  const auto managedMatrix = tpetraCrsMatrix->getLocalMatrix();
410  const auto managedGraph = managedMatrix.graph;
411 
412  // Create runtime unmanaged versions
413  using StaticCrsGraphType = typename LocalMatrixType::StaticCrsGraphType;
414  StaticCrsGraphType unmanagedGraph;
415  unmanagedGraph.entries = typename StaticCrsGraphType::entries_type(managedGraph.entries.data(),managedGraph.entries.extent(0));
416  unmanagedGraph.row_map = typename StaticCrsGraphType::row_map_type(managedGraph.row_map.data(),managedGraph.row_map.extent(0));
417  unmanagedGraph.row_block_offsets = typename StaticCrsGraphType::row_block_type(managedGraph.row_block_offsets.data(),managedGraph.row_block_offsets.extent(0));
418 
419  typename LocalMatrixType::values_type unmanagedValues(managedMatrix.values.data(),managedMatrix.values.extent(0));
420  LocalMatrixType unmanagedMatrix(managedMatrix.values.label(), managedMatrix.numCols(), unmanagedValues, unmanagedGraph);
421  new (&hostJacTpetraBlocks(row,col)) LocalMatrixType(unmanagedMatrix);
422  }
423 
424  hostBlockExistsInJac(row,col) = 1;
425  }
426  else {
427  hostBlockExistsInJac(row,col) = 0;
428  }
429  }
430  }
431  typename Kokkos::View<LocalMatrixType**,PHX::Device>
432  jacTpetraBlocks("panzer::ScatterResidual_BlockedTpetra<Jacobian>::jacTpetraBlocks",numFieldBlocks,numFieldBlocks);
433  Kokkos::deep_copy(jacTpetraBlocks,hostJacTpetraBlocks);
434  Kokkos::deep_copy(blockExistsInJac,hostBlockExistsInJac);
435  PHX::Device::fence();
436 
437  // worksetLIDs is larger for Jacobian than Residual fill. Need the
438  // entire set of field offsets for derivative indexing no matter
439  // which block row you are scattering. The residual only needs the
440  // lids for the sub-block that it is scattering to. The subviews
441  // below are to offset the LID blocks correctly.
442  const auto& globalIndexers = globalIndexer_->getFieldDOFManagers();
443  for (size_t block=0; block < globalIndexers.size(); ++block) {
444  const auto subviewOfBlockLIDs = Kokkos::subview(worksetLIDs_,Kokkos::ALL(), std::make_pair(blockOffsets_(block),blockOffsets_(block+1)));
445  globalIndexers[block]->getElementLIDs(localCellIds,subviewOfBlockLIDs);
446  }
447 
448  // Loop over scattered fields
449  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
450 
451  const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
452  typename Tpetra::Vector<double,LO,GO,PHX::Device>::dual_view_type::t_dev kokkosResidual;
453  if (haveResidual) {
454  const auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(blockRowIndex),true))->getTpetraVector());
455  kokkosResidual = tpetraResidual.template getLocalView<PHX::mem_space>();
456  }
457 
458  // Class data fields for lambda capture
459  const Kokkos::View<const int*,PHX::Device> fieldOffsets = fieldOffsets_[fieldIndex];
460  const Kokkos::View<const LO**,PHX::Device> worksetLIDs = worksetLIDs_;
461  const PHX::View<const ScalarT**> fieldValues = scatterFields_[fieldIndex].get_static_view();
462  const Kokkos::View<const LO*,PHX::Device> blockOffsets = blockOffsets_;
463 
464  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
465  LO cLIDs[256];
466  typename Sacado::ScalarType<ScalarT>::type vals[256];
467 
468  for(int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
469  typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
470  typename FieldType::array_type::reference_type tmpFieldVal = fieldValues(cell,basis);
471  const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
472 
473  if (haveResidual)
474  Kokkos::atomic_add(&kokkosResidual(rowLID,0), tmpFieldVal.val());
475 
476  for (int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
477  if (blockExistsInJac(blockRowIndex,blockColIndex)) {
478  const int start = blockOffsets(blockColIndex);
479  const int stop = blockOffsets(blockColIndex+1);
480  const int sensSize = stop-start;
481  // Views may be padded. Use contiguous arrays here
482  for (int i=0; i < sensSize; ++i) {
483  cLIDs[i] = worksetLIDs(cell,start+i);
484  vals[i] = tmpFieldVal.fastAccessDx(start+i);
485  }
486  jacTpetraBlocks(blockRowIndex,blockColIndex).sumIntoValues(rowLID,cLIDs,sensSize,vals,true,true);
487  }
488  }
489  }
490  });
491 
492  }
493 
494  // Placement delete on view of matrices
495  for (int row=0; row < numFieldBlocks; ++row) {
496  for (int col=0; col < numFieldBlocks; ++col) {
497  if (hostBlockExistsInJac(row,col)) {
498  hostJacTpetraBlocks(row,col).~CrsMatrix();
499  }
500  }
501  }
502 
503 }
504 
505 // **********************************************************************
506 
507 #endif
Pushes residual values into the residual vector for a Newton-based solve.
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
void postRegistrationSetup(typename TRAITS::SetupData, PHX::FieldManager< TRAITS > &)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
FieldType
The type of discretization to use for a field pattern.
Teuchos::RCP< LinearObjContainer > getGhostedLOC() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool nonnull(const boost::shared_ptr< T > &p)
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...
bool isType(const std::string &name) const
ScatterResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &)
Kokkos::View< const int *, PHX::Device > offsets