Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterResidual_Tpetra_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_TPETRA_IMPL_HPP
44#define PANZER_SCATTER_RESIDUAL_TPETRA_IMPL_HPP
45
46#include "Teuchos_RCP.hpp"
47#include "Teuchos_Assert.hpp"
48
49#include "Phalanx_DataLayout.hpp"
50
52#include "Panzer_PureBasis.hpp"
57
58#include "Phalanx_DataLayout_MDALayout.hpp"
59
60#include "Teuchos_FancyOStream.hpp"
61
62#include "Tpetra_Vector.hpp"
63#include "Tpetra_Map.hpp"
64#include "Tpetra_CrsMatrix.hpp"
65
66// **********************************************************************
67// Specialization: Residual
68// **********************************************************************
69
70template<typename TRAITS,typename LO,typename GO,typename NodeT>
72ScatterResidual_Tpetra(const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
73 const Teuchos::ParameterList& p)
74 : globalIndexer_(indexer)
75 , globalDataKey_("Residual Scatter Container")
76{
77 std::string scatterName = p.get<std::string>("Scatter Name");
78 scatterHolder_ =
79 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
80
81 // get names to be evaluated
82 const std::vector<std::string>& names =
83 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
84
85 // grab map from evaluated names to field names
86 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
87
88 Teuchos::RCP<PHX::DataLayout> dl =
89 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
90
91 // build the vector of fields that this is dependent on
92 scatterFields_.resize(names.size());
93 scratch_offsets_.resize(names.size());
94 for (std::size_t eq = 0; eq < names.size(); ++eq) {
95 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
96
97 // tell the field manager that we depend on this field
98 this->addDependentField(scatterFields_[eq]);
99 }
100
101 // this is what this evaluator provides
102 this->addEvaluatedField(*scatterHolder_);
103
104 if (p.isType<std::string>("Global Data Key"))
105 globalDataKey_ = p.get<std::string>("Global Data Key");
106
107 this->setName(scatterName+" Scatter Residual");
108}
109
110// **********************************************************************
111template<typename TRAITS,typename LO,typename GO,typename NodeT>
113postRegistrationSetup(typename TRAITS::SetupData d,
115{
116 fieldIds_.resize(scatterFields_.size());
117 const Workset & workset_0 = (*d.worksets_)[0];
118 std::string blockId = this->wda(workset_0).block_id;
119
120
121 // load required field numbers for fast use
122 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
123 // get field ID from DOF manager
124 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
125 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
126
127 const std::vector<int> & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldIds_[fd]);
128 scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
129 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
130 }
131 scratch_lids_ = PHX::View<LO**>("lids",scatterFields_[0].extent(0),
132 globalIndexer_->getElementBlockGIDCount(blockId));
133
134}
135
136// **********************************************************************
137template<typename TRAITS,typename LO,typename GO,typename NodeT>
139preEvaluate(typename TRAITS::PreEvalData d)
140{
142
143 // extract linear object container
144 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
145
146 if(tpetraContainer_==Teuchos::null) {
147 // extract linear object container
148 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
149 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
150 }
151}
152
153
154// **********************************************************************
155// Specialization: Tangent
156// **********************************************************************
157
158template<typename TRAITS,typename LO,typename GO,typename NodeT>
160ScatterResidual_Tpetra(const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
161 const Teuchos::ParameterList& p)
162 : globalIndexer_(indexer)
163 , globalDataKey_("Residual Scatter Container")
164{
165 std::string scatterName = p.get<std::string>("Scatter Name");
166 scatterHolder_ =
167 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
168
169 // get names to be evaluated
170 const std::vector<std::string>& names =
171 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
172
173 // grab map from evaluated names to field names
174 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
175
176 Teuchos::RCP<PHX::DataLayout> dl =
177 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
178
179 // build the vector of fields that this is dependent on
180 scatterFields_.resize(names.size());
181 for (std::size_t eq = 0; eq < names.size(); ++eq) {
182 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
183
184 // tell the field manager that we depend on this field
185 this->addDependentField(scatterFields_[eq]);
186 }
187
188 // this is what this evaluator provides
189 this->addEvaluatedField(*scatterHolder_);
190
191 if (p.isType<std::string>("Global Data Key"))
192 globalDataKey_ = p.get<std::string>("Global Data Key");
193
194 this->setName(scatterName+" Scatter Tangent");
195}
196
197// **********************************************************************
198template<typename TRAITS,typename LO,typename GO,typename NodeT>
200postRegistrationSetup(typename TRAITS::SetupData /* d */,
202{
203 fieldIds_.resize(scatterFields_.size());
204 // load required field numbers for fast use
205 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
206 // get field ID from DOF manager
207 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
208 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
209 }
210}
211
212// **********************************************************************
213template<typename TRAITS,typename LO,typename GO,typename NodeT>
215preEvaluate(typename TRAITS::PreEvalData d)
216{
217 using Teuchos::RCP;
218 using Teuchos::rcp_dynamic_cast;
219
221
222 // this is the list of parameters and their names that this scatter has to account for
223 std::vector<std::string> activeParameters =
224 rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
225
226 dfdp_vectors_.clear();
227 for(std::size_t i=0;i<activeParameters.size();i++) {
228 RCP<typename LOC::VectorType> vec =
229 rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
230 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
231 dfdp_vectors_.push_back(vec_array);
232 }
233}
234
235// **********************************************************************
236template<typename TRAITS,typename LO,typename GO,typename NodeT>
238evaluateFields(typename TRAITS::EvalData workset)
239{
240 // for convenience pull out some objects from workset
241 std::string blockId = this->wda(workset).block_id;
242 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
243
244 // NOTE: A reordering of these loops will likely improve performance
245 // The "getGIDFieldOffsets may be expensive. However the
246 // "getElementGIDs" can be cheaper. However the lookup for LIDs
247 // may be more expensive!
248
249 // scatter operation for each cell in workset
250 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
251 std::size_t cellLocalId = localCellIds[worksetCellIndex];
252
253 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
254
255 // loop over each field to be scattered
256 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
257 int fieldNum = fieldIds_[fieldIndex];
258 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
259
260 // loop over basis functions
261 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
262 int offset = elmtOffset[basis];
263 LO lid = LIDs[offset];
264 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
265 for(int d=0;d<value.size();d++)
266 dfdp_vectors_[d][lid] += value.fastAccessDx(d);
267 }
268 }
269 }
270}
271
272// **********************************************************************
273// Specialization: Jacobian
274// **********************************************************************
275
276template<typename TRAITS,typename LO,typename GO,typename NodeT>
278ScatterResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
279 const Teuchos::ParameterList& p)
280 : globalIndexer_(indexer)
281 , globalDataKey_("Residual Scatter Container")
282 , my_derivative_size_(0)
283 , other_derivative_size_(0)
284{
285 std::string scatterName = p.get<std::string>("Scatter Name");
286 scatterHolder_ =
287 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
288
289 // get names to be evaluated
290 const std::vector<std::string>& names =
291 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
292
293 // grab map from evaluated names to field names
294 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
295
296 Teuchos::RCP<PHX::DataLayout> dl =
297 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
298
299 // build the vector of fields that this is dependent on
300 scatterFields_.resize(names.size());
301 scratch_offsets_.resize(names.size());
302 for (std::size_t eq = 0; eq < names.size(); ++eq) {
303 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
304
305 // tell the field manager that we depend on this field
306 this->addDependentField(scatterFields_[eq]);
307 }
308
309 // this is what this evaluator provides
310 this->addEvaluatedField(*scatterHolder_);
311
312 if (p.isType<std::string>("Global Data Key"))
313 globalDataKey_ = p.get<std::string>("Global Data Key");
314
315 this->setName(scatterName+" Scatter Residual (Jacobian)");
316}
317
318// **********************************************************************
319template<typename TRAITS,typename LO,typename GO,typename NodeT>
321postRegistrationSetup(typename TRAITS::SetupData d,
323{
324 fieldIds_.resize(scatterFields_.size());
325
326 const Workset & workset_0 = (*d.worksets_)[0];
327 std::string blockId = this->wda(workset_0).block_id;
328
329 // load required field numbers for fast use
330 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
331 // get field ID from DOF manager
332 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
333 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
334
335 int fieldNum = fieldIds_[fd];
336 const std::vector<int> & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
337 scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
338 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
339 }
340
341 my_derivative_size_ = globalIndexer_->getElementBlockGIDCount(blockId);
342 if (Teuchos::nonnull(workset_0.other)) {
343 auto otherBlockId = workset_0.other->block_id;
344 other_derivative_size_ = globalIndexer_->getElementBlockGIDCount(otherBlockId);
345 }
346 scratch_lids_ = Kokkos::View<LO**, Kokkos::LayoutRight, PHX::Device>(
347 "lids", scatterFields_[0].extent(0), my_derivative_size_ + other_derivative_size_ );
348 scratch_vals_ = Kokkos::View<typename Sacado::ScalarType<ScalarT>::type**, Kokkos::LayoutRight, PHX::Device>(
349 "vals", scatterFields_[0].extent(0), my_derivative_size_ + other_derivative_size_ );
350}
351
352// **********************************************************************
353template<typename TRAITS,typename LO,typename GO,typename NodeT>
355preEvaluate(typename TRAITS::PreEvalData d)
356{
358
359 // extract linear object container
360 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
361
362 if(tpetraContainer_==Teuchos::null) {
363 // extract linear object container
364 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
365 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
366 }
367}
368
369
370// **********************************************************************
371namespace panzer {
372namespace {
373
374template <typename ScalarT,typename LO,typename GO,typename NodeT,typename LocalMatrixT>
375class ScatterResidual_Jacobian_Functor {
376public:
377 typedef typename PHX::Device execution_space;
378 typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
379
381 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
382 LocalMatrixT jac; // Kokkos jacobian type
383
384 Kokkos::View<const LO**, Kokkos::LayoutRight, PHX::Device> lids; // local indices for unknowns.
385 Kokkos::View<typename Sacado::ScalarType<ScalarT>::type**, Kokkos::LayoutRight, PHX::Device> vals;
386 PHX::View<const int*> offsets; // how to get a particular field
387 FieldType field;
388
389
390 KOKKOS_INLINE_FUNCTION
391 void operator()(const unsigned int cell) const
392 {
393 int numIds = lids.extent(1);
394
395 // loop over the basis functions (currently they are nodes)
396 for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
397 typename FieldType::array_type::reference_type scatterField = field(cell,basis);
398 int offset = offsets(basis);
399 LO lid = lids(cell,offset);
400
401 // Sum residual
402 if(fillResidual)
403 Kokkos::atomic_add(&r_data(lid,0), scatterField.val());
404
405 // loop over the sensitivity indices: all DOFs on a cell
406 for(int sensIndex=0;sensIndex<numIds;++sensIndex)
407 vals(cell,sensIndex) = scatterField.fastAccessDx(sensIndex);
408
409 // Sum Jacobian
410 jac.sumIntoValues(lid, &lids(cell,0), numIds, &vals(cell,0), true, true);
411 } // end basis
412 }
413};
414
415template <typename ScalarT,typename LO,typename GO,typename NodeT>
416class ScatterResidual_Residual_Functor {
417public:
418 typedef typename PHX::Device execution_space;
419 typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
420
421 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
422
423 PHX::View<const LO**> lids; // local indices for unknowns
424 PHX::View<const int*> offsets; // how to get a particular field
426
427
428 KOKKOS_INLINE_FUNCTION
429 void operator()(const unsigned int cell) const
430 {
431
432 // loop over the basis functions (currently they are nodes)
433 for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
434 int offset = offsets(basis);
435 LO lid = lids(cell,offset);
436 Kokkos::atomic_add(&r_data(lid,0), field(cell,basis));
437
438 } // end basis
439 }
440};
441
442}
443}
444
445// **********************************************************************
446template<typename TRAITS,typename LO,typename GO,typename NodeT>
448evaluateFields(typename TRAITS::EvalData workset)
449{
451
452 // for convenience pull out some objects from workset
453 std::string blockId = this->wda(workset).block_id;
454
455 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
456
457 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
458
459 ScatterResidual_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
460 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
461 functor.lids = scratch_lids_;
462
463 // for each field, do a parallel for loop
464 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
465 functor.offsets = scratch_offsets_[fieldIndex];
466 functor.field = scatterFields_[fieldIndex];
467
468 Kokkos::parallel_for(workset.num_cells,functor);
469 }
470}
471
472// **********************************************************************
473template<typename TRAITS,typename LO,typename GO,typename NodeT>
475evaluateFields(typename TRAITS::EvalData workset)
476{
478
479 typedef typename LOC::CrsMatrixType::local_matrix_device_type LocalMatrixT;
480
481 // for convenience pull out some objects from workset
482 std::string blockId = this->wda(workset).block_id;
483
484 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
485 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
486
487 // Cache scratch lids. For interface bc problems the derivative
488 // dimension extent spans two cells. Use subviews to get the self
489 // lids and the other lids.
490 if (Teuchos::nonnull(workset.other)) {
491 auto my_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(0,my_derivative_size_));
492 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,my_scratch_lids);
493 auto other_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(my_derivative_size_,my_derivative_size_ + other_derivative_size_));
494 globalIndexer_->getElementLIDs(workset.other->cell_local_ids_k,other_scratch_lids);
495 }
496 else {
497 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
498 }
499
500 ScatterResidual_Jacobian_Functor<ScalarT,LO,GO,NodeT,LocalMatrixT> functor;
501 functor.fillResidual = (r!=Teuchos::null);
502 if(functor.fillResidual)
503 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
504 functor.jac = Jac->getLocalMatrixDevice();
505 functor.lids = scratch_lids_;
506 functor.vals = scratch_vals_;
507
508 // for each field, do a parallel for loop
509 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
510 functor.offsets = scratch_offsets_[fieldIndex];
511 functor.field = scatterFields_[fieldIndex];
512
513 Kokkos::parallel_for(workset.num_cells,functor);
514 }
515
516}
517
518// **********************************************************************
519
520#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.
Kokkos::View< const LO **, Kokkos::LayoutRight, PHX::Device > lids
Kokkos::View< typename Sacado::ScalarType< ScalarT >::type **, Kokkos::LayoutRight, PHX::Device > vals
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > r_data
Pushes residual values into the residual vector for a Newton-based solve.
const Teuchos::RCP< VectorType > get_f() const
std::string block_id
DEPRECATED - use: getElementBlock()
Teuchos::RCP< WorksetDetails > other
FieldType
The type of discretization to use for a field pattern.