Panzer  Version of the Day
Panzer_ScatterDirichletResidual_BlockedEpetra_Hessian_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_ScatterDirichletResidual_BlockedEpetra_Hessian_impl_hpp__
44 #define __Panzer_ScatterDirichletResidual_BlockedEpetra_Hessian_impl_hpp__
45 
46 // only do this if required by the user
47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
48 
49 // the includes for this file come in as a result of the includes in the main
50 // blocked Epetra scatter dirichlet residual file
51 
52 
53 namespace panzer {
54 
55 // **************************************************************
56 // Hessian Specialization
57 // **************************************************************
58 template<typename TRAITS,typename LO,typename GO>
60 ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & rIndexers,
61  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & cIndexers,
62  const Teuchos::ParameterList& p,
63  bool useDiscreteAdjoint)
64  : rowIndexers_(rIndexers)
65  , colIndexers_(cIndexers)
66  , globalDataKey_("Residual Scatter Container")
67 {
68  std::string scatterName = p.get<std::string>("Scatter Name");
69  scatterHolder_ =
70  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
71 
72  // get names to be evaluated
73  const std::vector<std::string>& names =
74  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
75 
76  // grab map from evaluated names to field names
77  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
78 
80  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
81 
82  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
83  local_side_id_ = p.get<int>("Local Side ID");
84 
85  // build the vector of fields that this is dependent on
86  scatterFields_.resize(names.size());
87  for (std::size_t eq = 0; eq < names.size(); ++eq) {
88  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
89 
90  // tell the field manager that we depend on this field
91  this->addDependentField(scatterFields_[eq]);
92  }
93 
94  checkApplyBC_ = p.get<bool>("Check Apply BC");
95  if (checkApplyBC_) {
96  applyBC_.resize(names.size());
97  for (std::size_t eq = 0; eq < names.size(); ++eq) {
98  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
99  this->addDependentField(applyBC_[eq]);
100  }
101  }
102 
103  // this is what this evaluator provides
104  this->addEvaluatedField(*scatterHolder_);
105 
106  if (p.isType<std::string>("Global Data Key"))
107  globalDataKey_ = p.get<std::string>("Global Data Key");
108 
109  if(colIndexers_.size()==0)
110  colIndexers_ = rowIndexers_;
111 
112  this->setName(scatterName+" Scatter Dirichlet Residual : BlockedEpetra (Hessian)");
113 }
114 
115 template<typename TRAITS,typename LO,typename GO>
116 void
118 postRegistrationSetup(typename TRAITS::SetupData d,
120 {
121  indexerIds_.resize(scatterFields_.size());
122  subFieldIds_.resize(scatterFields_.size());
123 
124  // load required field numbers for fast use
125  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
126  // get field ID from DOF manager
127  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
128 
129  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
130  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
131 
132  // fill field data object
133  this->utils.setFieldData(scatterFields_[fd],fm);
134 
135  if (checkApplyBC_)
136  this->utils.setFieldData(applyBC_[fd],fm);
137  }
138 
139  // get the number of nodes (Should be renamed basis)
140  num_nodes = scatterFields_[0].dimension(1);
141  num_eq = scatterFields_.size();
142 }
143 
144 template<typename TRAITS,typename LO,typename GO>
145 void
146 ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
147 preEvaluate(typename TRAITS::PreEvalData d)
148 {
149  typedef BlockedEpetraLinearObjContainer BLOC;
150 
151  using Teuchos::rcp_dynamic_cast;
152 
153  // extract dirichlet counter from container
154  Teuchos::RCP<const BLOC> blockContainer
155  = rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject("Dirichlet Counter"),true);
156 
157  dirichletCounter_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
158  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
159 
160  // extract linear object container
161  blockContainer = rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject(globalDataKey_),true);
162  TEUCHOS_ASSERT(!Teuchos::is_null(blockContainer));
163 
164  Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockContainer->get_A());
165 }
166 
167 template<typename TRAITS,typename LO,typename GO>
168 void
170 evaluateFields(typename TRAITS::EvalData workset)
171 {
172  using Teuchos::RCP;
173  using Teuchos::ArrayRCP;
174  using Teuchos::ptrFromRef;
175  using Teuchos::rcp_dynamic_cast;
176 
177  using Thyra::SpmdVectorBase;
178 
179  // for convenience pull out some objects from workset
180  std::string blockId = this->wda(workset).block_id;
181  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
182 
183  int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
184 
185  std::vector<int> blockOffsets;
186  computeBlockOffsets(blockId,colIndexers_,blockOffsets);
187 
188  std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
189 
190  // NOTE: A reordering of these loops will likely improve performance
191  // The "getGIDFieldOffsets may be expensive. However the
192  // "getElementGIDs" can be cheaper. However the lookup for LIDs
193  // may be more expensive!
194 
195  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
196  int rowIndexer = indexerIds_[fieldIndex];
197  int subFieldNum = subFieldIds_[fieldIndex];
198 
199  // loop over each field to be scattered
200  Teuchos::ArrayRCP<double> local_dc;
201  rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
202  ->getNonconstLocalData(ptrFromRef(local_dc));
203 
204  auto subRowIndexer = rowIndexers_[rowIndexer];
205  auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
206  const std::vector<int> & subElmtOffset = subIndicePair.first;
207  const std::vector<int> & subBasisIdMap = subIndicePair.second;
208 
209  // scatter operation for each cell in workset
210  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
211  std::size_t cellLocalId = localCellIds[worksetCellIndex];
212 
213  const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
214 
215  // loop over basis functions
216  for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
217  int offset = subElmtOffset[basis];
218  int lid = rLIDs[offset];
219  if(lid<0) // not on this processor
220  continue;
221 
222  int basisId = subBasisIdMap[basis];
223 
224  if (checkApplyBC_)
225  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
226  continue;
227 
228  // zero out matrix row
229  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
230  int start = blockOffsets[colIndexer];
231  int end = blockOffsets[colIndexer+1];
232 
233  if(end-start<=0)
234  continue;
235 
236  // check hash table for jacobian sub block
237  std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
238  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
239  // if you didn't find one before, add it to the hash table
240  if(subJac==Teuchos::null) {
241  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
242 
243  // block operator is null, don't do anything (it is excluded)
244  if(Teuchos::is_null(tOp))
245  continue;
246 
247  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
248  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
249  jacEpetraBlocks[blockIndex] = subJac;
250  }
251 
252  int numEntries = 0;
253  int * rowIndices = 0;
254  double * rowValues = 0;
255 
256  subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
257 
258  for(int i=0;i<numEntries;i++)
259  rowValues[i] = 0.0;
260  }
261 
262  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
263 
264  local_dc[lid] = 1.0; // mark row as dirichlet
265 
266  // loop over the sensitivity indices: all DOFs on a cell
267  std::vector<double> jacRow(scatterField.size(),0.0);
268 
269  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
270  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
271 
272  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
273  int start = blockOffsets[colIndexer];
274  int end = blockOffsets[colIndexer+1];
275 
276  if(end-start<=0)
277  continue;
278 
279  auto subColIndexer = colIndexers_[colIndexer];
280  const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
281 
282  TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
283 
284  // check hash table for jacobian sub block
285  std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
286  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
287 
288  // if you didn't find one before, add it to the hash table
289  if(subJac==Teuchos::null) {
290  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
291 
292  // block operator is null, don't do anything (it is excluded)
293  if(Teuchos::is_null(tOp))
294  continue;
295 
296  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
297  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
298  jacEpetraBlocks[blockIndex] = subJac;
299  }
300 
301  // Sum Jacobian
302  int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
303  if(err!=0) {
304  std::stringstream ss;
305  ss << "Failed inserting row: " << " (" << lid << "): ";
306  for(int i=0;i<end-start;i++)
307  ss << cLIDs[i] << " ";
308  ss << std::endl;
309  ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
310 
311  ss << "scatter field = ";
312  scatterFields_[fieldIndex].print(ss);
313  ss << std::endl;
314 
315  TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
316  }
317 
318  }
319  }
320  }
321  }
322 }
323 
324 }
325 
326 // **************************************************************
327 #endif
328 
329 #endif
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis, std::vector< int > &blockOffsets)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
PHX::MDField< ScalarT > vector
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isType(const std::string &name) const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
#define TEUCHOS_ASSERT(assertion_test)
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)