Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterDirichletResidual_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_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP
44#define PANZER_SCATTER_DIRICHLET_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// **********************************************************************
63// Specialization: Residual
64// **********************************************************************
65
66
67template<typename TRAITS,typename LO,typename GO,typename NodeT>
69ScatterDirichletResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
70 const Teuchos::ParameterList& p)
71 : globalIndexer_(indexer)
72 , globalDataKey_("Residual Scatter Container")
73{
74 std::string scatterName = p.get<std::string>("Scatter Name");
75 scatterHolder_ =
76 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
77
78 // get names to be evaluated
79 const std::vector<std::string>& names =
80 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
81
82 // grab map from evaluated names to field names
83 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
84
85 // determine if we are scattering an initial condition
86 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
87
88 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
89 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
90 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional ;
91 if (!scatterIC_) {
92 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
93 local_side_id_ = p.get<int>("Local Side ID");
94 }
95
96 // build the vector of fields that this is dependent on
97 scatterFields_.resize(names.size());
98 for (std::size_t eq = 0; eq < names.size(); ++eq) {
99 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
100
101 // tell the field manager that we depend on this field
102 this->addDependentField(scatterFields_[eq]);
103 }
104
105 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
106 if (checkApplyBC_) {
107 applyBC_.resize(names.size());
108 for (std::size_t eq = 0; eq < names.size(); ++eq) {
109 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
110 this->addDependentField(applyBC_[eq]);
111 }
112 }
113
114 // this is what this evaluator provides
115 this->addEvaluatedField(*scatterHolder_);
116
117 if (p.isType<std::string>("Global Data Key"))
118 globalDataKey_ = p.get<std::string>("Global Data Key");
119
120 this->setName(scatterName+" Scatter Residual");
121}
122
123// **********************************************************************
124template<typename TRAITS,typename LO,typename GO,typename NodeT>
126postRegistrationSetup(typename TRAITS::SetupData /* d */,
128{
129 fieldIds_.resize(scatterFields_.size());
130
131 // load required field numbers for fast use
132 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
133 // get field ID from DOF manager
134 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
135 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
136 }
137
138 // get the number of nodes (Should be renamed basis)
139 num_nodes = scatterFields_[0].extent(1);
140}
141
142// **********************************************************************
143template<typename TRAITS,typename LO,typename GO,typename NodeT>
145preEvaluate(typename TRAITS::PreEvalData d)
146{
147 // extract linear object container
148 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
149
150 if(tpetraContainer_==Teuchos::null) {
151 // extract linear object container
152 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
153 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
154
155 dirichletCounter_ = Teuchos::null;
156 }
157 else {
158 // extract dirichlet counter from container
159 Teuchos::RCP<LOC> tpetraContainer
160 = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
161
162 dirichletCounter_ = tpetraContainer->get_f();
163 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
164 }
165}
166
167// **********************************************************************
168template<typename TRAITS,typename LO,typename GO,typename NodeT>
170evaluateFields(typename TRAITS::EvalData workset)
171{
172 std::vector<GO> GIDs;
173 std::vector<LO> LIDs;
174
175 // for convenience pull out some objects from workset
176 std::string blockId = this->wda(workset).block_id;
177 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
178
179 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
180 tpetraContainer_->get_f() :
181 tpetraContainer_->get_x();
182
183 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
184 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
185
186 // NOTE: A reordering of these loops will likely improve performance
187 // The "getGIDFieldOffsets may be expensive. However the
188 // "getElementGIDs" can be cheaper. However the lookup for LIDs
189 // may be more expensive!
190
191
192
193 // loop over each field to be scattered
194 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
195 int fieldNum = fieldIds_[fieldIndex];
196 auto scatterFields_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
197 Kokkos::deep_copy(scatterFields_h, scatterFields_[fieldIndex].get_static_view());
198
199 // scatter operation for each cell in workset
200 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
201 std::size_t cellLocalId = localCellIds[worksetCellIndex];
202
203 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
204
205 // caculate the local IDs for this element
206 LIDs.resize(GIDs.size());
207 for(std::size_t i=0;i<GIDs.size();i++)
208 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
209
210 if (!scatterIC_) {
211 // this call "should" get the right ordering according to the Intrepid2 basis
212 const std::pair<std::vector<int>,std::vector<int> > & indicePair
213 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
214 const std::vector<int> & elmtOffset = indicePair.first;
215 const std::vector<int> & basisIdMap = indicePair.second;
216
217 // loop over basis functions
218 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
219 int offset = elmtOffset[basis];
220 LO lid = LIDs[offset];
221 if(lid<0) // not on this processor!
222 continue;
223
224 int basisId = basisIdMap[basis];
225
226 if (checkApplyBC_)
227 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
228 continue;
229
230 r_array[lid] = scatterFields_h(worksetCellIndex,basisId);
231
232 // record that you set a dirichlet condition
233 dc_array[lid] = 1.0;
234 }
235 } else {
236 // this call "should" get the right ordering according to the Intrepid2 basis
237 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
238
239 // loop over basis functions
240 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
241 int offset = elmtOffset[basis];
242 LO lid = LIDs[offset];
243 if(lid<0) // not on this processor!
244 continue;
245
246 r_array[lid] = scatterFields_h(worksetCellIndex,basis);
247
248 // record that you set a dirichlet condition
249 dc_array[lid] = 1.0;
250 }
251 }
252 }
253 }
254}
255
256// **********************************************************************
257// Specialization: Tangent
258// **********************************************************************
259
260
261template<typename TRAITS,typename LO,typename GO,typename NodeT>
263ScatterDirichletResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
264 const Teuchos::ParameterList& p)
265 : globalIndexer_(indexer)
266 , globalDataKey_("Residual Scatter Container")
267{
268 std::string scatterName = p.get<std::string>("Scatter Name");
269 scatterHolder_ =
270 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
271
272 // get names to be evaluated
273 const std::vector<std::string>& names =
274 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
275
276 // grab map from evaluated names to field names
277 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
278
279 // determine if we are scattering an initial condition
280 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
281
282 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
283 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
284 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional ;
285 if (!scatterIC_) {
286 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
287 local_side_id_ = p.get<int>("Local Side ID");
288 }
289
290 // build the vector of fields that this is dependent on
291 scatterFields_.resize(names.size());
292 for (std::size_t eq = 0; eq < names.size(); ++eq) {
293 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
294
295 // tell the field manager that we depend on this field
296 this->addDependentField(scatterFields_[eq]);
297 }
298
299 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
300 if (checkApplyBC_) {
301 applyBC_.resize(names.size());
302 for (std::size_t eq = 0; eq < names.size(); ++eq) {
303 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
304 this->addDependentField(applyBC_[eq]);
305 }
306 }
307
308 // this is what this evaluator provides
309 this->addEvaluatedField(*scatterHolder_);
310
311 if (p.isType<std::string>("Global Data Key"))
312 globalDataKey_ = p.get<std::string>("Global Data Key");
313
314 this->setName(scatterName+" Scatter Tangent");
315}
316
317// **********************************************************************
318template<typename TRAITS,typename LO,typename GO,typename NodeT>
320postRegistrationSetup(typename TRAITS::SetupData /* d */,
322{
323 fieldIds_.resize(scatterFields_.size());
324
325 // load required field numbers for fast use
326 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
327 // get field ID from DOF manager
328 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
329 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
330 }
331
332 // get the number of nodes (Should be renamed basis)
333 num_nodes = scatterFields_[0].extent(1);
334}
335
336// **********************************************************************
337template<typename TRAITS,typename LO,typename GO,typename NodeT>
339preEvaluate(typename TRAITS::PreEvalData d)
340{
341 // extract linear object container
342 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
343
344 if(tpetraContainer_==Teuchos::null) {
345 // extract linear object container
346 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
347 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
348
349 dirichletCounter_ = Teuchos::null;
350 }
351 else {
352 // extract dirichlet counter from container
353 Teuchos::RCP<LOC> tpetraContainer
354 = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
355
356 dirichletCounter_ = tpetraContainer->get_f();
357 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
358 }
359
360 using Teuchos::RCP;
361 using Teuchos::rcp_dynamic_cast;
362
363 // this is the list of parameters and their names that this scatter has to account for
364 std::vector<std::string> activeParameters =
365 rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
366
367 // ETP 02/03/16: This code needs to be updated to properly handle scatterIC_
368 TEUCHOS_ASSERT(!scatterIC_);
369 dfdp_vectors_.clear();
370 for(std::size_t i=0;i<activeParameters.size();i++) {
371 RCP<typename LOC::VectorType> vec =
372 rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
373 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
374 dfdp_vectors_.push_back(vec_array);
375 }
376}
377
378// **********************************************************************
379template<typename TRAITS,typename LO,typename GO,typename NodeT>
381evaluateFields(typename TRAITS::EvalData workset)
382{
383 std::vector<GO> GIDs;
384 std::vector<LO> LIDs;
385
386 // for convenience pull out some objects from workset
387 std::string blockId = this->wda(workset).block_id;
388 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
389
390 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
391 tpetraContainer_->get_f() :
392 tpetraContainer_->get_x();
393
394 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
395 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
396
397 // NOTE: A reordering of these loops will likely improve performance
398 // The "getGIDFieldOffsets may be expensive. However the
399 // "getElementGIDs" can be cheaper. However the lookup for LIDs
400 // may be more expensive!
401
402
403 // scatter operation for each cell in workset
404 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
405 std::size_t cellLocalId = localCellIds[worksetCellIndex];
406
407 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
408
409 // caculate the local IDs for this element
410 LIDs.resize(GIDs.size());
411 for(std::size_t i=0;i<GIDs.size();i++)
412 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
413
414 // loop over each field to be scattered
415 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
416 int fieldNum = fieldIds_[fieldIndex];
417
418 if (!scatterIC_) {
419 // this call "should" get the right ordering according to the Intrepid2 basis
420 const std::pair<std::vector<int>,std::vector<int> > & indicePair
421 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
422 const std::vector<int> & elmtOffset = indicePair.first;
423 const std::vector<int> & basisIdMap = indicePair.second;
424
425 // loop over basis functions
426 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
427 int offset = elmtOffset[basis];
428 LO lid = LIDs[offset];
429 if(lid<0) // not on this processor!
430 continue;
431
432 int basisId = basisIdMap[basis];
433
434 if (checkApplyBC_)
435 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
436 continue;
437
438 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
439 //r_array[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
440
441 // then scatter the sensitivity vectors
442 if(value.size()==0)
443 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
444 dfdp_vectors_[d][lid] = 0.0;
445 else
446 for(int d=0;d<value.size();d++) {
447 dfdp_vectors_[d][lid] = value.fastAccessDx(d);
448 }
449
450 // record that you set a dirichlet condition
451 dc_array[lid] = 1.0;
452 }
453 } else {
454 // this call "should" get the right ordering according to the Intrepid2 basis
455 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
456
457 // loop over basis functions
458 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
459 int offset = elmtOffset[basis];
460 LO lid = LIDs[offset];
461 if(lid<0) // not on this processor!
462 continue;
463
464 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
465 //r_array[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
466
467 // then scatter the sensitivity vectors
468 if(value.size()==0)
469 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
470 dfdp_vectors_[d][lid] = 0.0;
471 else
472 for(int d=0;d<value.size();d++) {
473 dfdp_vectors_[d][lid] = value.fastAccessDx(d);
474 }
475
476 // record that you set a dirichlet condition
477 dc_array[lid] = 1.0;
478 }
479 }
480 }
481 }
482}
483
484// **********************************************************************
485// Specialization: Jacobian
486// **********************************************************************
487
488template<typename TRAITS,typename LO,typename GO,typename NodeT>
490ScatterDirichletResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
491 const Teuchos::ParameterList& p)
492 : globalIndexer_(indexer)
493 , globalDataKey_("Residual Scatter Container")
494{
495 std::string scatterName = p.get<std::string>("Scatter Name");
496 scatterHolder_ =
497 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
498
499 // get names to be evaluated
500 const std::vector<std::string>& names =
501 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
502
503 // grab map from evaluated names to field names
504 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
505
506 Teuchos::RCP<PHX::DataLayout> dl =
507 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
508
509 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
510 local_side_id_ = p.get<int>("Local Side ID");
511
512 // build the vector of fields that this is dependent on
513 scatterFields_.resize(names.size());
514 for (std::size_t eq = 0; eq < names.size(); ++eq) {
515 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
516
517 // tell the field manager that we depend on this field
518 this->addDependentField(scatterFields_[eq]);
519 }
520
521 checkApplyBC_ = p.get<bool>("Check Apply BC");
522 if (checkApplyBC_) {
523 applyBC_.resize(names.size());
524 for (std::size_t eq = 0; eq < names.size(); ++eq) {
525 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
526 this->addDependentField(applyBC_[eq]);
527 }
528 }
529
530 // this is what this evaluator provides
531 this->addEvaluatedField(*scatterHolder_);
532
533 if (p.isType<std::string>("Global Data Key"))
534 globalDataKey_ = p.get<std::string>("Global Data Key");
535
536 this->setName(scatterName+" Scatter Residual (Jacobian)");
537}
538
539// **********************************************************************
540template<typename TRAITS,typename LO,typename GO,typename NodeT>
542postRegistrationSetup(typename TRAITS::SetupData /* d */,
544{
545 fieldIds_.resize(scatterFields_.size());
546 // load required field numbers for fast use
547 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
548 // get field ID from DOF manager
549 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
550 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
551 }
552
553 // get the number of nodes (Should be renamed basis)
554 num_nodes = scatterFields_[0].extent(1);
555 num_eq = scatterFields_.size();
556}
557
558// **********************************************************************
559template<typename TRAITS,typename LO,typename GO,typename NodeT>
561preEvaluate(typename TRAITS::PreEvalData d)
562{
563 // extract linear object container
564 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
565
566 if(tpetraContainer_==Teuchos::null) {
567 // extract linear object container
568 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
569 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
570
571 dirichletCounter_ = Teuchos::null;
572 }
573 else {
574 // extract dirichlet counter from container
575 Teuchos::RCP<LOC> tpetraContainer
576 = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
577
578 dirichletCounter_ = tpetraContainer->get_f();
579 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
580 }
581}
582
583// **********************************************************************
584template<typename TRAITS,typename LO,typename GO,typename NodeT>
586evaluateFields(typename TRAITS::EvalData workset)
587{
588 std::vector<GO> GIDs;
589
590 // for convenience pull out some objects from workset
591 std::string blockId = this->wda(workset).block_id;
592 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
593
594 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
595 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
596
597 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
598 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
599
600 // NOTE: A reordering of these loops will likely improve performance
601 // The "getGIDFieldOffsets may be expensive. However the
602 // "getElementGIDs" can be cheaper. However the lookup for LIDs
603 // may be more expensive!
604
605 // scatter operation for each cell in workset
606 auto LIDs = globalIndexer_->getLIDs();
607 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
608 Kokkos::deep_copy(LIDs_h, LIDs);
609 // loop over each field to be scattered
610 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
611 int fieldNum = fieldIds_[fieldIndex];
612 auto scatterFields_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
613 Kokkos::deep_copy(scatterFields_h, scatterFields_[fieldIndex].get_static_view());
614 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
615 std::size_t cellLocalId = localCellIds[worksetCellIndex];
616
617 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
618
619 // this call "should" get the right ordering according to the Intrepid2 basis
620 const std::pair<std::vector<int>,std::vector<int> > & indicePair
621 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
622 const std::vector<int> & elmtOffset = indicePair.first;
623 const std::vector<int> & basisIdMap = indicePair.second;
624
625 // loop over basis functions
626 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
627 int offset = elmtOffset[basis];
628 int lid = LIDs_h(cellLocalId, offset);
629 if(lid<0) // not on this processor
630 continue;
631
632 int basisId = basisIdMap[basis];
633
634 if (checkApplyBC_)
635 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
636 continue;
637
638 // zero out matrix row
639 {
640 std::size_t sz = Jac->getNumEntriesInLocalRow(lid);
641 std::size_t numEntries = 0;
642 typename LOC::CrsMatrixType::nonconst_local_inds_host_view_type rowIndices("indices", sz);
643 typename LOC::CrsMatrixType::nonconst_values_host_view_type rowValues("values", sz);
644
645 // Jac->getLocalRowView(lid,numEntries,rowValues,rowIndices);
646 Jac->getLocalRowCopy(lid,rowIndices,rowValues,numEntries);
647
648 for(std::size_t i=0;i<numEntries;i++)
649 rowValues(i) = 0.0;
650
651 Jac->replaceLocalValues(lid,rowIndices,rowValues);
652 }
653
654 GO gid = GIDs[offset];
655 const ScalarT scatterField = scatterFields_h(worksetCellIndex,basisId);
656
657 r_array[lid] = scatterField.val();
658 dc_array[lid] = 1.0; // mark row as dirichlet
659
660 // loop over the sensitivity indices: all DOFs on a cell
661 std::vector<double> jacRow(scatterField.size(),0.0);
662
663 for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
664 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
665 TEUCHOS_ASSERT(jacRow.size()==GIDs.size());
666
667 Jac->replaceGlobalValues(gid, GIDs, jacRow);
668 }
669 }
670 }
671}
672
673// **********************************************************************
674
675#endif
Pushes residual values into the residual vector for a Newton-based solve.