Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_GatherSolution_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_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
44#define PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
45
46#include "Teuchos_Assert.hpp"
47#include "Phalanx_DataLayout.hpp"
48
51#include "Panzer_PureBasis.hpp"
52#include "Panzer_TpetraLinearObjFactory.hpp"
56
57#include "Teuchos_FancyOStream.hpp"
58
59#include "Thyra_SpmdVectorBase.hpp"
60#include "Thyra_ProductVectorBase.hpp"
61
62#include "Tpetra_Map.hpp"
63
64template <typename EvalT,typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
67 const Teuchos::RCP<const BlockedDOFManager> & indexer,
68 const Teuchos::ParameterList& p)
69{
70 const std::vector<std::string>& names =
71 *(p.get< Teuchos::RCP< std::vector<std::string> > >("DOF Names"));
72
73 Teuchos::RCP<panzer::PureBasis> basis =
74 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis");
75
76 for (std::size_t fd = 0; fd < names.size(); ++fd) {
77 PHX::MDField<ScalarT,Cell,NODE> field = PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
78 this->addEvaluatedField(field.fieldTag());
79 }
80
81 this->setName("Gather Solution");
82}
83
84// **********************************************************************
85// Specialization: Residual
86// **********************************************************************
87
88template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
91 const Teuchos::RCP<const BlockedDOFManager> & indexer,
92 const Teuchos::ParameterList& p)
93 : globalIndexer_(indexer)
94 , has_tangent_fields_(false)
95{
96 typedef std::vector< std::vector<std::string> > vvstring;
97
99 input.setParameterList(p);
100
101 const std::vector<std::string> & names = input.getDofNames();
102 Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
103 const vvstring & tangent_field_names = input.getTangentNames();
104
105 indexerNames_ = input.getIndexerNames();
106 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
107 globalDataKey_ = input.getGlobalDataKey();
108
109 // allocate fields
110 gatherFields_.resize(names.size());
111 for (std::size_t fd = 0; fd < names.size(); ++fd) {
112 gatherFields_[fd] =
113 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
114 this->addEvaluatedField(gatherFields_[fd]);
115 }
116
117 // Setup dependent tangent fields if requested
118 if (tangent_field_names.size()>0) {
119 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
120
121 has_tangent_fields_ = true;
122 tangentFields_.resize(gatherFields_.size());
123 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
124 tangentFields_[fd].resize(tangent_field_names[fd].size());
125 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
126 tangentFields_[fd][i] =
127 PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
128 this->addDependentField(tangentFields_[fd][i]);
129 }
130 }
131 }
132
133 // figure out what the first active name is
134 std::string firstName = "<none>";
135 if(names.size()>0)
136 firstName = names[0];
137
138 std::string n = "GatherSolution (BlockedTpetra): "+firstName+" (Residual)";
139 this->setName(n);
140}
141
142// **********************************************************************
143template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
145postRegistrationSetup(typename TRAITS::SetupData d,
147{
148 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
149
150 const Workset & workset_0 = (*d.worksets_)[0];
151 const std::string blockId = this->wda(workset_0).block_id;
152
153 fieldIds_.resize(gatherFields_.size());
154 fieldOffsets_.resize(gatherFields_.size());
155 fieldGlobalIndexers_.resize(gatherFields_.size());
156 productVectorBlockIndex_.resize(gatherFields_.size());
157 int maxElementBlockGIDCount = -1;
158 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
159 // get field ID from DOF manager
160 const std::string& fieldName = indexerNames_[fd];
161 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
162 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
163 fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
164 fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName); // Field number in the sub-global-indexer
165
166 const std::vector<int>& offsets = fieldGlobalIndexers_[fd]->getGIDFieldOffsets(blockId,fieldIds_[fd]);
167 fieldOffsets_[fd] = PHX::View<int*>("GatherSolution_BlockedTpetra(Residual):fieldOffsets",offsets.size());
168 auto hostFieldOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
169 for(std::size_t i=0; i < offsets.size(); ++i)
170 hostFieldOffsets(i) = offsets[i];
171 Kokkos::deep_copy(fieldOffsets_[fd],hostFieldOffsets);
172
173 maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
174 }
175
176 // We will use one workset lid view for all fields, but has to be
177 // sized big enough to hold the largest elementBlockGIDCount in the
178 // ProductVector.
179 worksetLIDs_ = PHX::View<LO**>("GatherSolution_BlockedTpetra(Residual):worksetLIDs",
180 gatherFields_[0].extent(0),
181 maxElementBlockGIDCount);
182
183 indexerNames_.clear(); // Don't need this anymore
184}
185
186// **********************************************************************
187template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
189preEvaluate(typename TRAITS::PreEvalData d)
190{
191 // extract linear object container
192 blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
193}
194
195// **********************************************************************
196template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
198evaluateFields(typename TRAITS::EvalData workset)
199{
200 using Teuchos::RCP;
201 using Teuchos::rcp_dynamic_cast;
202 using Thyra::VectorBase;
204
205 const PHX::View<const int*>& localCellIds = this->wda(workset).cell_local_ids_k;
206
207 RCP<ProductVectorBase<ScalarT>> thyraBlockSolution;
208 if (useTimeDerivativeSolutionVector_)
209 thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_dxdt(),true);
210 else
211 thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_x(),true);
212
213 // Loop over gathered fields
214 int currentWorksetLIDSubBlock = -1;
215 for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
216 // workset LIDs only change for different sub blocks
217 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
218 const std::string blockId = this->wda(workset).block_id;
219 const int num_dofs = fieldGlobalIndexers_[fieldIndex]->getElementBlockGIDCount(blockId);
220 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_,num_dofs);
221 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
222 }
223
224 const auto& tpetraSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<ScalarT,LO,GO,NodeT>>(thyraBlockSolution->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),true))->getTpetraVector());
225 const auto& kokkosSolution = tpetraSolution.getLocalViewDevice(Tpetra::Access::ReadOnly);
226
227 // Class data fields for lambda capture
228 const auto& fieldOffsets = fieldOffsets_[fieldIndex];
229 const auto& worksetLIDs = worksetLIDs_;
230 const auto& fieldValues = gatherFields_[fieldIndex];
231
232 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
233 for(int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
234 const int lid = worksetLIDs(cell,fieldOffsets(basis));
235 fieldValues(cell,basis) = kokkosSolution(lid,0);
236 }
237 });
238 }
239
240}
241
242// **********************************************************************
243// Specialization: Tangent
244// **********************************************************************
245
246template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
249 const Teuchos::RCP<const BlockedDOFManager> & indexer,
250 const Teuchos::ParameterList& p)
251 : gidIndexer_(indexer)
252 , has_tangent_fields_(false)
253{
254 typedef std::vector< std::vector<std::string> > vvstring;
255
257 input.setParameterList(p);
258
259 const std::vector<std::string> & names = input.getDofNames();
260 Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
261 const vvstring & tangent_field_names = input.getTangentNames();
262
263 indexerNames_ = input.getIndexerNames();
264 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
265 globalDataKey_ = input.getGlobalDataKey();
266
267 // allocate fields
268 gatherFields_.resize(names.size());
269 for (std::size_t fd = 0; fd < names.size(); ++fd) {
270 gatherFields_[fd] =
271 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
272 this->addEvaluatedField(gatherFields_[fd]);
273 }
274
275 // Setup dependent tangent fields if requested
276 if (tangent_field_names.size()>0) {
277 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
278
279 has_tangent_fields_ = true;
280 tangentFields_.resize(gatherFields_.size());
281 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
282 tangentFields_[fd].resize(tangent_field_names[fd].size());
283 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
284 tangentFields_[fd][i] =
285 PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
286 this->addDependentField(tangentFields_[fd][i]);
287 }
288 }
289 }
290
291 // figure out what the first active name is
292 std::string firstName = "<none>";
293 if(names.size()>0)
294 firstName = names[0];
295
296 std::string n = "GatherSolution (BlockedTpetra): "+firstName+" (Tangent)";
297 this->setName(n);
298}
299
300// **********************************************************************
301template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
303postRegistrationSetup(typename TRAITS::SetupData /* d */,
305{
306 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
307
308 fieldIds_.resize(gatherFields_.size());
309
310 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
311 // get field ID from DOF manager
312 const std::string& fieldName = indexerNames_[fd];
313 fieldIds_[fd] = gidIndexer_->getFieldNum(fieldName);
314 }
315
316 indexerNames_.clear(); // Don't need this anymore
317}
318
319// **********************************************************************
320template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
322preEvaluate(typename TRAITS::PreEvalData d)
323{
324 // extract linear object container
325 blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
326}
327
328// **********************************************************************
329template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
331evaluateFields(typename TRAITS::EvalData workset)
332{
333 using Teuchos::RCP;
334 using Teuchos::ArrayRCP;
335 using Teuchos::ptrFromRef;
336 using Teuchos::rcp_dynamic_cast;
337
338 using Thyra::VectorBase;
339 using Thyra::SpmdVectorBase;
341
342 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
343 out.setShowProcRank(true);
344 out.setOutputToRootOnly(-1);
345
346 std::vector<std::pair<int,GO> > GIDs;
347 std::vector<LO> LIDs;
348
349 // for convenience pull out some objects from workset
350 std::string blockId = this->wda(workset).block_id;
351 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
352
353 Teuchos::RCP<ProductVectorBase<double> > x;
354 if (useTimeDerivativeSolutionVector_)
355 x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
356 else
357 x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
358
359 // gather operation for each cell in workset
360 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
361 LO cellLocalId = localCellIds[worksetCellIndex];
362
363 gidIndexer_->getElementGIDsPair(cellLocalId,GIDs,blockId);
364
365 // caculate the local IDs for this element
366 LIDs.resize(GIDs.size());
367 for(std::size_t i=0;i<GIDs.size();i++) {
368 // used for doing local ID lookups
369 RCP<const MapType> x_map = blockedContainer_->getMapForBlock(GIDs[i].first);
370
371 LIDs[i] = x_map->getLocalElement(GIDs[i].second);
372 }
373
374 // loop over the fields to be gathered
375 Teuchos::ArrayRCP<const double> local_x;
376 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
377 int fieldNum = fieldIds_[fieldIndex];
378 int indexerId = gidIndexer_->getFieldBlock(fieldNum);
379
380 // grab local data for inputing
381 RCP<SpmdVectorBase<double> > block_x = rcp_dynamic_cast<SpmdVectorBase<double> >(x->getNonconstVectorBlock(indexerId));
382 block_x->getLocalData(ptrFromRef(local_x));
383
384 const std::vector<int> & elmtOffset = gidIndexer_->getGIDFieldOffsets(blockId,fieldNum);
385
386 // loop over basis functions and fill the fields
387 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
388 int offset = elmtOffset[basis];
389 int lid = LIDs[offset];
390
391 if (!has_tangent_fields_)
392 (gatherFields_[fieldIndex])(worksetCellIndex,basis) = local_x[lid];
393 else {
394 (gatherFields_[fieldIndex])(worksetCellIndex,basis).val() = local_x[lid];
395 for (std::size_t i=0; i<tangentFields_[fieldIndex].size(); ++i)
396 (gatherFields_[fieldIndex])(worksetCellIndex,basis).fastAccessDx(i) =
397 tangentFields_[fieldIndex][i](worksetCellIndex,basis).val();
398 }
399 }
400 }
401 }
402}
403
404// **********************************************************************
405// Specialization: Jacobian
406// **********************************************************************
407
408template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
411 const Teuchos::RCP<const BlockedDOFManager> & indexer,
412 const Teuchos::ParameterList& p)
413 : globalIndexer_(indexer)
414{
416 input.setParameterList(p);
417
418 const std::vector<std::string> & names = input.getDofNames();
419 Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
420
421 indexerNames_ = input.getIndexerNames();
422 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
423 globalDataKey_ = input.getGlobalDataKey();
424
425 disableSensitivities_ = !input.firstSensitivitiesAvailable();
426
427 gatherFields_.resize(names.size());
428 for (std::size_t fd = 0; fd < names.size(); ++fd) {
429 PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->functional);
430 gatherFields_[fd] = f;
431 this->addEvaluatedField(gatherFields_[fd]);
432 }
433
434 // figure out what the first active name is
435 std::string firstName = "<none>";
436 if(names.size()>0)
437 firstName = names[0];
438
439 // print out convenience
440 if(disableSensitivities_) {
441 std::string n = "GatherSolution (BlockedTpetra, No Sensitivities): "+firstName+" (Jacobian)";
442 this->setName(n);
443 }
444 else {
445 std::string n = "GatherSolution (BlockedTpetra): "+firstName+" ("+PHX::print<EvalT>()+") ";
446 this->setName(n);
447 }
448}
449
450// **********************************************************************
451template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
453postRegistrationSetup(typename TRAITS::SetupData d,
455{
456 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
457
458 const Workset & workset_0 = (*d.worksets_)[0];
459 const std::string blockId = this->wda(workset_0).block_id;
460
461 fieldIds_.resize(gatherFields_.size());
462 fieldOffsets_.resize(gatherFields_.size());
463 productVectorBlockIndex_.resize(gatherFields_.size());
464 int maxElementBlockGIDCount = -1;
465 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
466
467 const std::string fieldName = indexerNames_[fd];
468 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
469 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
470 const auto& subGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
471 fieldIds_[fd] = subGlobalIndexer->getFieldNum(fieldName); // Field number in the sub-global-indexer
472
473 const std::vector<int>& offsets = subGlobalIndexer->getGIDFieldOffsets(blockId,fieldIds_[fd]);
474 fieldOffsets_[fd] = PHX::View<int*>("GatherSolution_BlockedTpetra(Jacobian):fieldOffsets",offsets.size());
475 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
476 for (std::size_t i=0; i < offsets.size(); ++i)
477 hostOffsets(i) = offsets[i];
478 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
479 maxElementBlockGIDCount = std::max(subGlobalIndexer->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
480 }
481
482 // We will use one workset lid view for all fields, but has to be
483 // sized big enough to hold the largest elementBlockGIDCount in the
484 // ProductVector.
485 worksetLIDs_ = PHX::View<LO**>("ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
486 gatherFields_[0].extent(0),
487 maxElementBlockGIDCount);
488
489 // Compute the block offsets
490 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
491 const int numBlocks = static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
492 blockOffsets_ = PHX::View<LO*>("GatherSolution_BlockedTpetra(Jacobian):blockOffsets_",
493 numBlocks+1); // Number of blocks, plus a sentinel
494 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
495 for (int blk=0;blk<numBlocks;++blk) {
496 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
497 hostBlockOffsets(blk) = blockOffset;
498 }
499 hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
500 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
501
502 indexerNames_.clear(); // Don't need this anymore
503}
504
505template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
507preEvaluate(typename TRAITS::PreEvalData d)
508{
509 // extract linear object container
510 blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
511}
512
513// **********************************************************************
514template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
516evaluateFields(typename TRAITS::EvalData workset)
517{
518 using Teuchos::RCP;
519 using Teuchos::rcp_dynamic_cast;
520 using Thyra::VectorBase;
522
523 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
524
525 RealType seedValue = RealType(0.0);
526 RCP<ProductVectorBase<double>> blockedSolution;
527 if (useTimeDerivativeSolutionVector_) {
528 blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
529 seedValue = workset.alpha;
530 }
531 else {
532 blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
533 seedValue = workset.beta;
534 }
535
536 // turn off sensitivies: this may be faster if we don't expand the term
537 // but I suspect not because anywhere it is used the full complement of
538 // sensitivies will be needed anyway.
539 if(disableSensitivities_)
540 seedValue = 0.0;
541
542 // Loop over fields to gather
543 int currentWorksetLIDSubBlock = -1;
544 for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
545 // workset LIDs only change if in different sub blocks
546 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
547 const auto& blockIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]];
548 const std::string blockId = this->wda(workset).block_id;
549 const int num_dofs = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]]->getElementBlockGIDCount(blockId);
550 blockIndexer->getElementLIDs(localCellIds,worksetLIDs_,num_dofs);
551 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
552 }
553
554 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
555 const auto& subblockSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(blockedSolution->getNonconstVectorBlock(blockRowIndex),true))->getTpetraVector());
556 const auto kokkosSolution = subblockSolution.getLocalViewDevice(Tpetra::Access::ReadOnly);
557
558 // Class data fields for lambda capture
559 const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
560 const PHX::View<const LO**> worksetLIDs = worksetLIDs_;
561 const PHX::View<ScalarT**> fieldValues = gatherFields_[fieldIndex].get_static_view();
562 const PHX::View<const LO*> blockOffsets = blockOffsets_;
563 auto blockOffsets_h = Kokkos::create_mirror_view(blockOffsets);
564 Kokkos::deep_copy(blockOffsets_h, blockOffsets);
565 const int blockStart = blockOffsets_h(blockRowIndex);
566
567 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
568 for (int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
569 const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
570 fieldValues(cell,basis).zero();
571 fieldValues(cell,basis).val() = kokkosSolution(rowLID,0);
572 fieldValues(cell,basis).fastAccessDx(blockStart+fieldOffsets(basis)) = seedValue;
573 }
574 });
575
576 }
577}
578
579// **********************************************************************
580
581#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.
Gathers solution values from the Newton solution vector into the nodal fields of the field manager.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
GatherSolution_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager > &indexer)
std::string getGlobalDataKey() const
Name of the global evaluation data container to use for the source vector (all types)
void setParameterList(const Teuchos::ParameterList &pl)
const std::vector< std::vector< std::string > > & getTangentNames() const
Get the name of the tangent fields (tangent only)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
const std::vector< std::string > & getIndexerNames() const
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
std::string block_id
DEPRECATED - use: getElementBlock()