Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterResidual_Epetra_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_EPETRA_IMPL_HPP
44#define PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP
45
46#include "Teuchos_RCP.hpp"
47#include "Teuchos_Assert.hpp"
48
49#include "Phalanx_DataLayout.hpp"
50
51#include "Epetra_Map.h"
52#include "Epetra_Vector.h"
53#include "Epetra_CrsMatrix.h"
54
56#include "Panzer_PureBasis.hpp"
62
63#include "Phalanx_DataLayout_MDALayout.hpp"
64
65#include "Teuchos_FancyOStream.hpp"
66
67// **********************************************************************
68// Specialization: Residual
69// **********************************************************************
70
71template<typename TRAITS,typename LO,typename GO>
73ScatterResidual_Epetra(const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
74 const Teuchos::RCP<const panzer::GlobalIndexer> & /* cIndexer */,
75 const Teuchos::ParameterList& p,
76 bool useDiscreteAdjoint)
77 : globalIndexer_(indexer)
78 , globalDataKey_("Residual Scatter Container")
79 , useDiscreteAdjoint_(useDiscreteAdjoint)
80{
81 std::string scatterName = p.get<std::string>("Scatter Name");
82 scatterHolder_ =
83 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
84
85 // get names to be evaluated
86 const std::vector<std::string>& names =
87 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
88
89 // grab map from evaluated names to field names
90 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
91
92 Teuchos::RCP<PHX::DataLayout> dl =
93 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
94
95 // build the vector of fields that this is dependent on
96 scatterFields_.resize(names.size());
97 for (std::size_t eq = 0; eq < names.size(); ++eq) {
98 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
99
100 // tell the field manager that we depend on this field
101 this->addDependentField(scatterFields_[eq]);
102 }
103
104 // this is what this evaluator provides
105 this->addEvaluatedField(*scatterHolder_);
106
107 if (p.isType<std::string>("Global Data Key"))
108 globalDataKey_ = p.get<std::string>("Global Data Key");
109
110 this->setName(scatterName+" Scatter Residual");
111}
112
113// **********************************************************************
114template<typename TRAITS,typename LO,typename GO>
116postRegistrationSetup(typename TRAITS::SetupData /* d */,
118{
119 fieldIds_.resize(scatterFields_.size());
120 // load required field numbers for fast use
121 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
122 // get field ID from DOF manager
123 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
124 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
125 }
126}
127
128// **********************************************************************
129template<typename TRAITS,typename LO,typename GO>
131preEvaluate(typename TRAITS::PreEvalData d)
132{
133 // extract linear object container
134 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
135
136 if(epetraContainer_==Teuchos::null) {
137 // extract linear object container
138 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
139 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
140 }
141}
142
143// **********************************************************************
144template<typename TRAITS,typename LO,typename GO>
146evaluateFields(typename TRAITS::EvalData workset)
147{
148 // for convenience pull out some objects from workset
149 std::string blockId = this->wda(workset).block_id;
150 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
151
152 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
153
154 // NOTE: A reordering of these loops will likely improve performance
155 // The "getGIDFieldOffsets may be expensive. However the
156 // "getElementGIDs" can be cheaper. However the lookup for LIDs
157 // may be more expensive!
158
159 // scatter operation for each cell in workset
160 auto LIDs = globalIndexer_->getLIDs();
161 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
162 Kokkos::deep_copy(LIDs_h, LIDs);
163 // loop over each field to be scattered
164 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
165 int fieldNum = fieldIds_[fieldIndex];
166 auto field = PHX::as_view(scatterFields_[fieldIndex]);
167 auto field_h = Kokkos::create_mirror_view(field);
168 Kokkos::deep_copy(field_h, field);
169
170 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
171 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
172 std::size_t cellLocalId = localCellIds[worksetCellIndex];
173
174 // loop over basis functions
175 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
176 int offset = elmtOffset[basis];
177 int lid = LIDs_h(cellLocalId, offset);
178 (*r)[lid] += field_h(worksetCellIndex,basis);
179 }
180 }
181 }
182}
183
184// **********************************************************************
185// Specialization: Tangent
186// **********************************************************************
187
188template<typename TRAITS,typename LO,typename GO>
190ScatterResidual_Epetra(const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
191 const Teuchos::RCP<const panzer::GlobalIndexer> & /* cIndexer */,
192 const Teuchos::ParameterList& p,
193 bool /* useDiscreteAdjoint */)
194 : globalIndexer_(indexer)
195{
196 std::string scatterName = p.get<std::string>("Scatter Name");
197 scatterHolder_ =
198 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
199
200 // get names to be evaluated
201 const std::vector<std::string>& names =
202 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
203
204 // grab map from evaluated names to field names
205 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
206
207 Teuchos::RCP<PHX::DataLayout> dl =
208 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
209
210 // build the vector of fields that this is dependent on
211 scatterFields_.resize(names.size());
212 for (std::size_t eq = 0; eq < names.size(); ++eq) {
213 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
214
215 // tell the field manager that we depend on this field
216 this->addDependentField(scatterFields_[eq]);
217 }
218
219 // this is what this evaluator provides
220 this->addEvaluatedField(*scatterHolder_);
221
222 this->setName(scatterName+" Scatter Tangent");
223}
224
225// **********************************************************************
226template<typename TRAITS,typename LO,typename GO>
228postRegistrationSetup(typename TRAITS::SetupData /* d */,
230{
231 fieldIds_.resize(scatterFields_.size());
232 // load required field numbers for fast use
233 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
234 // get field ID from DOF manager
235 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
236 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
237 }
238}
239
240// **********************************************************************
241template<typename TRAITS,typename LO,typename GO>
243preEvaluate(typename TRAITS::PreEvalData d)
244{
245 using Teuchos::RCP;
246 using Teuchos::rcp_dynamic_cast;
247
248 // this is the list of parameters and their names that this scatter has to account for
249 std::vector<std::string> activeParameters =
250 rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
251
252 dfdp_vectors_.clear();
253 for(std::size_t i=0;i<activeParameters.size();i++) {
254 RCP<Epetra_Vector> vec = rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
255 dfdp_vectors_.push_back(vec);
256 }
257}
258
259// **********************************************************************
260template<typename TRAITS,typename LO,typename GO>
262evaluateFields(typename TRAITS::EvalData workset)
263{
264 // for convenience pull out some objects from workset
265 std::string blockId = this->wda(workset).block_id;
266 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
267
268 // NOTE: A reordering of these loops will likely improve performance
269 // The "getGIDFieldOffsets may be expensive. However the
270 // "getElementGIDs" can be cheaper. However the lookup for LIDs
271 // may be more expensive!
272
273 // loop over each field to be scattered
274 auto LIDs = globalIndexer_->getLIDs();
275 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
276 Kokkos::deep_copy(LIDs_h, LIDs);
277 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
278 int fieldNum = fieldIds_[fieldIndex];
279 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
280 auto scatterField_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
281 Kokkos::deep_copy(scatterField_h, scatterFields_[fieldIndex].get_static_view());
282 // scatter operation for each cell in workset
283 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
284 std::size_t cellLocalId = localCellIds[worksetCellIndex];
285
286 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
287 // loop over basis functions
288 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
289 int offset = elmtOffset[basis];
290 int lid = LIDs_h(cellLocalId, offset);
291
292 // scatter the sensitivity vectors
293 ScalarT value = scatterField_h(worksetCellIndex,basis);
294 for(int d=0;d<value.size();d++)
295 (*dfdp_vectors_[d])[lid] += value.fastAccessDx(d);
296 }
297 }
298 }
299}
300
301// **********************************************************************
302// Specialization: Jacobian
303// **********************************************************************
304
305template<typename TRAITS,typename LO,typename GO>
307ScatterResidual_Epetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
308 const Teuchos::RCP<const panzer::GlobalIndexer> & cIndexer,
309 const Teuchos::ParameterList& p,
310 bool useDiscreteAdjoint)
311 : globalIndexer_(indexer)
312 , colGlobalIndexer_(cIndexer)
313 , globalDataKey_("Residual Scatter Container")
314 , useDiscreteAdjoint_(useDiscreteAdjoint)
315{
316 std::string scatterName = p.get<std::string>("Scatter Name");
317 scatterHolder_ =
318 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
319
320 // get names to be evaluated
321 const std::vector<std::string>& names =
322 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
323
324 // grab map from evaluated names to field names
325 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
326
327 Teuchos::RCP<PHX::DataLayout> dl =
328 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
329
330 // build the vector of fields that this is dependent on
331 scatterFields_.resize(names.size());
332 for (std::size_t eq = 0; eq < names.size(); ++eq) {
333 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
334
335 // tell the field manager that we depend on this field
336 this->addDependentField(scatterFields_[eq]);
337 }
338
339 // this is what this evaluator provides
340 this->addEvaluatedField(*scatterHolder_);
341
342 if (p.isType<std::string>("Global Data Key"))
343 globalDataKey_ = p.get<std::string>("Global Data Key");
344 if (p.isType<bool>("Use Discrete Adjoint"))
345 useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
346
347 // discrete adjoint does not work with non-square matrices
348 if(useDiscreteAdjoint)
349 { TEUCHOS_ASSERT(colGlobalIndexer_==globalIndexer_); }
350
351 this->setName(scatterName+" Scatter Residual Epetra (Jacobian)");
352}
353
354// **********************************************************************
355template<typename TRAITS,typename LO,typename GO>
357postRegistrationSetup(typename TRAITS::SetupData /* d */,
359{
360 // globalIndexer_ = d.globalIndexer_;
361
362 fieldIds_.resize(scatterFields_.size());
363 // load required field numbers for fast use
364 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
365 // get field ID from DOF manager
366 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
367 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
368 }
369}
370
371// **********************************************************************
372template<typename TRAITS,typename LO,typename GO>
374preEvaluate(typename TRAITS::PreEvalData d)
375{
376 // extract linear object container
377 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
378
379 if(epetraContainer_==Teuchos::null) {
380 // extract linear object container
381 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
382 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
383 }
384}
385
386// **********************************************************************
387template<typename TRAITS,typename LO,typename GO>
389evaluateFields(typename TRAITS::EvalData workset)
390{
391 std::vector<double> jacRow;
392
393 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
394
395 // for convenience pull out some objects from workset
396 std::string blockId = this->wda(workset).block_id;
397 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
398
399 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
400 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
401
402 const Teuchos::RCP<const panzer::GlobalIndexer>&
403 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
404
405 // NOTE: A reordering of these loops will likely improve performance
406 // The "getGIDFieldOffsets" may be expensive. However the
407 // "getElementGIDs" can be cheaper. However the lookup for LIDs
408 // may be more expensive!
409
410 // scatter operation for each cell in workset
411 auto LIDs = globalIndexer_->getLIDs();
412 auto colLIDs = colGlobalIndexer->getLIDs();
413 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
414 auto colLIDs_h = Kokkos::create_mirror_view(colLIDs);
415 Kokkos::deep_copy(LIDs_h, LIDs);
416 Kokkos::deep_copy(colLIDs_h, colLIDs);
417 std::vector<typename decltype(scatterFields_[0].get_static_view())::HostMirror> scatterFields_h;
418 for ( std::size_t i=0; i< scatterFields_.size(); ++i) {
419 scatterFields_h.push_back(Kokkos::create_mirror_view(scatterFields_[i].get_static_view()));
420 Kokkos::deep_copy(scatterFields_h[i], scatterFields_[i].get_static_view());
421 }
422
423 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
424 std::size_t cellLocalId = localCellIds[worksetCellIndex];
425
426 std::vector<int> cLIDs;
427 for (int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
428 cLIDs.push_back(colLIDs_h(cellLocalId, i));
429 if (Teuchos::nonnull(workset.other)) {
430 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
431 for (int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
432 cLIDs.push_back(colLIDs_h(other_cellLocalId, i));
433 }
434
435 // loop over each field to be scattered
436 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
437 int fieldNum = fieldIds_[fieldIndex];
438 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
439
440 // loop over the basis functions (currently they are nodes)
441 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
442 const ScalarT scatterField = (scatterFields_h[fieldIndex])(worksetCellIndex,rowBasisNum);
443 int rowOffset = elmtOffset[rowBasisNum];
444 int row = LIDs_h(cellLocalId,rowOffset);
445
446 // Sum residual
447 if(r!=Teuchos::null)
448 r->SumIntoMyValue(row,0,scatterField.val());
449
450 // Sum Jacobian
451 if(useDiscreteAdjoint_) {
452 // loop over the sensitivity indices: all DOFs on a cell
453 jacRow.resize(scatterField.size());
454
455 for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
456 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
457
458 for(std::size_t c=0;c<cLIDs.size();c++) {
459 int err = Jac->SumIntoMyValues(cLIDs[c], 1, &jacRow[c],&row);
460 TEUCHOS_ASSERT_EQUALITY(err,0);
461 }
462 }
463 else {
464 int err = Jac->SumIntoMyValues(
465 row,
466 std::min(cLIDs.size(), static_cast<size_t>(scatterField.size())),
467 scatterField.dx(),
469
470 TEUCHOS_ASSERT_EQUALITY(err,0);
471 }
472 } // end rowBasisNum
473 } // end fieldIndex
474 }
475}
476
477// **********************************************************************
478
479#endif
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.
T * ptrFromStlVector(std::vector< T > &v)