Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterResidual_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_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP
44#define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_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
57#include "Panzer_PureBasis.hpp"
60#include "Panzer_HashUtils.hpp"
62
63#include "Thyra_SpmdVectorBase.hpp"
64#include "Thyra_ProductVectorBase.hpp"
65#include "Thyra_DefaultProductVector.hpp"
66#include "Thyra_BlockedLinearOpBase.hpp"
67#include "Thyra_DefaultBlockedLinearOp.hpp"
68#include "Thyra_get_Epetra_Operator.hpp"
69
70#include "Phalanx_DataLayout_MDALayout.hpp"
71
72#include "Teuchos_FancyOStream.hpp"
73
74// **********************************************************************
75// Specialization: Residual
76// **********************************************************************
77
78template<typename TRAITS,typename LO,typename GO>
80ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
81 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
82 const Teuchos::ParameterList& p,
83 bool /* useDiscreteAdjoint */)
84 : rowIndexers_(rIndexers)
85 , colIndexers_(cIndexers)
86 , globalDataKey_("Residual Scatter Container")
87{
88 std::string scatterName = p.get<std::string>("Scatter Name");
89 scatterHolder_ =
90 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
91
92 // get names to be evaluated
93 const std::vector<std::string>& names =
94 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
95
96 // grab map from evaluated names to field names
97 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
98
99 Teuchos::RCP<PHX::DataLayout> dl =
100 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
101
102 // build the vector of fields that this is dependent on
103 scatterFields_.resize(names.size());
104 for (std::size_t eq = 0; eq < names.size(); ++eq) {
105 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
106
107 // tell the field manager that we depend on this field
108 this->addDependentField(scatterFields_[eq]);
109 }
110
111 // this is what this evaluator provides
112 this->addEvaluatedField(*scatterHolder_);
113
114 if (p.isType<std::string>("Global Data Key"))
115 globalDataKey_ = p.get<std::string>("Global Data Key");
116
117 this->setName(scatterName+" Scatter Residual");
118}
119
120// **********************************************************************
121template<typename TRAITS,typename LO,typename GO>
123postRegistrationSetup(typename TRAITS::SetupData /* d */,
125{
126 indexerIds_.resize(scatterFields_.size());
127 subFieldIds_.resize(scatterFields_.size());
128
129 // load required field numbers for fast use
130 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
131 // get field ID from DOF manager
132 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
133
134 indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
135 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
136 }
137}
138
139// **********************************************************************
140template<typename TRAITS,typename LO,typename GO>
142preEvaluate(typename TRAITS::PreEvalData d)
143{
146
147 // extract linear object container
148 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
149 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
150
151 // if its blocked do this
152 if(blockedContainer!=Teuchos::null)
153 r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
154 else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
155 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
156
157 TEUCHOS_ASSERT(r_!=Teuchos::null);
158}
159
160// **********************************************************************
161template<typename TRAITS,typename LO,typename GO>
163evaluateFields(typename TRAITS::EvalData workset)
164{
165 using Teuchos::RCP;
166 using Teuchos::ptrFromRef;
167 using Teuchos::rcp_dynamic_cast;
168
169 using Thyra::VectorBase;
170 using Thyra::SpmdVectorBase;
172
173 // for convenience pull out some objects from workset
174 std::string blockId = this->wda(workset).block_id;
175 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
176
177 // NOTE: A reordering of these loops will likely improve performance
178 // The "getGIDFieldOffsets may be expensive. However the
179 // "getElementGIDs" can be cheaper. However the lookup for LIDs
180 // may be more expensive!
181
182 // loop over each field to be scattered
183 Teuchos::ArrayRCP<double> local_r;
184 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
185 int indexerId = indexerIds_[fieldIndex];
186 int subFieldNum = subFieldIds_[fieldIndex];
187
188 // grab local data for inputing
189 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
190
191 auto subRowIndexer = rowIndexers_[indexerId];
192 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
193
194 auto field = PHX::as_view(scatterFields_[fieldIndex]);
195 auto field_h = Kokkos::create_mirror_view(field);
196 Kokkos::deep_copy(field_h, field);
197
198 // scatter operation for each cell in workset
199 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
200 std::size_t cellLocalId = localCellIds[worksetCellIndex];
201
202 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
203 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
204 Kokkos::deep_copy(LIDs_h, LIDs);
205
206 // loop over basis functions
207 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
208 int offset = elmtOffset[basis];
209 int lid = LIDs_h[offset];
210 local_r[lid] += field_h(worksetCellIndex,basis);
211 }
212 }
213 }
214}
215
216// **********************************************************************
217// Specialization: Tangent
218// **********************************************************************
219
220template<typename TRAITS,typename LO,typename GO>
222ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
223 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
224 const Teuchos::ParameterList& p,
225 bool /* useDiscreteAdjoint */)
226 : rowIndexers_(rIndexers)
227 , colIndexers_(cIndexers)
228 , globalDataKey_("Residual Scatter Container")
229{
230 std::string scatterName = p.get<std::string>("Scatter Name");
231 scatterHolder_ =
232 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
233
234 // get names to be evaluated
235 const std::vector<std::string>& names =
236 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
237
238 // grab map from evaluated names to field names
239 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
240
241 Teuchos::RCP<PHX::DataLayout> dl =
242 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
243
244 // build the vector of fields that this is dependent on
245 scatterFields_.resize(names.size());
246 for (std::size_t eq = 0; eq < names.size(); ++eq) {
247 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
248
249 // tell the field manager that we depend on this field
250 this->addDependentField(scatterFields_[eq]);
251 }
252
253 // this is what this evaluator provides
254 this->addEvaluatedField(*scatterHolder_);
255
256 if (p.isType<std::string>("Global Data Key"))
257 globalDataKey_ = p.get<std::string>("Global Data Key");
258
259 this->setName(scatterName+" Scatter Tangent");
260}
261
262// **********************************************************************
263template<typename TRAITS,typename LO,typename GO>
265postRegistrationSetup(typename TRAITS::SetupData /* d */,
267{
268 indexerIds_.resize(scatterFields_.size());
269 subFieldIds_.resize(scatterFields_.size());
270
271 // load required field numbers for fast use
272 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
273 // get field ID from DOF manager
274 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
275
276 indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
277 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
278 }
279}
280
281// **********************************************************************
282template<typename TRAITS,typename LO,typename GO>
284preEvaluate(typename TRAITS::PreEvalData d)
285{
288
289 // extract linear object container
290 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
291 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
292
293 // if its blocked do this
294 if(blockedContainer!=Teuchos::null)
295 r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
296 else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
297 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
298
299 TEUCHOS_ASSERT(r_!=Teuchos::null);
300}
301
302// **********************************************************************
303template<typename TRAITS,typename LO,typename GO>
305evaluateFields(typename TRAITS::EvalData workset)
306{
307 TEUCHOS_ASSERT(false);
308
309 using Teuchos::RCP;
310 using Teuchos::ptrFromRef;
311 using Teuchos::rcp_dynamic_cast;
312
313 using Thyra::VectorBase;
314 using Thyra::SpmdVectorBase;
316
317 // for convenience pull out some objects from workset
318 std::string blockId = this->wda(workset).block_id;
319 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
320
321 // NOTE: A reordering of these loops will likely improve performance
322 // The "getGIDFieldOffsets may be expensive. However the
323 // "getElementGIDs" can be cheaper. However the lookup for LIDs
324 // may be more expensive!
325
326 // loop over each field to be scattered
327 Teuchos::ArrayRCP<double> local_r;
328 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
329 int indexerId = indexerIds_[fieldIndex];
330 int subFieldNum = subFieldIds_[fieldIndex];
331
332 // grab local data for inputing
333 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
334
335 auto subRowIndexer = rowIndexers_[indexerId];
336 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
337
338 // scatter operation for each cell in workset
339 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
340 std::size_t cellLocalId = localCellIds[worksetCellIndex];
341
342 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
343
344 // loop over basis functions
345 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
346 int offset = elmtOffset[basis];
347 int lid = LIDs[offset];
348 local_r[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
349 }
350 }
351 }
352}
353
354// **********************************************************************
355// Specialization: Jacobian
356// **********************************************************************
357
358template<typename TRAITS,typename LO,typename GO>
360ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
361 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
362 const Teuchos::ParameterList& p,
363 bool useDiscreteAdjoint)
364 : rowIndexers_(rIndexers)
365 , colIndexers_(cIndexers)
366 , globalDataKey_("Residual Scatter Container")
367 , useDiscreteAdjoint_(useDiscreteAdjoint)
368{
369 std::string scatterName = p.get<std::string>("Scatter Name");
370 scatterHolder_ =
371 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
372
373 // get names to be evaluated
374 const std::vector<std::string>& names =
375 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
376
377 // grab map from evaluated names to field names
378 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
379
380 Teuchos::RCP<PHX::DataLayout> dl =
381 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
382
383 // build the vector of fields that this is dependent on
384 scatterFields_.resize(names.size());
385 for (std::size_t eq = 0; eq < names.size(); ++eq) {
386 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
387
388 // tell the field manager that we depend on this field
389 this->addDependentField(scatterFields_[eq]);
390 }
391
392 // this is what this evaluator provides
393 this->addEvaluatedField(*scatterHolder_);
394
395 if (p.isType<std::string>("Global Data Key"))
396 globalDataKey_ = p.get<std::string>("Global Data Key");
397 if (p.isType<bool>("Use Discrete Adjoint"))
398 useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
399
400 // discrete adjoint does not work with non-square matrices
401 if(useDiscreteAdjoint)
402 { TEUCHOS_ASSERT(colIndexers_.size()==0); }
403
404 if(colIndexers_.size()==0)
405 colIndexers_ = rowIndexers_;
406
407 this->setName(scatterName+" Scatter Residual BlockedEpetra (Jacobian)");
408}
409
410// **********************************************************************
411template<typename TRAITS,typename LO,typename GO>
413postRegistrationSetup(typename TRAITS::SetupData /* d */,
415{
416 indexerIds_.resize(scatterFields_.size());
417 subFieldIds_.resize(scatterFields_.size());
418
419 // load required field numbers for fast use
420 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
421 // get field ID from DOF manager
422 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
423
424 indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
425 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
426 }
427}
428
429// **********************************************************************
430template<typename TRAITS,typename LO,typename GO>
432preEvaluate(typename TRAITS::PreEvalData d)
433{
434 using Teuchos::RCP;
435 using Teuchos::rcp_dynamic_cast;
436
439
440 // extract linear object container
441 RCP<const BLOC> blockedContainer = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
442 RCP<const ELOC> epetraContainer = rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
443
444 // if its blocked do this
445 if(blockedContainer!=Teuchos::null) {
446 r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
447 Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockedContainer->get_A());
448 }
449 else if(epetraContainer!=Teuchos::null) {
450 // if its straight up epetra do this
451 if(epetraContainer->get_f_th()!=Teuchos::null)
452 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
453
454 // convert it into a blocked operator
455 RCP<Thyra::LinearOpBase<double> > J = blockedContainer->get_A_th();
456 Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(Thyra::nonconstBlock1x1(J));
457 }
458
459 TEUCHOS_ASSERT(Jac_!=Teuchos::null);
460}
461
462// **********************************************************************
463template<typename TRAITS,typename LO,typename GO>
465evaluateFields(typename TRAITS::EvalData workset)
466{
467 using Teuchos::RCP;
468 using Teuchos::ArrayRCP;
469 using Teuchos::ptrFromRef;
470 using Teuchos::rcp_dynamic_cast;
471
472 using Thyra::VectorBase;
473 using Thyra::SpmdVectorBase;
476
477 std::vector<double> jacRow;
478
479 // for convenience pull out some objects from workset
480 std::string blockId = this->wda(workset).block_id;
481 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
482
483 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
484
485 std::vector<int> blockOffsets;
486 computeBlockOffsets(blockId,colIndexers_,blockOffsets);
487
488 std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
489
490 // loop over each field to be scattered
491 Teuchos::ArrayRCP<double> local_r;
492 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
493 int rowIndexer = indexerIds_[fieldIndex];
494 int subFieldNum = subFieldIds_[fieldIndex];
495
496 // grab local data for inputing
497 if(r_!=Teuchos::null)
498 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))->getNonconstLocalData(ptrFromRef(local_r));
499
500 auto subRowIndexer = rowIndexers_[rowIndexer];
501 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
502
503 auto field = scatterFields_[fieldIndex].get_view();
504 auto field_h = Kokkos::create_mirror_view(field);
505 Kokkos::deep_copy(field_h, field);
506
507 auto rLIDs = subRowIndexer->getLIDs();
508 auto rLIDs_h = Kokkos::create_mirror_view(rLIDs);
509 Kokkos::deep_copy(rLIDs_h, rLIDs);
510
511 // scatter operation for each cell in workset
512 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
513 std::size_t cellLocalId = localCellIds[worksetCellIndex];
514
515 // loop over the basis functions (currently they are nodes)
516 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
517 const ScalarT scatterField = field_h(worksetCellIndex,rowBasisNum);
518 int rowOffset = elmtOffset[rowBasisNum];
519 int r_lid = rLIDs_h(cellLocalId, rowOffset);
520
521 // Sum residual
522 if(local_r!=Teuchos::null)
523 local_r[r_lid] += (scatterField.val());
524
525 // loop over the sensitivity indices: all DOFs on a cell
526 jacRow.resize(scatterField.size());
527
528 // For Neumann conditions with no dependence on degrees of freedom, there should be no Jacobian contribution
529 if(scatterField.size() == 0)
530 continue;
531
532 for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
533 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
534
535 // scatter the row to each block
536 for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
537 int start = blockOffsets[colIndexer];
538 int end = blockOffsets[colIndexer+1];
539
540 if(end-start<=0)
541 continue;
542
543 auto subColIndexer = colIndexers_[colIndexer];
544 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
545 auto cLIDs_h = Kokkos::create_mirror_view(cLIDs);
546 Kokkos::deep_copy(cLIDs_h, cLIDs);
547
548 TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
549
550 // check hash table for jacobian sub block
551 std::pair<int,int> blockIndex = useDiscreteAdjoint_ ? std::make_pair(colIndexer,rowIndexer)
552 : std::make_pair(rowIndexer,colIndexer);
553 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
554
555 // if you didn't find one before, add it to the hash table
556 if(subJac==Teuchos::null) {
557 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
558
559 // block operator is null, don't do anything (it is excluded)
560 if(Teuchos::is_null(tOp))
561 continue;
562
563 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
564 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
565 jacEpetraBlocks[blockIndex] = subJac;
566 }
567
568 // Sum Jacobian
569 if(!useDiscreteAdjoint_) {
570 int err = subJac->SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs_h[0]);
571 if(err!=0) {
572
573 std::stringstream ss;
574 ss << "Failed inserting row: " << "LID = " << r_lid << "): ";
575 for(int i=start;i<end;i++)
576 ss << cLIDs_h[i] << " ";
577 ss << std::endl;
578 ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
579
580 ss << "scatter field = ";
581 scatterFields_[fieldIndex].print(ss);
582 ss << std::endl;
583
584 TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
585 }
586 }
587 else {
588 for(std::size_t c=0;c<cLIDs.size();c++) {
589 int err = subJac->SumIntoMyValues(cLIDs_h[c], 1, &jacRow[start+c],&r_lid);
590 TEUCHOS_ASSERT_EQUALITY(err,0);
591 }
592 }
593 }
594 } // end rowBasisNum
595 } // end fieldIndex
596 }
597}
598
599// **********************************************************************
600
601#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)