Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterDirichletResidual_BlockedEpetra_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
43#ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_IMPL_HPP
44#define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_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
58#include "Panzer_PureBasis.hpp"
61
62#include "Phalanx_DataLayout_MDALayout.hpp"
63
64#include "Thyra_SpmdVectorBase.hpp"
65#include "Thyra_ProductVectorBase.hpp"
66#include "Thyra_DefaultProductVector.hpp"
67#include "Thyra_BlockedLinearOpBase.hpp"
68#include "Thyra_get_Epetra_Operator.hpp"
69
70#include "Teuchos_FancyOStream.hpp"
71
72#include <unordered_map>
73
74// **********************************************************************
75// Specialization: Residual
76// **********************************************************************
77
78
79template<typename TRAITS,typename LO,typename GO>
81ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
82 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
83 const Teuchos::ParameterList& p,
84 bool /* useDiscreteAdjoint */)
85 : rowIndexers_(rIndexers)
86 , colIndexers_(cIndexers)
87 , globalDataKey_("Residual Scatter Container")
88{
89 std::string scatterName = p.get<std::string>("Scatter Name");
90 scatterHolder_ =
91 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
92
93 // get names to be evaluated
94 const std::vector<std::string>& names =
95 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
96
97 // grab map from evaluated names to field names
98 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
99
100 // determine if we are scattering an initial condition
101 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
102
103 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
104 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
105 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
106 if (!scatterIC_) {
107 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
108 local_side_id_ = p.get<int>("Local Side ID");
109 }
110
111 // build the vector of fields that this is dependent on
112 scatterFields_.resize(names.size());
113 for (std::size_t eq = 0; eq < names.size(); ++eq) {
114 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
115
116 // tell the field manager that we depend on this field
117 this->addDependentField(scatterFields_[eq]);
118 }
119
120 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
121 if (checkApplyBC_) {
122 applyBC_.resize(names.size());
123 for (std::size_t eq = 0; eq < names.size(); ++eq) {
124 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
125 this->addDependentField(applyBC_[eq]);
126 }
127 }
128
129 // this is what this evaluator provides
130 this->addEvaluatedField(*scatterHolder_);
131
132 if (p.isType<std::string>("Global Data Key"))
133 globalDataKey_ = p.get<std::string>("Global Data Key");
134
135 this->setName(scatterName+" Scatter Residual");
136}
137
138// **********************************************************************
139template<typename TRAITS,typename LO,typename GO>
141postRegistrationSetup(typename TRAITS::SetupData /* d */,
143{
144 indexerIds_.resize(scatterFields_.size());
145 subFieldIds_.resize(scatterFields_.size());
146
147 // load required field numbers for fast use
148 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
149 // get field ID from DOF manager
150 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
151
152 indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
153 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
154 }
155
156 // get the number of nodes (Should be renamed basis)
157 num_nodes = scatterFields_[0].extent(1);
158}
159
160// **********************************************************************
161template<typename TRAITS,typename LO,typename GO>
163preEvaluate(typename TRAITS::PreEvalData d)
164{
167
168 using Teuchos::rcp_dynamic_cast;
170
171 // extract dirichlet counter from container
172 Teuchos::RCP<BLOC> blockContainer
173 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
174
175 dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
176 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
177
178 // extract linear object container
179 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
180 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
181
182 // if its blocked do this
183 if(blockedContainer!=Teuchos::null)
184 r_ = (!scatterIC_) ?
185 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
186 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
187 else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
188 r_ = (!scatterIC_) ?
189 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
190 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
191
192 TEUCHOS_ASSERT(r_!=Teuchos::null);
193}
194
195// **********************************************************************
196template<typename TRAITS,typename LO,typename GO>
198evaluateFields(typename TRAITS::EvalData workset)
199{
200 using Teuchos::RCP;
201 using Teuchos::ArrayRCP;
202 using Teuchos::ptrFromRef;
203 using Teuchos::rcp_dynamic_cast;
204
205 using Thyra::VectorBase;
206 using Thyra::SpmdVectorBase;
208
209 // for convenience pull out some objects from workset
210 std::string blockId = this->wda(workset).block_id;
211 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
212
213 // NOTE: A reordering of these loops will likely improve performance
214 // The "getGIDFieldOffsets may be expensive. However the
215 // "getElementGIDs" can be cheaper. However the lookup for LIDs
216 // may be more expensive!
217
218 // loop over each field to be scattered
219 Teuchos::ArrayRCP<double> local_r;
220 Teuchos::ArrayRCP<double> local_dc;
221 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
222 int rowIndexer = indexerIds_[fieldIndex];
223 int subFieldNum = subFieldIds_[fieldIndex];
224
225 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
226 ->getNonconstLocalData(ptrFromRef(local_dc));
227
228 // grab local data for inputing
229 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
230 ->getNonconstLocalData(ptrFromRef(local_r));
231
232 auto subRowIndexer = rowIndexers_[rowIndexer];
233 auto LIDs = subRowIndexer->getLIDs();
234 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
235 Kokkos::deep_copy(LIDs_h, LIDs);
236
237 auto field = scatterFields_[fieldIndex].get_view();
238 auto field_h = Kokkos::create_mirror_view(field);
239 Kokkos::deep_copy(field_h, field);
240
241 BCFieldType::array_type::HostMirror applyBC_h;
242 if(checkApplyBC_){
243 auto applyBC = applyBC_[fieldIndex].get_static_view();
244 applyBC_h = Kokkos::create_mirror_view(applyBC);
245 Kokkos::deep_copy(applyBC_h, applyBC);
246 }
247
248 // scatter operation for each cell in workset
249 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
250 std::size_t cellLocalId = localCellIds[worksetCellIndex];
251
252 if (!scatterIC_) {
253 // this call "should" get the right ordering according to the Intrepid2 basis
254 const std::pair<std::vector<int>,std::vector<int> > & indicePair
255 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
256 const std::vector<int> & elmtOffset = indicePair.first;
257 const std::vector<int> & basisIdMap = indicePair.second;
258
259 // loop over basis functions
260 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
261 int offset = elmtOffset[basis];
262 int lid = LIDs_h(cellLocalId, offset);
263 if(lid<0) // not on this processor!
264 continue;
265
266 int basisId = basisIdMap[basis];
267
268 if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
269 continue;
270
271 local_r[lid] = field_h(worksetCellIndex,basisId);
272
273 // record that you set a dirichlet condition
274 local_dc[lid] = 1.0;
275 }
276 } else {
277 // this call "should" get the right ordering according to the Intrepid2 basis
278 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
279
280 // loop over basis functions
281 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
282 int offset = elmtOffset[basis];
283 int lid = LIDs_h(cellLocalId, offset);
284 if(lid<0) // not on this processor!
285 continue;
286
287 local_r[lid] = field_h(worksetCellIndex,basis);
288
289 // record that you set a dirichlet condition
290 local_dc[lid] = 1.0;
291 }
292 }
293 }
294 }
295}
296
297// **********************************************************************
298// Specialization: Tangent
299// **********************************************************************
300
301
302template<typename TRAITS,typename LO,typename GO>
304ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
305 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
306 const Teuchos::ParameterList& p,
307 bool /* useDiscreteAdjoint */)
308 : rowIndexers_(rIndexers)
309 , colIndexers_(cIndexers)
310 , globalDataKey_("Residual Scatter Container")
311{
312 std::string scatterName = p.get<std::string>("Scatter Name");
313 scatterHolder_ =
314 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
315
316 // get names to be evaluated
317 const std::vector<std::string>& names =
318 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
319
320 // grab map from evaluated names to field names
321 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
322
323 // determine if we are scattering an initial condition
324 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
325
326 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
327 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
328 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
329 if (!scatterIC_) {
330 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
331 local_side_id_ = p.get<int>("Local Side ID");
332 }
333
334 // build the vector of fields that this is dependent on
335 scatterFields_.resize(names.size());
336 for (std::size_t eq = 0; eq < names.size(); ++eq) {
337 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
338
339 // tell the field manager that we depend on this field
340 this->addDependentField(scatterFields_[eq]);
341 }
342
343 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
344 if (checkApplyBC_) {
345 applyBC_.resize(names.size());
346 for (std::size_t eq = 0; eq < names.size(); ++eq) {
347 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
348 this->addDependentField(applyBC_[eq]);
349 }
350 }
351
352 // this is what this evaluator provides
353 this->addEvaluatedField(*scatterHolder_);
354
355 if (p.isType<std::string>("Global Data Key"))
356 globalDataKey_ = p.get<std::string>("Global Data Key");
357
358 this->setName(scatterName+" Scatter Tangent");
359}
360
361// **********************************************************************
362template<typename TRAITS,typename LO,typename GO>
364postRegistrationSetup(typename TRAITS::SetupData /* d */,
366{
367 indexerIds_.resize(scatterFields_.size());
368 subFieldIds_.resize(scatterFields_.size());
369
370 // load required field numbers for fast use
371 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
372 // get field ID from DOF manager
373 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
374
375 indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
376 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
377 }
378
379 // get the number of nodes (Should be renamed basis)
380 num_nodes = scatterFields_[0].extent(1);
381}
382
383// **********************************************************************
384template<typename TRAITS,typename LO,typename GO>
386preEvaluate(typename TRAITS::PreEvalData d)
387{
390
391 using Teuchos::rcp_dynamic_cast;
393
394 // extract dirichlet counter from container
395 Teuchos::RCP<BLOC> blockContainer
396 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
397
398 dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
399 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
400
401 // extract linear object container
402 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
403 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
404
405 // if its blocked do this
406 if(blockedContainer!=Teuchos::null)
407 r_ = (!scatterIC_) ?
408 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
409 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
410 else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
411 r_ = (!scatterIC_) ?
412 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
413 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
414
415 TEUCHOS_ASSERT(r_!=Teuchos::null);
416}
417
418// **********************************************************************
419template<typename TRAITS,typename LO,typename GO>
421evaluateFields(typename TRAITS::EvalData workset)
422{
423 TEUCHOS_ASSERT(false);
424
425 using Teuchos::RCP;
426 using Teuchos::ArrayRCP;
427 using Teuchos::ptrFromRef;
428 using Teuchos::rcp_dynamic_cast;
429
430 using Thyra::VectorBase;
431 using Thyra::SpmdVectorBase;
433
434 // for convenience pull out some objects from workset
435 std::string blockId = this->wda(workset).block_id;
436 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
437
438 // NOTE: A reordering of these loops will likely improve performance
439 // The "getGIDFieldOffsets may be expensive. However the
440 // "getElementGIDs" can be cheaper. However the lookup for LIDs
441 // may be more expensive!
442
443 // loop over each field to be scattered
444 Teuchos::ArrayRCP<double> local_r;
445 Teuchos::ArrayRCP<double> local_dc;
446 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
447 int rowIndexer = indexerIds_[fieldIndex];
448 int subFieldNum = subFieldIds_[fieldIndex];
449
450 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
451 ->getNonconstLocalData(ptrFromRef(local_dc));
452
453 // grab local data for inputing
454 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
455 ->getNonconstLocalData(ptrFromRef(local_r));
456
457 auto subRowIndexer = rowIndexers_[rowIndexer];
458
459 auto LIDs = subRowIndexer->getLIDs();
460 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
461 Kokkos::deep_copy(LIDs_h, LIDs);
462
463 auto field = scatterFields_[fieldIndex].get_view();
464 auto field_h = Kokkos::create_mirror_view(field);
465 Kokkos::deep_copy(field_h, field);
466
467 BCFieldType::array_type::HostMirror applyBC_h;
468 if(checkApplyBC_){
469 auto applyBC = applyBC_[fieldIndex].get_static_view();
470 applyBC_h = Kokkos::create_mirror_view(applyBC);
471 Kokkos::deep_copy(applyBC_h, applyBC);
472 }
473
474 // scatter operation for each cell in workset
475 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
476 std::size_t cellLocalId = localCellIds[worksetCellIndex];
477
478 if (!scatterIC_) {
479 // this call "should" get the right ordering according to the Intrepid2 basis
480 const std::pair<std::vector<int>,std::vector<int> > & indicePair
481 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
482 const std::vector<int> & elmtOffset = indicePair.first;
483 const std::vector<int> & basisIdMap = indicePair.second;
484
485 // loop over basis functions
486 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
487 int offset = elmtOffset[basis];
488 int lid = LIDs_h(cellLocalId, offset);
489 if(lid<0) // not on this processor!
490 continue;
491
492 int basisId = basisIdMap[basis];
493
494 if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
495 continue;
496
497 local_r[lid] = field_h(worksetCellIndex,basisId).val();
498
499 // record that you set a dirichlet condition
500 local_dc[lid] = 1.0;
501 }
502 } else {
503 // this call "should" get the right ordering according to the Intrepid2 basis
504 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
505
506 // loop over basis functions
507 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
508 int offset = elmtOffset[basis];
509 int lid = LIDs_h(cellLocalId, offset);
510 if(lid<0) // not on this processor!
511 continue;
512
513 local_r[lid] = field_h(worksetCellIndex,basis).val();
514
515 // record that you set a dirichlet condition
516 local_dc[lid] = 1.0;
517 }
518 }
519 }
520 }
521}
522
523// **********************************************************************
524// Specialization: Jacobian
525// **********************************************************************
526
527template<typename TRAITS,typename LO,typename GO>
529ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
530 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
531 const Teuchos::ParameterList& p,
532 bool /* useDiscreteAdjoint */)
533 : rowIndexers_(rIndexers)
534 , colIndexers_(cIndexers)
535 , globalDataKey_("Residual Scatter Container")
536{
537 std::string scatterName = p.get<std::string>("Scatter Name");
538 scatterHolder_ =
539 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
540
541 // get names to be evaluated
542 const std::vector<std::string>& names =
543 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
544
545 // grab map from evaluated names to field names
546 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
547
548 Teuchos::RCP<PHX::DataLayout> dl =
549 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
550
551 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
552 local_side_id_ = p.get<int>("Local Side ID");
553
554 // build the vector of fields that this is dependent on
555 scatterFields_.resize(names.size());
556 for (std::size_t eq = 0; eq < names.size(); ++eq) {
557 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
558
559 // tell the field manager that we depend on this field
560 this->addDependentField(scatterFields_[eq]);
561 }
562
563 checkApplyBC_ = p.get<bool>("Check Apply BC");
564 if (checkApplyBC_) {
565 applyBC_.resize(names.size());
566 for (std::size_t eq = 0; eq < names.size(); ++eq) {
567 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
568 this->addDependentField(applyBC_[eq]);
569 }
570 }
571
572 // this is what this evaluator provides
573 this->addEvaluatedField(*scatterHolder_);
574
575 if (p.isType<std::string>("Global Data Key"))
576 globalDataKey_ = p.get<std::string>("Global Data Key");
577
578 if(colIndexers_.size()==0)
579 colIndexers_ = rowIndexers_;
580
581 this->setName(scatterName+" Scatter Residual (Jacobian)");
582}
583
584// **********************************************************************
585template<typename TRAITS,typename LO,typename GO>
587postRegistrationSetup(typename TRAITS::SetupData /* d */,
589{
590 indexerIds_.resize(scatterFields_.size());
591 subFieldIds_.resize(scatterFields_.size());
592
593 // load required field numbers for fast use
594 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
595 // get field ID from DOF manager
596 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
597
598 indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
599 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
600 }
601
602 // get the number of nodes (Should be renamed basis)
603 num_nodes = scatterFields_[0].extent(1);
604 num_eq = scatterFields_.size();
605}
606
607// **********************************************************************
608template<typename TRAITS,typename LO,typename GO>
610preEvaluate(typename TRAITS::PreEvalData d)
611{
613
614 using Teuchos::rcp_dynamic_cast;
615
616 // extract dirichlet counter from container
617 Teuchos::RCP<const BLOC> blockContainer
618 = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
619
620 dirichletCounter_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
621 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
622
623 // extract linear object container
624 blockContainer = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_),true);
625 TEUCHOS_ASSERT(!Teuchos::is_null(blockContainer));
626
627 r_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f());
628 Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockContainer->get_A());
629}
630
631// **********************************************************************
632template<typename TRAITS,typename LO,typename GO>
634evaluateFields(typename TRAITS::EvalData workset)
635{
636 using Teuchos::RCP;
637 using Teuchos::ArrayRCP;
638 using Teuchos::ptrFromRef;
639 using Teuchos::rcp_dynamic_cast;
640
641 using Thyra::SpmdVectorBase;
642
643 // for convenience pull out some objects from workset
644 std::string blockId = this->wda(workset).block_id;
645 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
646
647 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
648
649 std::vector<int> blockOffsets;
650 computeBlockOffsets(blockId,colIndexers_,blockOffsets);
651
652 std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
653
654 // NOTE: A reordering of these loops will likely improve performance
655 // The "getGIDFieldOffsets may be expensive. However the
656 // "getElementGIDs" can be cheaper. However the lookup for LIDs
657 // may be more expensive!
658
659 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
660 int rowIndexer = indexerIds_[fieldIndex];
661 int subFieldNum = subFieldIds_[fieldIndex];
662
663 // loop over each field to be scattered
664 Teuchos::ArrayRCP<double> local_dc;
665 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
666 ->getNonconstLocalData(ptrFromRef(local_dc));
667
668 // grab local data for inputing
669 Teuchos::ArrayRCP<double> local_r;
670 if(r_!=Teuchos::null)
671 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
672 ->getNonconstLocalData(ptrFromRef(local_r));
673
674 auto subRowIndexer = rowIndexers_[rowIndexer];
675 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
676 const std::vector<int> & subElmtOffset = subIndicePair.first;
677 const std::vector<int> & subBasisIdMap = subIndicePair.second;
678
679 auto rLIDs = subRowIndexer->getLIDs();
680 auto rLIDs_h = Kokkos::create_mirror_view(rLIDs);
681 Kokkos::deep_copy(rLIDs_h, rLIDs);
682
683 auto field = scatterFields_[fieldIndex].get_view();
684 auto field_h = Kokkos::create_mirror_view(field);
685 Kokkos::deep_copy(field_h, field);
686
687 BCFieldType::array_type::HostMirror applyBC_h;
688 if(checkApplyBC_){
689 auto applyBC = applyBC_[fieldIndex].get_static_view();
690 applyBC_h = Kokkos::create_mirror_view(applyBC);
691 Kokkos::deep_copy(applyBC_h, applyBC);
692 }
693
694 // scatter operation for each cell in workset
695 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
696 std::size_t cellLocalId = localCellIds[worksetCellIndex];
697
698 // loop over basis functions
699 for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
700 int offset = subElmtOffset[basis];
701 int lid = rLIDs_h(cellLocalId, offset);
702 if(lid<0) // not on this processor
703 continue;
704
705 int basisId = subBasisIdMap[basis];
706
707 if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
708 continue;
709
710 // zero out matrix row
711 for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
712 int start = blockOffsets[colIndexer];
713 int end = blockOffsets[colIndexer+1];
714
715 if(end-start<=0)
716 continue;
717
718 // check hash table for jacobian sub block
719 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
720 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
721 // if you didn't find one before, add it to the hash table
722 if(subJac==Teuchos::null) {
723 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
724
725 // block operator is null, don't do anything (it is excluded)
726 if(Teuchos::is_null(tOp))
727 continue;
728
729 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
730 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
731 jacEpetraBlocks[blockIndex] = subJac;
732 }
733
734 int numEntries = 0;
735 int * rowIndices = 0;
736 double * rowValues = 0;
737
738 subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
739
740 for(int i=0;i<numEntries;i++)
741 rowValues[i] = 0.0;
742 }
743
744 const ScalarT scatterField = field_h(worksetCellIndex,basisId);
745
746 if(r_!=Teuchos::null)
747 local_r[lid] = scatterField.val();
748 local_dc[lid] = 1.0; // mark row as dirichlet
749
750 // loop over the sensitivity indices: all DOFs on a cell
751 std::vector<double> jacRow(scatterField.size(),0.0);
752
753 for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
754 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
755
756 for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
757 int start = blockOffsets[colIndexer];
758 int end = blockOffsets[colIndexer+1];
759
760 if(end-start<=0)
761 continue;
762
763 auto subColIndexer = colIndexers_[colIndexer];
764
765 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
766 auto cLIDs_h = Kokkos::create_mirror_view(cLIDs);
767 Kokkos::deep_copy(cLIDs_h, cLIDs);
768
769 TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
770
771 // check hash table for jacobian sub block
772 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
773 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
774
775 // if you didn't find one before, add it to the hash table
776 if(subJac==Teuchos::null) {
777 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
778
779 // block operator is null, don't do anything (it is excluded)
780 if(Teuchos::is_null(tOp))
781 continue;
782
783 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
784 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
785 jacEpetraBlocks[blockIndex] = subJac;
786 }
787
788 // Sum Jacobian
789 int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs_h[0]);
790 if(err!=0) {
791 std::stringstream ss;
792 ss << "Failed inserting row: " << " (" << lid << "): ";
793 for(int i=0;i<end-start;i++)
794 ss << cLIDs_h[i] << " ";
795 ss << std::endl;
796 ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
797
798 ss << "scatter field = ";
799 scatterFields_[fieldIndex].print(ss);
800 ss << std::endl;
801
802 TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
803 }
804
805 }
806 }
807 }
808 }
809}
810
811// **********************************************************************
812
813#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.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer > > &ugis, std::vector< int > &blockOffsets)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer > > &ugis)