Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_GatherSolution_BlockedEpetra_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_GatherSolution_BlockedEpetra_impl_hpp__
44#define __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
45
47//
48// Include Files
49//
51
52// Panzer
56#include "Panzer_PureBasis.hpp"
60
61// Phalanx
62#include "Phalanx_DataLayout.hpp"
63
64// Teuchos
65#include "Teuchos_Assert.hpp"
66#include "Teuchos_FancyOStream.hpp"
67
68// Thyra
69#include "Thyra_ProductVectorBase.hpp"
70#include "Thyra_SpmdVectorBase.hpp"
71
73//
74// Initializing Constructor (Residual Specialization)
75//
77template<typename TRAITS, typename LO, typename GO>
81 const std::vector<Teuchos::RCP<const GlobalIndexer>>&
82 indexers,
83 const Teuchos::ParameterList& p)
84 :
85 indexers_(indexers),
86 hasTangentFields_(false)
87{
89 using PHX::MDField;
90 using PHX::print;
91 using std::size_t;
92 using std::string;
93 using std::vector;
94 using Teuchos::RCP;
95 using vvstring = std::vector<std::vector<std::string>>;
97 input.setParameterList(p);
98 const vector<string>& names = input.getDofNames();
99 RCP<const PureBasis> basis = input.getBasis();
100 const vvstring& tangentFieldNames = input.getTangentNames();
101 indexerNames_ = input.getIndexerNames();
102 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
103 globalDataKey_ = input.getGlobalDataKey();
104
105 // Allocate the fields.
106 int numFields(names.size());
107 gatherFields_.resize(numFields);
108 for (int fd(0); fd < numFields; ++fd)
109 {
110 gatherFields_[fd] =
111 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
112 this->addEvaluatedField(gatherFields_[fd]);
113 } // end loop over names
114
115 // Setup the dependent tangent fields, if requested.
116 if (tangentFieldNames.size() > 0)
117 {
118 TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size());
119 hasTangentFields_ = true;
120 tangentFields_.resize(numFields);
121 for (int fd(0); fd < numFields; ++fd)
122 {
123 int numTangentFields(tangentFieldNames[fd].size());
124 tangentFields_[fd].resize(numTangentFields);
125 for (int i(0); i < numTangentFields; ++i)
126 {
127 tangentFields_[fd][i] =
128 MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
129 basis->functional);
130 this->addDependentField(tangentFields_[fd][i]);
131 } // end loop over tangentFieldNames[fd]
132 } // end loop over gatherFields_
133 } // end if we have tangent fields
134
135 // Figure out what the first active name is.
136 string firstName("<none>");
137 if (numFields > 0)
138 firstName = names[0];
139 string n("GatherSolution (BlockedEpetra): " + firstName + " (" +
140 print<EvalT>() + ")");
141 this->setName(n);
142} // end of Initializing Constructor (Residual Specialization)
143
145//
146// postRegistrationSetup() (Residual Specialization)
147//
149template<typename TRAITS, typename LO, typename GO>
150void
154 typename TRAITS::SetupData /* d */,
156{
157 using std::size_t;
158 using std::string;
159 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
160 int numFields(gatherFields_.size());
161 indexerIds_.resize(numFields);
162 subFieldIds_.resize(numFields);
163 for (int fd(0); fd < numFields; ++fd)
164 {
165 // Get the field ID from the DOF manager.
166 const string& fieldName(indexerNames_[fd]);
167 indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
168 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
169 TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
170 } // end loop over gatherFields_
171 indexerNames_.clear();
172} // end of postRegistrationSetup() (Residual Specialization)
173
175//
176// preEvaluate() (Residual Specialization)
177//
179template<typename TRAITS, typename LO, typename GO>
180void
181panzer::
182GatherSolution_BlockedEpetra<panzer::Traits::Residual, TRAITS, LO, GO>::
183preEvaluate(
184 typename TRAITS::PreEvalData d)
185{
186 using std::logic_error;
187 using std::string;
188 using Teuchos::RCP;
189 using Teuchos::rcp_dynamic_cast;
190 using Teuchos::typeName;
195
196 // First try the refactored ReadOnly container.
197 RCP<GED> ged;
198 string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
199 if (d.gedc->containsDataObject(globalDataKey_ + post))
200 {
201 ged = d.gedc->getDataObject(globalDataKey_ + post);
202 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
203 return;
204 } // end of the refactored ReadOnly way
205
206 // Now try the old path.
207 {
208 ged = d.gedc->getDataObject(globalDataKey_);
209
210 // Extract the linear object container.
211 auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
212 auto beLoc = rcp_dynamic_cast<const BELOC>(ged);
213 if (not roGed.is_null())
214 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
215 else if (not beLoc.is_null())
216 {
217 if (useTimeDerivativeSolutionVector_)
218 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
219 else // if (not useTimeDerivativeSolutionVector_)
220 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
221 TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
222 "Residual: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
223 "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
224 typeName(*ged));
225 } // end if we have a roGed or beLoc
226 } // end of the old path
227} // end of preEvaluate() (Residual Specialization)
228
230//
231// evaluateFields() (Residual Specialization)
232//
234template<typename TRAITS, typename LO, typename GO>
235void
239 typename TRAITS::EvalData workset)
240{
241 using PHX::MDField;
242 using std::size_t;
243 using std::string;
244 using std::vector;
245 using Teuchos::ArrayRCP;
246 using Teuchos::ptrFromRef;
247 using Teuchos::RCP;
248 using Teuchos::rcp_dynamic_cast;
250 using Thyra::SpmdVectorBase;
251 using Thyra::VectorBase;
252
253 // For convenience, pull out some objects from the workset.
254 string blockId(this->wda(workset).block_id);
255 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
256 int numFields(gatherFields_.size()), numCells(localCellIds.size());
257
258 // Loop over the fields to be gathered.
259 for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
260 {
261 MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
262 auto field_h = Kokkos::create_mirror_view(field.get_static_view());
263
264 int indexerId(indexerIds_[fieldInd]),
265 subFieldNum(subFieldIds_[fieldInd]);
266
267 // Grab the local data for inputing.
268 ArrayRCP<const double> x;
269 Teuchos::RCP<const ReadOnlyVector_GlobalEvaluationData> xEvRoGed;
270
271 if(not x_.is_null()) {
272 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
273 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
274 }
275 else {
276 xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
277 }
278
279 auto subRowIndexer = indexers_[indexerId];
280 const vector<int>& elmtOffset =
281 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
282 int numBases(elmtOffset.size());
283
284 auto LIDs = subRowIndexer->getLIDs();
285 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
286 Kokkos::deep_copy(LIDs_h, LIDs);
287
288 // Gather operation for each cell in the workset.
289 for (int cell(0); cell < numCells; ++cell)
290 {
291 LO cellLocalId = localCellIds[cell];
292
293 // Loop over the basis functions and fill the fields.
294 for (int basis(0); basis < numBases; ++basis)
295 {
296 int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
297 if(x_.is_null())
298 field_h(cell, basis) = (*xEvRoGed)[lid];
299 else
300 field_h(cell, basis) = x[lid];
301 } // end loop over the basis functions
302 } // end loop over localCellIds
303 Kokkos::deep_copy(field.get_static_view(), field_h);
304 } // end loop over the fields to be gathered
305} // end of evaluateFields() (Residual Specialization)
306
308//
309// Initializing Constructor (Tangent Specialization)
310//
312template<typename TRAITS, typename LO, typename GO>
315 const std::vector<Teuchos::RCP<const GlobalIndexer>>&
316 indexers,
317 const Teuchos::ParameterList& p)
318 :
319 indexers_(indexers),
320 hasTangentFields_(false)
321{
322 using panzer::PureBasis;
323 using PHX::MDField;
324 using PHX::print;
325 using std::size_t;
326 using std::string;
327 using std::vector;
328 using Teuchos::RCP;
329 using vvstring = std::vector<std::vector<std::string>>;
331 input.setParameterList(p);
332 const vector<string>& names = input.getDofNames();
333 RCP<const PureBasis> basis = input.getBasis();
334 const vvstring& tangentFieldNames = input.getTangentNames();
335 indexerNames_ = input.getIndexerNames();
336 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
337 globalDataKey_ = input.getGlobalDataKey();
338
339 // Allocate the fields.
340 int numFields(names.size());
341 gatherFields_.resize(numFields);
342 for (int fd(0); fd < numFields; ++fd)
343 {
344 gatherFields_[fd] =
345 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
346 this->addEvaluatedField(gatherFields_[fd]);
347 } // end loop over names
348
349 // Set up the dependent tangent fields, if requested.
350 if (tangentFieldNames.size() > 0)
351 {
352 TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size());
353 hasTangentFields_ = true;
354 tangentFields_.resize(numFields);
355 for (int fd(0); fd < numFields; ++fd)
356 {
357 int numTangentFields(tangentFieldNames[fd].size());
358 tangentFields_[fd].resize(numTangentFields);
359 for (int i(0); i < numTangentFields; ++i)
360 {
361 tangentFields_[fd][i] =
362 MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
363 basis->functional);
364 this->addDependentField(tangentFields_[fd][i]);
365 } // end loop over tangentFieldNames
366 } // end loop over gatherFields_
367 } // end if we have tangent fields
368
369 // Figure out what the first active name is.
370 string firstName("<none>");
371 if (numFields > 0)
372 firstName = names[0];
373 string n("GatherSolution Tangent (BlockedEpetra): " + firstName + " (" +
374 print<EvalT>() + ")");
375 this->setName(n);
376} // end of Initializing Constructor (Tangent Specialization)
377
379//
380// postRegistrationSetup() (Tangent Specialization)
381//
383template<typename TRAITS, typename LO, typename GO>
384void
387 typename TRAITS::SetupData /* d */,
389{
390 using std::size_t;
391 using std::string;
392 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
393 int numFields(gatherFields_.size());
394 indexerIds_.resize(numFields);
395 subFieldIds_.resize(numFields);
396 for (int fd(0); fd < numFields; ++fd)
397 {
398 // Get the field ID from the DOF manager.
399 const string& fieldName(indexerNames_[fd]);
400 indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
401 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
402 TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
403 } // end loop over gatherFields_
404 indexerNames_.clear();
405} // end of postRegistrationSetup() (Tangent Specialization)
406
408//
409// preEvaluate() (Tangent Specialization)
410//
412template<typename TRAITS, typename LO, typename GO>
413void
416 typename TRAITS::PreEvalData d)
417{
418 using std::logic_error;
419 using std::string;
420 using Teuchos::RCP;
421 using Teuchos::rcp_dynamic_cast;
422 using Teuchos::typeName;
427
428 // First try the refactored ReadOnly container.
429 RCP<GED> ged;
430 string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
431 if (d.gedc->containsDataObject(globalDataKey_ + post))
432 {
433 ged = d.gedc->getDataObject(globalDataKey_ + post);
434 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
435 return;
436 } // end of the refactored ReadOnly way
437
438 // Now try the old path.
439 {
440 ged = d.gedc->getDataObject(globalDataKey_);
441
442 // Extract the linear object container.
443 auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
444 auto beLoc = rcp_dynamic_cast<const BELOC>(ged);
445 if (not roGed.is_null())
446 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
447 else if (not beLoc.is_null())
448 {
449 if (useTimeDerivativeSolutionVector_)
450 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
451 else // if (not useTimeDerivativeSolutionVector_)
452 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
453 TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
454 "Tangent: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
455 "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
456 typeName(*ged));
457 } // end if we have a roGed or beLoc
458 } // end of the old path
459} // end of preEvaluate() (Tangent Specialization)
460
462//
463// evaluateFields() (Tangent Specialization)
464//
466template<typename TRAITS, typename LO, typename GO>
467void
470 typename TRAITS::EvalData workset)
471{
472 using PHX::MDField;
473 using std::pair;
474 using std::size_t;
475 using std::string;
476 using std::vector;
477 using Teuchos::ArrayRCP;
478 using Teuchos::ptrFromRef;
479 using Teuchos::RCP;
480 using Teuchos::rcp_dynamic_cast;
482 using Thyra::SpmdVectorBase;
483 using Thyra::VectorBase;
484
485 // For convenience, pull out some objects from the workset.
486 vector<pair<int, GO>> GIDs;
487 string blockId(this->wda(workset).block_id);
488 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
489 int numFields(gatherFields_.size()), numCells(localCellIds.size());
490
491 if (x_.is_null())
492 {
493 // Loop over the fields to be gathered.
494 for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
495 {
496 MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
497 int indexerId(indexerIds_[fieldInd]),
498 subFieldNum(subFieldIds_[fieldInd]);
499
500 // Grab the local data for inputing.
501 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
502 auto subRowIndexer = indexers_[indexerId];
503 const vector<int>& elmtOffset =
504 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
505 int numBases(elmtOffset.size());
506
507 // Gather operation for each cell in the workset.
508 for (int cell(0); cell < numCells; ++cell)
509 {
510 LO cellLocalId = localCellIds[cell];
511 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
512
513 // Loop over the basis functions and fill the fields.
514 for (int basis(0); basis < numBases; ++basis)
515 {
516 int offset(elmtOffset[basis]), lid(LIDs[offset]);
517 field(cell, basis) = (*xEvRoGed)[lid];
518 } // end loop over the basis functions
519 } // end loop over localCellIds
520 } // end loop over the fields to be gathered
521 }
522 else // if (not x_.is_null())
523 {
524 // Loop over the fields to be gathered.
525 for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
526 {
527 MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
528 int indexerId(indexerIds_[fieldInd]),
529 subFieldNum(subFieldIds_[fieldInd]);
530
531 // Grab the local data for inputing.
532 ArrayRCP<const double> x;
533 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
534 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
535 auto subRowIndexer = indexers_[indexerId];
536 const vector<int>& elmtOffset =
537 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
538 int numBases(elmtOffset.size());
539
540 // Gather operation for each cell in the workset.
541 for (int cell(0); cell < numCells; ++cell)
542 {
543 LO cellLocalId = localCellIds[cell];
544 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
545
546 // Loop over the basis functions and fill the fields.
547 for (int basis(0); basis < numBases; ++basis)
548 {
549 int offset(elmtOffset[basis]), lid(LIDs[offset]);
550 field(cell, basis) = x[lid];
551 } // end loop over the basis functions
552 } // end loop over localCellIds
553 } // end loop over the fields to be gathered
554 } // end if (x_.is_null()) or not
555
556 // Deal with the tangent fields.
557 if (hasTangentFields_)
558 {
559 // Loop over the fields to be gathered.
560 for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
561 {
562 MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
563 int indexerId(indexerIds_[fieldInd]),
564 subFieldNum(subFieldIds_[fieldInd]);
565 auto subRowIndexer = indexers_[indexerId];
566 const vector<int>& elmtOffset =
567 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
568 int numBases(elmtOffset.size());
569
570 // Gather operation for each cell in the workset.
571 for (int cell(0); cell < numCells; ++cell)
572 {
573 // Loop over the basis functions and fill the fields.
574 for (int basis(0); basis < numBases; ++basis)
575 {
576 int numTangentFields(tangentFields_[fieldInd].size());
577 for (int i(0); i < numTangentFields; ++i)
578 field(cell, basis).fastAccessDx(i) =
579 tangentFields_[fieldInd][i](cell, basis).val();
580 } // end loop over the basis functions
581 } // end loop over localCellIds
582 } // end loop over the fields to be gathered
583 } // end if (hasTangentFields_)
584} // end of evaluateFields() (Tangent Specialization)
585
587//
588// Initializing Constructor (Jacobian Specialization)
589//
591template<typename TRAITS, typename LO, typename GO>
595 const std::vector<Teuchos::RCP<const GlobalIndexer>>&
596 indexers,
597 const Teuchos::ParameterList& p)
598 :
599 indexers_(indexers)
600{
601 using panzer::PureBasis;
602 using PHX::MDField;
603 using PHX::print;
604 using std::size_t;
605 using std::string;
606 using std::vector;
607 using Teuchos::RCP;
609 input.setParameterList(p);
610 const vector<string>& names = input.getDofNames();
611 RCP<const PureBasis> basis = input.getBasis();
612 indexerNames_ = input.getIndexerNames();
613 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
614 globalDataKey_ = input.getGlobalDataKey();
615 gatherSeedIndex_ = input.getGatherSeedIndex();
616 sensitivitiesName_ = input.getSensitivitiesName();
617 disableSensitivities_ = not input.firstSensitivitiesAvailable();
618
619 // Allocate the fields.
620 int numFields(names.size());
621 gatherFields_.resize(numFields);
622 for (int fd(0); fd < numFields; ++fd)
623 {
624 MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
625 gatherFields_[fd] = f;
626 this->addEvaluatedField(gatherFields_[fd]);
627 } // end loop over names
628
629 // Figure out what the first active name is.
630 string firstName("<none>"), n("GatherSolution (BlockedEpetra");
631 if (numFields > 0)
632 firstName = names[0];
633 if (disableSensitivities_)
634 n += ", No Sensitivities";
635 n += "): " + firstName + " (" + print<EvalT>() + ")";
636 this->setName(n);
637} // end of Initializing Constructor (Jacobian Specialization)
638
640//
641// postRegistrationSetup() (Jacobian Specialization)
642//
644template<typename TRAITS, typename LO, typename GO>
645void
649 typename TRAITS::SetupData /* d */,
651{
652 using std::size_t;
653 using std::string;
654 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
655 int numFields(gatherFields_.size());
656 indexerIds_.resize(numFields);
657 subFieldIds_.resize(numFields);
658 for (int fd(0); fd < numFields; ++fd)
659 {
660 // Get the field ID from the DOF manager.
661 const string& fieldName(indexerNames_[fd]);
662 indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
663 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
664 TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
665 } // end loop over gatherFields_
666 indexerNames_.clear();
667} // end of postRegistrationSetup() (Jacobian Specialization)
668
670//
671// preEvaluate() (Jacobian Specialization)
672//
674template<typename TRAITS, typename LO, typename GO>
675void
676panzer::
677GatherSolution_BlockedEpetra<panzer::Traits::Jacobian, TRAITS, LO, GO>::
678preEvaluate(
679 typename TRAITS::PreEvalData d)
680{
681 using std::endl;
682 using std::logic_error;
683 using std::string;
684 using Teuchos::RCP;
685 using Teuchos::rcp_dynamic_cast;
686 using Teuchos::typeName;
691
692 // Manage sensitivities.
693 applySensitivities_ = false;
694 if ((not disableSensitivities_ ) and
695 (d.first_sensitivities_name == sensitivitiesName_))
696 applySensitivities_ = true;
697
698 // First try the refactored ReadOnly container.
699 RCP<GED> ged;
700 string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
701 if (d.gedc->containsDataObject(globalDataKey_ + post))
702 {
703 ged = d.gedc->getDataObject(globalDataKey_ + post);
704 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
705 return;
706 } // end of the refactored ReadOnly way
707
708 // Now try the old path.
709 {
710 ged = d.gedc->getDataObject(globalDataKey_);
711
712 // Extract the linear object container.
713 auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
714 auto beLoc = rcp_dynamic_cast<const BELOC>(ged);
715 if (not roGed.is_null())
716 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
717 else if (not beLoc.is_null())
718 {
719 if (useTimeDerivativeSolutionVector_)
720 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
721 else // if (not useTimeDerivativeSolutionVector_)
722 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
723 TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
724 "Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ <<
725 "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
726 typeName(*ged));
727 } // end if we have a roGed or beLoc
728 } // end of the old path
729} // end of preEvaluate() (Jacobian Specialization)
730
732//
733// evaluateFields() (Jacobian Specialization)
734//
736template<typename TRAITS, typename LO, typename GO>
737void
741 typename TRAITS::EvalData workset)
742{
743 using PHX::MDField;
744 using std::size_t;
745 using std::string;
746 using std::vector;
747 using Teuchos::ArrayRCP;
748 using Teuchos::ptrFromRef;
749 using Teuchos::RCP;
750 using Teuchos::rcp_dynamic_cast;
752 using Thyra::SpmdVectorBase;
753 using Thyra::VectorBase;
754
755 // For convenience, pull out some objects from the workset.
756 string blockId(this->wda(workset).block_id);
757 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
758 int numFields(gatherFields_.size()), numCells(localCellIds.size());
759 double seedValue(0);
760 if (applySensitivities_)
761 {
762 if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
763 seedValue = workset.alpha;
764 else if (gatherSeedIndex_ < 0)
765 seedValue = workset.beta;
766 else if (not useTimeDerivativeSolutionVector_)
767 seedValue = workset.gather_seeds[gatherSeedIndex_];
768 else
769 TEUCHOS_ASSERT(false);
770 } // end if (applySensitivities_)
771
772 // Turn off sensitivies: This may be faster if we don't expand the term, but
773 // I suspect not, because anywhere it is used the full complement of
774 // sensitivies will be needed anyway.
775 if (not applySensitivities_)
776 seedValue = 0.0;
777
778 vector<int> blockOffsets;
779 computeBlockOffsets(blockId, indexers_, blockOffsets);
780 int numDerivs(blockOffsets[blockOffsets.size() - 1]);
781
782 // NOTE: A reordering of these loops will likely improve performance. The
783 // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
784 // can be cheaper. However the lookup for LIDs may be more expensive!
785
786 // Loop over the fields to be gathered.
787 for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
788 {
789 MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
790 auto field_h = Kokkos::create_mirror_view(field.get_view());
791
792 int indexerId(indexerIds_[fieldInd]), subFieldNum(subFieldIds_[fieldInd]);
793
794 // Grab the local data for inputing.
795 ArrayRCP<const double> x;
796 Teuchos::RCP<const ReadOnlyVector_GlobalEvaluationData> xEvRoGed;
797 if(not x_.is_null()) {
798 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
799 }else {
800 xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
801 }
802
803 auto subRowIndexer = indexers_[indexerId];
804 const vector<int>& elmtOffset =
805 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
806 int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
807
808 auto LIDs = subRowIndexer->getLIDs();
809 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
810 Kokkos::deep_copy(LIDs_h, LIDs);
811
812 // Gather operation for each cell in the workset.
813 for (int cell(0); cell < numCells; ++cell)
814 {
815 LO cellLocalId = localCellIds[cell];
816
817 // Loop over the basis functions and fill the fields.
818 for (int basis(0); basis < numBases; ++basis)
819 {
820 // Set the value and seed the FAD object.
821 int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
822 if(x_.is_null())
823 field_h(cell, basis) = ScalarT(numDerivs, (*xEvRoGed)[lid]);
824 else
825 field_h(cell, basis) = ScalarT(numDerivs, x[lid]);
826
827 field_h(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
828 } // end loop over the basis functions
829 } // end loop over localCellIds
830 Kokkos::deep_copy(field.get_static_view(), field_h);
831 } // end loop over the fields to be gathered
832} // end of evaluateFields() (Jacobian Specialization)
833
834#endif // __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
int numFields
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.
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors.
Gathers solution values from the Newton solution vector into the nodal fields of the field manager.
void evaluateFields(typename TRAITS::EvalData d)
Evaluate Fields.
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
Constructor.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Post-Registration Setup.
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)
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
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
std::string getSensitivitiesName() const
The name of the sensitivities. Enables sensitivities at "preEvaluate" time (Jacobian and Hessian)
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)
Description and data layouts associated with a particular basis.
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer > > &ugis, std::vector< int > &blockOffsets)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer > > &ugis)