Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_ScatterResidual_Epetra_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_Epetra_Hessian_impl_hpp__
44 #define __Panzer_ScatterResidual_Epetra_Hessian_impl_hpp__
45 
46 // only do this if required by the user
47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
48 
49 #include "Teuchos_RCP.hpp"
50 #include "Teuchos_Assert.hpp"
51 
52 #include "Phalanx_DataLayout.hpp"
53 
54 #include "Epetra_Map.h"
55 #include "Epetra_Vector.h"
56 #include "Epetra_CrsMatrix.h"
57 
60 #include "Panzer_PureBasis.hpp"
65 
66 #include "Phalanx_DataLayout_MDALayout.hpp"
67 
68 #include "Teuchos_FancyOStream.hpp"
69 
70 // **********************************************************************
71 // Specialization: Hessian
72 // **********************************************************************
73 
74 template<typename TRAITS,typename LO,typename GO>
77  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & cIndexer,
78  const Teuchos::ParameterList& p,
79  bool useDiscreteAdjoint)
80  : globalIndexer_(indexer)
81  , colGlobalIndexer_(cIndexer)
82  , globalDataKey_("Residual Scatter Container")
83  , useDiscreteAdjoint_(useDiscreteAdjoint)
84 {
85  std::string scatterName = p.get<std::string>("Scatter Name");
86  scatterHolder_ =
87  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
88 
89  // get names to be evaluated
90  const std::vector<std::string>& names =
91  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
92 
93  // grab map from evaluated names to field names
94  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
95 
97  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
98 
99  // build the vector of fields that this is dependent on
100  scatterFields_.resize(names.size());
101  for (std::size_t eq = 0; eq < names.size(); ++eq) {
102  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
103 
104  // tell the field manager that we depend on this field
105  this->addDependentField(scatterFields_[eq]);
106  }
107 
108  // this is what this evaluator provides
109  this->addEvaluatedField(*scatterHolder_);
110 
111  if (p.isType<std::string>("Global Data Key"))
112  globalDataKey_ = p.get<std::string>("Global Data Key");
113  if (p.isType<bool>("Use Discrete Adjoint"))
114  useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
115 
116  // discrete adjoint does not work with non-square matrices
117  if(useDiscreteAdjoint)
118  { TEUCHOS_ASSERT(colGlobalIndexer_==globalIndexer_); }
119 
120  this->setName(scatterName+" Scatter Residual Epetra (Jacobian)");
121 }
122 
123 // **********************************************************************
124 template<typename TRAITS,typename LO,typename GO>
126 postRegistrationSetup(typename TRAITS::SetupData /* d */,
127  PHX::FieldManager<TRAITS>& /* fm */)
128 {
129  fieldIds_.resize(scatterFields_.size());
130  // load required field numbers for fast use
131  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
132  // get field ID from DOF manager
133  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
134  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
135  }
136 }
137 
138 // **********************************************************************
139 template<typename TRAITS,typename LO,typename GO>
141 preEvaluate(typename TRAITS::PreEvalData d)
142 {
143  // extract linear object container
144  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
145 
146  if(epetraContainer_==Teuchos::null) {
147  // extract linear object container
148  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
149  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
150  }
151 }
152 
153 // **********************************************************************
154 template<typename TRAITS,typename LO,typename GO>
156 evaluateFields(typename TRAITS::EvalData workset)
157 {
158  using Kokkos::View;
159  using PHX::Device;
160  using std::vector;
161 
162  std::vector<double> jacRow;
163 
164  bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
165 
166  // for convenience pull out some objects from workset
167  std::string blockId = this->wda(workset).block_id;
168  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
169 
170  Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
171  Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
172 
174  colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
175 
176  // NOTE: A reordering of these loops will likely improve performance
177  // The "getGIDFieldOffsets" may be expensive. However the
178  // "getElementGIDs" can be cheaper. However the lookup for LIDs
179  // may be more expensive!
180 
181  // scatter operation for each cell in workset
182  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
183  std::size_t cellLocalId = localCellIds[worksetCellIndex];
184 
185  auto rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
186  auto initial_cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
187  vector<int> cLIDs;
188  for (int i(0); i < static_cast<int>(initial_cLIDs.extent(0)); ++i)
189  cLIDs.push_back(initial_cLIDs(i));
190  if (Teuchos::nonnull(workset.other)) {
191  const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
192  auto other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
193  for (int i(0); i < static_cast<int>(other_cLIDs.extent(0)); ++i)
194  cLIDs.push_back(other_cLIDs(i));
195  }
196 
197  // loop over each field to be scattered
198  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
199  int fieldNum = fieldIds_[fieldIndex];
200  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
201 
202  // loop over the basis functions (currently they are nodes)
203  for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
204  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
205  int rowOffset = elmtOffset[rowBasisNum];
206  int row = rLIDs[rowOffset];
207 
208  // loop over the sensitivity indices: all DOFs on a cell
209  jacRow.resize(scatterField.size());
210 
211  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
212  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
213 
214  {
215  int err = Jac->SumIntoMyValues(
216  row,
217  std::min(cLIDs.size(), static_cast<size_t>(scatterField.size())),
218  panzer::ptrFromStlVector(jacRow),
219  panzer::ptrFromStlVector(cLIDs));
221  }
222  } // end rowBasisNum
223  } // end fieldIndex
224  }
225 }
226 
227 // **********************************************************************
228 
229 #endif
230 
231 #endif
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T * ptrFromStlVector(std::vector< T > &v)
Pushes residual values into the residual vector for a Newton-based solve.
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)