Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_ScatterResidual_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_ScatterResidual_BlockedEpetra_Hessian_impl_hpp__
44 #define __Panzer_ScatterResidual_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 // Epetra scatter residual file
51 
52 namespace panzer {
53 
54 // **************************************************************
55 // Hessian Specialization
56 // **************************************************************
57 template<typename TRAITS,typename LO,typename GO>
59 ScatterResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer<LO,int> > > & rIndexers,
60  const std::vector<Teuchos::RCP<const GlobalIndexer<LO,int> > > & cIndexers,
61  const Teuchos::ParameterList& p,
62  bool useDiscreteAdjoint)
63  : rowIndexers_(rIndexers)
64  , colIndexers_(cIndexers)
65  , globalDataKey_("Residual Scatter Container")
66  , useDiscreteAdjoint_(useDiscreteAdjoint)
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<const panzer::PureBasis> >("Basis")->functional;
81 
82  // build the vector of fields that this is dependent on
83  scatterFields_.resize(names.size());
84  for (std::size_t eq = 0; eq < names.size(); ++eq) {
85  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
86 
87  // tell the field manager that we depend on this field
88  this->addDependentField(scatterFields_[eq]);
89  }
90 
91  // this is what this evaluator provides
92  this->addEvaluatedField(*scatterHolder_);
93 
94  if (p.isType<std::string>("Global Data Key"))
95  globalDataKey_ = p.get<std::string>("Global Data Key");
96  if (p.isType<bool>("Use Discrete Adjoint"))
97  useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
98 
99  // discrete adjoint does not work with non-square matrices
100  if(useDiscreteAdjoint)
101  { TEUCHOS_ASSERT(colIndexers_.size()==0); }
102 
103  if(colIndexers_.size()==0)
104  colIndexers_ = rowIndexers_;
105 
106  this->setName(scatterName+" Scatter Residual BlockedEpetra (Hessian)");
107 }
108 
109 template<typename TRAITS,typename LO,typename GO>
110 void
112 postRegistrationSetup(typename TRAITS::SetupData /* d */,
113  PHX::FieldManager<TRAITS>& /* fm */)
114 {
115  indexerIds_.resize(scatterFields_.size());
116  subFieldIds_.resize(scatterFields_.size());
117 
118  // load required field numbers for fast use
119  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
120  // get field ID from DOF manager
121  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
122 
123  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
124  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
125  }
126 }
127 
128 template<typename TRAITS,typename LO,typename GO>
129 void
131 preEvaluate(typename TRAITS::PreEvalData d)
132 {
133  using Teuchos::RCP;
134  using Teuchos::rcp_dynamic_cast;
135 
136  typedef BlockedEpetraLinearObjContainer BLOC;
137  typedef BlockedEpetraLinearObjContainer ELOC;
138 
139  // extract linear object container
140  RCP<const BLOC> blockedContainer = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
141  RCP<const ELOC> epetraContainer = rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
142 
143  // if its blocked do this
144  if(blockedContainer!=Teuchos::null) {
145  Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockedContainer->get_A());
146  }
147  else if(epetraContainer!=Teuchos::null) {
148  // convert it into a blocked operator
149  RCP<Thyra::LinearOpBase<double> > J = blockedContainer->get_A_th();
150  Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(Thyra::nonconstBlock1x1(J));
151  }
152 
153  TEUCHOS_ASSERT(Jac_!=Teuchos::null);
154 }
155 
156 template<typename TRAITS,typename LO,typename GO>
157 void
159 evaluateFields(typename TRAITS::EvalData workset)
160 {
161  using Teuchos::RCP;
162  using Teuchos::ArrayRCP;
163  using Teuchos::ptrFromRef;
164  using Teuchos::rcp_dynamic_cast;
165 
166  using Thyra::VectorBase;
167  using Thyra::SpmdVectorBase;
170 
171  std::vector<double> jacRow;
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  int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
178 
179  std::vector<int> blockOffsets;
180  computeBlockOffsets(blockId,colIndexers_,blockOffsets);
181 
182  std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
183 
184  // loop over each field to be scattered
185  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
186  int rowIndexer = indexerIds_[fieldIndex];
187  int subFieldNum = subFieldIds_[fieldIndex];
188 
189  auto subRowIndexer = rowIndexers_[rowIndexer];
190  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
191 
192  // scatter operation for each cell in workset
193  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
194  std::size_t cellLocalId = localCellIds[worksetCellIndex];
195 
196  auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
197 
198  // loop over the basis functions (currently they are nodes)
199  for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
200  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
201  int rowOffset = elmtOffset[rowBasisNum];
202  int r_lid = rLIDs[rowOffset];
203 
204  // loop over the sensitivity indices: all DOFs on a cell
205  jacRow.resize(scatterField.size());
206 
207  // For Neumann conditions with no dependence on degrees of freedom, there should be no Jacobian contribution
208  if(scatterField.size() == 0)
209  continue;
210 
211  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
212  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
213 
214  // scatter the row to each block
215  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
216  int start = blockOffsets[colIndexer];
217  int end = blockOffsets[colIndexer+1];
218 
219  if(end-start<=0)
220  continue;
221 
222  auto subColIndexer = colIndexers_[colIndexer];
223  auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
224 
225  TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
226 
227  // check hash table for jacobian sub block
228  std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
229  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
230 
231  // if you didn't find one before, add it to the hash table
232  if(subJac==Teuchos::null) {
233  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
234 
235  // block operator is null, don't do anything (it is excluded)
236  if(Teuchos::is_null(tOp))
237  continue;
238 
239  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
240  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
241  jacEpetraBlocks[blockIndex] = subJac;
242  }
243 
244  // Sum Jacobian
245  {
246  int err = subJac->SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs[0]);
247  if(err!=0) {
248 
249  std::stringstream ss;
250  ss << "Failed inserting row: " << "LID = " << r_lid << ": ";
251  for(int i=0;i<end-start;i++)
252  ss << cLIDs[i] << " ";
253  ss << std::endl;
254  ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
255 
256  ss << "scatter field = ";
257  scatterFields_[fieldIndex].print(ss);
258  ss << std::endl;
259 
260  ss << "values = ";
261  for(int i=start;i<end;i++)
262  ss << jacRow[i] << " ";
263  ss << std::endl;
264 
265  std::cout << ss.str() << std::endl;
266 
267  TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
268  }
269  }
270  }
271  } // end rowBasisNum
272  } // end fieldIndex
273  }
274 }
275 
276 }
277 
278 // **************************************************************
279 #endif
280 
281 #endif
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer >> &ugis, std::vector< int > &blockOffsets)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer >> &ugis)
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Pushes residual values into the residual vector for a Newton-based solve.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
TransListIter end