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// Panzer: A partial differential equation assembly
4// engine for strongly coupled complex multiphysics systems
5//
6// Copyright 2011 NTESS and the Panzer contributors.
7// SPDX-License-Identifier: BSD-3-Clause
8// *****************************************************************************
9// @HEADER
10
11#ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP
12#define PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP
13
14#include "Teuchos_RCP.hpp"
15#include "Teuchos_Assert.hpp"
16
17#include "Phalanx_DataLayout.hpp"
18
19#include "Panzer_GlobalIndexer.hpp"
20#include "Panzer_PureBasis.hpp"
21#include "Panzer_TpetraLinearObjContainer.hpp"
22#include "Panzer_LOCPair_GlobalEvaluationData.hpp"
23#include "Panzer_ParameterList_GlobalEvaluationData.hpp"
24#include "Panzer_GlobalEvaluationDataContainer.hpp"
25
26#include "Phalanx_DataLayout_MDALayout.hpp"
27
28#include "Teuchos_FancyOStream.hpp"
29
30// **********************************************************************
31// Specialization: Residual
32// **********************************************************************
33
34
35template<typename TRAITS,typename LO,typename GO,typename NodeT>
36panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Residual, TRAITS,LO,GO,NodeT>::
37ScatterDirichletResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
38 const Teuchos::ParameterList& p)
39 : globalIndexer_(indexer)
40 , globalDataKey_("Residual Scatter Container")
41{
42 std::string scatterName = p.get<std::string>("Scatter Name");
43 scatterHolder_ =
44 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
45
46 // get names to be evaluated
47 const std::vector<std::string>& names =
48 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
49
50 // grab map from evaluated names to field names
51 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
52
53 // determine if we are scattering an initial condition
54 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
55
56 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
57 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
58 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional ;
59 if (!scatterIC_) {
60 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
61 local_side_id_ = p.get<int>("Local Side ID");
62 scratch_basisIds_.resize(names.size());
63 }
64
65 // build the vector of fields that this is dependent on
66 scatterFields_.resize(names.size());
67 scratch_offsets_.resize(names.size());
68 for (std::size_t eq = 0; eq < names.size(); ++eq) {
69 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
70
71 // tell the field manager that we depend on this field
72 this->addDependentField(scatterFields_[eq]);
73 }
74
75 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
76 if (checkApplyBC_) {
77 applyBC_.resize(names.size());
78 for (std::size_t eq = 0; eq < names.size(); ++eq) {
79 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
80 this->addDependentField(applyBC_[eq]);
81 }
82 }
83
84 // this is what this evaluator provides
85 this->addEvaluatedField(*scatterHolder_);
86
87 if (p.isType<std::string>("Global Data Key"))
88 globalDataKey_ = p.get<std::string>("Global Data Key");
89
90 this->setName(scatterName+" Scatter Residual");
91}
92
93// **********************************************************************
94template<typename TRAITS,typename LO,typename GO,typename NodeT>
95void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Residual, TRAITS,LO,GO,NodeT>::
96postRegistrationSetup(typename TRAITS::SetupData d,
97 PHX::FieldManager<TRAITS>& /* fm */)
98{
99 fieldIds_.resize(scatterFields_.size());
100 const Workset & workset_0 = (*d.worksets_)[0];
101 std::string blockId = this->wda(workset_0).block_id;
102
103 // load required field numbers for fast use
104 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
105 // get field ID from DOF manager
106 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
107 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
108
109 if (!scatterIC_) {
110 const std::pair<std::vector<int>,std::vector<int> > & indicePair
111 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldIds_[fd], side_subcell_dim_, local_side_id_);
112 const std::vector<int> & offsets = indicePair.first;
113 const std::vector<int> & basisIdMap = indicePair.second;
114
115 scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
116 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
117
118 scratch_basisIds_[fd] = PHX::View<int*>("basisIds",basisIdMap.size());
119 Kokkos::deep_copy(scratch_basisIds_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(basisIdMap.data(), basisIdMap.size()));
120
121 } else {
122 const std::vector<int> & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldIds_[fd]);
123 scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
124 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
125 }
126 }
127
128 scratch_lids_ = PHX::View<LO**>("lids",scatterFields_[0].extent(0),
129 globalIndexer_->getElementBlockGIDCount(blockId));
130}
131
132// **********************************************************************
133template<typename TRAITS,typename LO,typename GO,typename NodeT>
134void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Residual, TRAITS,LO,GO,NodeT>::
135preEvaluate(typename TRAITS::PreEvalData d)
136{
137 // extract linear object container
138 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
139
140 if(tpetraContainer_==Teuchos::null) {
141 // extract linear object container
142 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
143 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
144
145 dirichletCounter_ = Teuchos::null;
146 }
147 else {
148 // extract dirichlet counter from container
149 Teuchos::RCP<LOC> tpetraContainer
150 = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
151
152 dirichletCounter_ = tpetraContainer->get_f();
153 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
154 }
155}
156
157// **********************************************************************
158namespace panzer {
159namespace {
160
161template <typename ScalarT,typename LO,typename GO,typename NodeT>
162class ScatterDirichletResidual_Residual_Functor {
163public:
164 typedef typename PHX::Device execution_space;
165 typedef PHX::MDField<const ScalarT,Cell,NODE> ScalarFieldType;
166 typedef PHX::MDField<const bool,Cell,NODE> BoolFieldType;
167
168 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
169 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> dirichlet_counter;
170
171 PHX::View<const LO**> lids; // local indices for unknowns
172 PHX::View<const int*> offsets; // how to get a particular field
173 PHX::View<const int*> basisIds;
174 ScalarFieldType field;
175 BoolFieldType applyBC;
176
177 bool checkApplyBC;
178
179 KOKKOS_INLINE_FUNCTION
180 void operator()(const unsigned int cell) const
181 {
182
183 // loop over the basis functions (currently they are nodes)
184 for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
185 int offset = offsets(basis);
186 LO lid = lids(cell,offset);
187 if (lid<0) continue; // not on this processor
188
189 int basisId = basisIds(basis);
190 if (checkApplyBC)
191 if(!applyBC(cell,basisId)) continue;
192
193 r_data(lid,0) = field(cell,basisId);
194
195 // record that you set a dirichlet condition
196 dirichlet_counter(lid,0) = 1.0;
197
198 } // end basis
199 }
200};
201
202template <typename ScalarT,typename LO,typename GO,typename NodeT>
203class ScatterDirichletResidualIC_Residual_Functor {
204public:
205 typedef typename PHX::Device execution_space;
206 typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
207
208 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
209 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> dirichlet_counter;
210
211 PHX::View<const LO**> lids; // local indices for unknowns
212 PHX::View<const int*> offsets; // how to get a particular field
213 FieldType field;
214
215 KOKKOS_INLINE_FUNCTION
216 void operator()(const unsigned int cell) const
217 {
218
219 // loop over the basis functions (currently they are nodes)
220 for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
221 int offset = offsets(basis);
222 LO lid = lids(cell,offset);
223 if (lid<0) continue; // not on this processor
224
225 r_data(lid,0) = field(cell,basis);
226
227 // record that you set a dirichlet condition
228 dirichlet_counter(lid,0) = 1.0;
229
230 } // end basis
231 }
232};
233}
234}
235
236// **********************************************************************
237template<typename TRAITS,typename LO,typename GO,typename NodeT>
238void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Residual, TRAITS,LO,GO,NodeT>::
239evaluateFields(typename TRAITS::EvalData workset)
240{
241 std::vector<GO> GIDs;
242 std::vector<LO> LIDs;
243
244 // for convenience pull out some objects from workset
245 std::string blockId = this->wda(workset).block_id;
246
247 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
248
249 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
250 tpetraContainer_->get_f() :
251 tpetraContainer_->get_x();
252
253 if (scatterIC_) {
254 ScatterDirichletResidualIC_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
255 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
256 functor.lids = scratch_lids_;
257 functor.dirichlet_counter = dirichletCounter_->getLocalViewDevice(Tpetra::Access::ReadWrite);
258
259 // for each field, do a parallel for loop
260 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
261 functor.offsets = scratch_offsets_[fieldIndex];
262 functor.field = scatterFields_[fieldIndex];
263
264 Kokkos::parallel_for(workset.num_cells,functor);
265 }
266 } else {
267 ScatterDirichletResidual_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
268 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
269 functor.lids = scratch_lids_;
270 functor.dirichlet_counter = dirichletCounter_->getLocalViewDevice(Tpetra::Access::ReadWrite);
271
272 // for each field, do a parallel for loop
273 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
274 functor.offsets = scratch_offsets_[fieldIndex];
275 functor.field = scatterFields_[fieldIndex];
276 if (checkApplyBC_) functor.applyBC = applyBC_[fieldIndex];
277 functor.checkApplyBC = checkApplyBC_;
278 functor.basisIds = scratch_basisIds_[fieldIndex];
279
280 Kokkos::parallel_for(workset.num_cells,functor);
281 }
282 }
283
284}
285
286// **********************************************************************
287// Specialization: Tangent
288// **********************************************************************
289
290
291template<typename TRAITS,typename LO,typename GO,typename NodeT>
292panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Tangent, TRAITS,LO,GO,NodeT>::
293ScatterDirichletResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
294 const Teuchos::ParameterList& p)
295 : globalIndexer_(indexer)
296 , globalDataKey_("Residual Scatter Container")
297{
298 std::string scatterName = p.get<std::string>("Scatter Name");
299 scatterHolder_ =
300 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
301
302 // get names to be evaluated
303 const std::vector<std::string>& names =
304 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
305
306 // grab map from evaluated names to field names
307 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
308
309 // determine if we are scattering an initial condition
310 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
311
312 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
313 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
314 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional ;
315 if (!scatterIC_) {
316 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
317 local_side_id_ = p.get<int>("Local Side ID");
318 scratch_basisIds_.resize(names.size());
319 }
320
321 // build the vector of fields that this is dependent on
322 scatterFields_.resize(names.size());
323 scratch_offsets_.resize(names.size());
324 for (std::size_t eq = 0; eq < names.size(); ++eq) {
325 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
326
327 // tell the field manager that we depend on this field
328 this->addDependentField(scatterFields_[eq]);
329 }
330
331 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
332 if (checkApplyBC_) {
333 applyBC_.resize(names.size());
334 for (std::size_t eq = 0; eq < names.size(); ++eq) {
335 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
336 this->addDependentField(applyBC_[eq]);
337 }
338 }
339
340 // this is what this evaluator provides
341 this->addEvaluatedField(*scatterHolder_);
342
343 if (p.isType<std::string>("Global Data Key"))
344 globalDataKey_ = p.get<std::string>("Global Data Key");
345
346 this->setName(scatterName+" Scatter Tangent");
347}
348
349// **********************************************************************
350template<typename TRAITS,typename LO,typename GO,typename NodeT>
351void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Tangent, TRAITS,LO,GO,NodeT>::
352postRegistrationSetup(typename TRAITS::SetupData d,
353 PHX::FieldManager<TRAITS>& /* fm */)
354{
355 fieldIds_.resize(scatterFields_.size());
356 const Workset & workset_0 = (*d.worksets_)[0];
357 std::string blockId = this->wda(workset_0).block_id;
358
359 // load required field numbers for fast use
360 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
361 // get field ID from DOF manager
362 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
363 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
364
365 if (!scatterIC_) {
366 const std::pair<std::vector<int>,std::vector<int> > & indicePair
367 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldIds_[fd], side_subcell_dim_, local_side_id_);
368 const std::vector<int> & offsets = indicePair.first;
369 const std::vector<int> & basisIdMap = indicePair.second;
370
371 scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
372 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
373
374 scratch_basisIds_[fd] = PHX::View<int*>("basisIds",basisIdMap.size());
375 Kokkos::deep_copy(scratch_basisIds_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(basisIdMap.data(), basisIdMap.size()));
376
377 } else {
378 const std::vector<int> & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldIds_[fd]);
379 scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
380 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
381 }
382 }
383
384 scratch_lids_ = PHX::View<LO**>("lids",scatterFields_[0].extent(0),
385 globalIndexer_->getElementBlockGIDCount(blockId));
386
387}
388
389// **********************************************************************
390template<typename TRAITS,typename LO,typename GO,typename NodeT>
391void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Tangent, TRAITS,LO,GO,NodeT>::
392preEvaluate(typename TRAITS::PreEvalData d)
393{
394 // extract linear object container
395 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
396
397 if(tpetraContainer_==Teuchos::null) {
398 // extract linear object container
399 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
400 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
401
402 dirichletCounter_ = Teuchos::null;
403 }
404 else {
405 // extract dirichlet counter from container
406 Teuchos::RCP<LOC> tpetraContainer
407 = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
408
409 dirichletCounter_ = tpetraContainer->get_f();
410 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
411 }
412
413 using Teuchos::RCP;
414 using Teuchos::rcp_dynamic_cast;
415
416 // this is the list of parameters and their names that this scatter has to account for
417 std::vector<std::string> activeParameters =
418 rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
419
420 dfdpFieldsVoV_.initialize("ScatterResidual_Tpetra<Tangent>::dfdpFieldsVoV_",activeParameters.size());
421
422 for(std::size_t i=0;i<activeParameters.size();i++) {
423 RCP<typename LOC::VectorType> vec =
424 rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
425 auto dfdp_view = vec->getLocalViewDevice(Tpetra::Access::ReadWrite);
426
427 dfdpFieldsVoV_.addView(dfdp_view,i);
428 }
429
430 dfdpFieldsVoV_.syncHostToDevice();
431
432}
433
434// **********************************************************************
435namespace panzer {
436namespace {
437
438template <typename ScalarT,typename LO,typename GO,typename NodeT>
439class ScatterDirichletResidual_Tangent_Functor {
440public:
441 typedef typename PHX::Device execution_space;
442 typedef PHX::MDField<const ScalarT,Cell,NODE> ScalarFieldType;
443 typedef PHX::MDField<const bool,Cell,NODE> BoolFieldType;
444
445 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
446 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> dirichlet_counter;
447
448 Kokkos::View<Kokkos::View<double**,Kokkos::LayoutLeft,PHX::Device>*> dfdp_fields; // tangent fields
449 double num_params;
450
451 PHX::View<const LO**> lids; // local indices for unknowns
452 PHX::View<const int*> offsets; // how to get a particular field
453 PHX::View<const int*> basisIds;
454 ScalarFieldType field;
455 BoolFieldType applyBC;
456
457 bool checkApplyBC;
458
459 KOKKOS_INLINE_FUNCTION
460 void operator()(const unsigned int cell) const
461 {
462
463 // loop over the basis functions (currently they are nodes)
464 for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
465 int offset = offsets(basis);
466 LO lid = lids(cell,offset);
467 if (lid<0) continue; // not on this processor
468
469 int basisId = basisIds(basis);
470 if (checkApplyBC)
471 if(!applyBC(cell,basisId)) continue;
472
473 r_data(lid,0) = field(cell,basisId).val();
474
475 // loop over the tangents
476 for(int i_param=0; i_param<num_params; i_param++)
477 dfdp_fields(i_param)(lid,0) = field(cell,basisId).fastAccessDx(i_param);
478
479 // record that you set a dirichlet condition
480 dirichlet_counter(lid,0) = 1.0;
481
482 } // end basis
483 }
484};
485
486template <typename ScalarT,typename LO,typename GO,typename NodeT>
487class ScatterDirichletResidualIC_Tangent_Functor {
488public:
489 typedef typename PHX::Device execution_space;
490 typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
491
492 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
493 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> dirichlet_counter;
494
495 Kokkos::View<Kokkos::View<double**,Kokkos::LayoutLeft,PHX::Device>*> dfdp_fields; // tangent fields
496 double num_params;
497
498 PHX::View<const LO**> lids; // local indices for unknowns
499 PHX::View<const int*> offsets; // how to get a particular field
500 FieldType field;
501
502 KOKKOS_INLINE_FUNCTION
503 void operator()(const unsigned int cell) const
504 {
505
506 // loop over the basis functions (currently they are nodes)
507 for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
508 int offset = offsets(basis);
509 LO lid = lids(cell,offset);
510 if (lid<0) continue; // not on this processor
511
512 r_data(lid,0) = field(cell,basis).val();
513
514 // loop over the tangents
515 for(int i_param=0; i_param<num_params; i_param++)
516 dfdp_fields(i_param)(lid,0) = field(cell,basis).fastAccessDx(i_param);
517
518 // record that you set a dirichlet condition
519 dirichlet_counter(lid,0) = 1.0;
520
521 } // end basis
522 }
523};
524}
525}
526
527// **********************************************************************
528template<typename TRAITS,typename LO,typename GO,typename NodeT>
529void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Tangent, TRAITS,LO,GO,NodeT>::
530evaluateFields(typename TRAITS::EvalData workset)
531{
532 std::vector<GO> GIDs;
533 std::vector<LO> LIDs;
534
535 // for convenience pull out some objects from workset
536 std::string blockId = this->wda(workset).block_id;
537
538 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
539
540 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
541 tpetraContainer_->get_f() :
542 tpetraContainer_->get_x();
543
544 if (scatterIC_) {
545 ScatterDirichletResidualIC_Tangent_Functor<ScalarT,LO,GO,NodeT> functor;
546 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
547 functor.lids = scratch_lids_;
548 functor.dirichlet_counter = dirichletCounter_->getLocalViewDevice(Tpetra::Access::ReadWrite);
549 functor.dfdp_fields = dfdpFieldsVoV_.getViewDevice();
550
551 // for each field, do a parallel for loop
552 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
553 functor.offsets = scratch_offsets_[fieldIndex];
554 functor.field = scatterFields_[fieldIndex];
555 functor.num_params = Kokkos::dimension_scalar(scatterFields_[fieldIndex].get_view())-1;
556
557 Kokkos::parallel_for(workset.num_cells,functor);
558 }
559 } else {
560 ScatterDirichletResidual_Tangent_Functor<ScalarT,LO,GO,NodeT> functor;
561 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
562 functor.lids = scratch_lids_;
563 functor.dirichlet_counter = dirichletCounter_->getLocalViewDevice(Tpetra::Access::ReadWrite);
564 functor.dfdp_fields = dfdpFieldsVoV_.getViewDevice();
565
566 // for each field, do a parallel for loop
567 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
568 functor.offsets = scratch_offsets_[fieldIndex];
569 functor.field = scatterFields_[fieldIndex];
570 if (checkApplyBC_) functor.applyBC = applyBC_[fieldIndex];
571 functor.checkApplyBC = checkApplyBC_;
572 functor.basisIds = scratch_basisIds_[fieldIndex];
573 functor.num_params = Kokkos::dimension_scalar(scatterFields_[fieldIndex].get_view())-1;
574
575 Kokkos::parallel_for(workset.num_cells,functor);
576 }
577 }
578
579}
580
581// **********************************************************************
582// Specialization: Jacobian
583// **********************************************************************
584
585template<typename TRAITS,typename LO,typename GO,typename NodeT>
586panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Jacobian, TRAITS,LO,GO,NodeT>::
587ScatterDirichletResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
588 const Teuchos::ParameterList& p)
589 : globalIndexer_(indexer)
590 , globalDataKey_("Residual Scatter Container")
591{
592 std::string scatterName = p.get<std::string>("Scatter Name");
593 scatterHolder_ =
594 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
595
596 // get names to be evaluated
597 const std::vector<std::string>& names =
598 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
599
600 // grab map from evaluated names to field names
601 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
602
603 Teuchos::RCP<PHX::DataLayout> dl =
604 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
605
606 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
607 local_side_id_ = p.get<int>("Local Side ID");
608
609 // build the vector of fields that this is dependent on
610 scatterFields_.resize(names.size());
611 for (std::size_t eq = 0; eq < names.size(); ++eq) {
612 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
613
614 // tell the field manager that we depend on this field
615 this->addDependentField(scatterFields_[eq]);
616 }
617
618 checkApplyBC_ = p.get<bool>("Check Apply BC");
619 if (checkApplyBC_) {
620 applyBC_.resize(names.size());
621 for (std::size_t eq = 0; eq < names.size(); ++eq) {
622 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
623 this->addDependentField(applyBC_[eq]);
624 }
625 }
626
627 // this is what this evaluator provides
628 this->addEvaluatedField(*scatterHolder_);
629
630 if (p.isType<std::string>("Global Data Key"))
631 globalDataKey_ = p.get<std::string>("Global Data Key");
632
633 this->setName(scatterName+" Scatter Residual (Jacobian)");
634}
635
636// **********************************************************************
637template<typename TRAITS,typename LO,typename GO,typename NodeT>
638void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Jacobian, TRAITS,LO,GO,NodeT>::
639postRegistrationSetup(typename TRAITS::SetupData /* d */,
640 PHX::FieldManager<TRAITS>& /* fm */)
641{
642 fieldIds_.resize(scatterFields_.size());
643 // load required field numbers for fast use
644 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
645 // get field ID from DOF manager
646 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
647 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
648 }
649
650 // get the number of nodes (Should be renamed basis)
651 num_nodes = scatterFields_[0].extent(1);
652 num_eq = scatterFields_.size();
653}
654
655// **********************************************************************
656template<typename TRAITS,typename LO,typename GO,typename NodeT>
657void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Jacobian, TRAITS,LO,GO,NodeT>::
658preEvaluate(typename TRAITS::PreEvalData d)
659{
660 // extract linear object container
661 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
662
663 if(tpetraContainer_==Teuchos::null) {
664 // extract linear object container
665 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
666 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
667
668 dirichletCounter_ = Teuchos::null;
669 }
670 else {
671 // extract dirichlet counter from container
672 Teuchos::RCP<LOC> tpetraContainer
673 = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
674
675 dirichletCounter_ = tpetraContainer->get_f();
676 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
677 }
678}
679
680// **********************************************************************
681template<typename TRAITS,typename LO,typename GO,typename NodeT>
682void panzer::ScatterDirichletResidual_Tpetra<panzer::Traits::Jacobian, TRAITS,LO,GO,NodeT>::
683evaluateFields(typename TRAITS::EvalData workset)
684{
685 std::vector<GO> GIDs;
686
687 // for convenience pull out some objects from workset
688 std::string blockId = this->wda(workset).block_id;
689 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
690
691 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
692 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
693
694 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
695 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
696
697 // NOTE: A reordering of these loops will likely improve performance
698 // The "getGIDFieldOffsets may be expensive. However the
699 // "getElementGIDs" can be cheaper. However the lookup for LIDs
700 // may be more expensive!
701
702 // scatter operation for each cell in workset
703 auto LIDs = globalIndexer_->getLIDs();
704 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
705 Kokkos::deep_copy(LIDs_h, LIDs);
706 // loop over each field to be scattered
707 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
708 int fieldNum = fieldIds_[fieldIndex];
709 auto scatterFields_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
710 Kokkos::deep_copy(scatterFields_h, scatterFields_[fieldIndex].get_static_view());
711 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
712 std::size_t cellLocalId = localCellIds[worksetCellIndex];
713
714 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
715
716 // this call "should" get the right ordering according to the Intrepid2 basis
717 const std::pair<std::vector<int>,std::vector<int> > & indicePair
718 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
719 const std::vector<int> & elmtOffset = indicePair.first;
720 const std::vector<int> & basisIdMap = indicePair.second;
721
722 // loop over basis functions
723 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
724 int offset = elmtOffset[basis];
725 int lid = LIDs_h(cellLocalId, offset);
726 if(lid<0) // not on this processor
727 continue;
728
729 int basisId = basisIdMap[basis];
730
731 if (checkApplyBC_)
732 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
733 continue;
734
735 // zero out matrix row
736 {
737 std::size_t sz = Jac->getNumEntriesInLocalRow(lid);
738 std::size_t numEntries = 0;
739 typename LOC::CrsMatrixType::nonconst_local_inds_host_view_type rowIndices("indices", sz);
740 typename LOC::CrsMatrixType::nonconst_values_host_view_type rowValues("values", sz);
741
742 // Jac->getLocalRowView(lid,numEntries,rowValues,rowIndices);
743 Jac->getLocalRowCopy(lid,rowIndices,rowValues,numEntries);
744
745 for(std::size_t i=0;i<numEntries;i++)
746 rowValues(i) = 0.0;
747
748 Jac->replaceLocalValues(lid,rowIndices,rowValues);
749 }
750
751 GO gid = GIDs[offset];
752 const ScalarT scatterField = scatterFields_h(worksetCellIndex,basisId);
753
754 r_array[lid] = scatterField.val();
755 dc_array[lid] = 1.0; // mark row as dirichlet
756
757 // loop over the sensitivity indices: all DOFs on a cell
758 std::vector<double> jacRow(scatterField.size(),0.0);
759
760 for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
761 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
762 TEUCHOS_ASSERT(jacRow.size()==GIDs.size());
763
764 Jac->replaceGlobalValues(gid, GIDs, jacRow);
765 }
766 }
767 }
768}
769
770// **********************************************************************
771
772#endif
offsets
PHX::View< const int * > offsets
Definition Panzer_DOFCurl_impl.hpp:148
Panzer_GlobalEvaluationDataContainer.hpp
Panzer_GlobalIndexer.hpp
field
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.
Definition Panzer_Integrator_CurlBasisDotVector_impl.hpp:559
Panzer_LOCPair_GlobalEvaluationData.hpp
Panzer_ParameterList_GlobalEvaluationData.hpp
Panzer_PureBasis.hpp
checkApplyBC
bool checkApplyBC
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:177
num_params
double num_params
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:449
dfdp_fields
Kokkos::View< Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > * > dfdp_fields
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:448
dirichlet_counter
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > dirichlet_counter
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:169
lids
PHX::View< const LO ** > lids
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:171
applyBC
BoolFieldType applyBC
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:175
basisIds
PHX::View< const int * > basisIds
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:173
r_data
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > r_data
Definition Panzer_ScatterDirichletResidual_Tpetra_impl.hpp:168
Panzer_TpetraLinearObjContainer.hpp
panzer::ScatterDirichletResidual_Tpetra< panzer::Traits::Jacobian, TRAITS, LO, GO, NodeT >::ScalarT
panzer::Traits::Jacobian::ScalarT ScalarT
Definition Panzer_ScatterDirichletResidual_Tpetra_decl.hpp:228
panzer::ScatterDirichletResidual_Tpetra
Pushes residual values into the residual vector for a Newton-based solve.
Definition Panzer_ScatterDirichletResidual_Tpetra_decl.hpp:44
panzer::WorksetDetails::block_id
std::string block_id
DEPRECATED - use: getElementBlock()
Definition Panzer_Workset.hpp:98
panzer::FieldType
FieldType
The type of discretization to use for a field pattern.
Definition Panzer_FieldType.hpp:18
Generated by doxygen 1.9.8