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