Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_GatherSolution_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_GatherSolution_Epetra_Hessian_impl_hpp__
44 #define __Panzer_GatherSolution_Epetra_Hessian_impl_hpp__
45 
46 // Only do this if required by the user.
47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
48 
50 //
51 // Include Files
52 //
54 
55 // Epetra
56 #include "Epetra_Map.h"
57 #include "Epetra_Vector.h"
58 
59 // Panzer
64 #include "Panzer_PureBasis.hpp"
66 
67 // Phalanx
68 #include "Phalanx_DataLayout.hpp"
69 
70 // Teuchos
71 #include "Teuchos_Assert.hpp"
72 #include "Teuchos_FancyOStream.hpp"
73 
74 // Thyra
75 #include "Thyra_SpmdVectorBase.hpp"
76 
78 //
79 // Initializing Constructor (Hessian Specialization)
80 //
82 template<typename TRAITS, typename LO, typename GO>
86  const Teuchos::ParameterList& p)
87  :
88  globalIndexer_(indexer)
89 {
90  using panzer::PureBasis;
91  using PHX::MDField;
92  using PHX::typeAsString;
93  using std::size_t;
94  using std::string;
95  using std::vector;
96  using Teuchos::RCP;
98  input.setParameterList(p);
99  const vector<string>& names = input.getDofNames();
100  RCP<const PureBasis> basis = input.getBasis();
101  indexerNames_ = input.getIndexerNames();
102  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
103  globalDataKey_ = input.getGlobalDataKey();
104  gatherSeedIndex_ = input.getGatherSeedIndex();
105  sensitivitiesName_ = input.getSensitivitiesName();
106  firstSensitivitiesAvailable_ = input.firstSensitivitiesAvailable();
107  secondSensitivitiesAvailable_ = input.secondSensitivitiesAvailable();
108  sensitivities2ndPrefix_ = input.getSecondSensitivityDataKeyPrefix();
109 
110  // Allocate the fields.
111  int numFields(names.size());
112  gatherFields_.resize(numFields);
113  for (int fd(0); fd < numFields; ++fd)
114  {
115  gatherFields_[fd] =
116  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
117  this->addEvaluatedField(gatherFields_[fd]);
118  } // end loop over names
119 
120  // Figure out what the first active name is.
121  string firstName("<none>"), n("Gather Solution (Epetra");
122  if (numFields > 0)
123  firstName = names[0];
124  if (not firstSensitivitiesAvailable_)
125  n += ", No First Sensitivities";
126  n += "): " + firstName + " (" + typeAsString<EvalT>() + ") ";
127  this->setName(n);
128 } // end of Initializing Constructor (Hessian Specialization)
129 
131 //
132 // postRegistrationSetup() (Hessian Specialization)
133 //
135 template<typename TRAITS, typename LO, typename GO>
136 void
139  typename TRAITS::SetupData /* d */,
140  PHX::FieldManager<TRAITS>& /* fm */)
141 {
142  using std::logic_error;
143  using std::size_t;
144  using std::string;
145  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
146  int numFields(gatherFields_.size());
147  fieldIds_.resize(numFields);
148  for (int fd(0); fd < numFields; ++fd)
149  {
150  // Get the field ID from the DOF manager.
151  const string& fieldName(indexerNames_[fd]);
152  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
153 
154  // This is the error return code; raise the alarm.
155  TEUCHOS_TEST_FOR_EXCEPTION(fieldIds_[fd] == -1, logic_error,
156  "GatherSolution_Epetra<Hessian>: Could not find field \"" + fieldName
157  + "\" in the global indexer. ")
158  } // end loop over gatherFields_
159  indexerNames_.clear();
160 } // end of postRegistrationSetup() (Hessian Specialization)
161 
163 //
164 // preEvaluate() (Hessian Specialization)
165 //
167 template<typename TRAITS, typename LO, typename GO>
168 void
171  typename TRAITS::PreEvalData d)
172 {
173  using std::logic_error;
174  using std::string;
175  using Teuchos::RCP;
176  using Teuchos::rcp_dynamic_cast;
179  using GED = panzer::GlobalEvaluationData;
180  using LOC = panzer::LinearObjContainer;
182 
183  // Manage sensitivities.
184  firstApplySensitivities_ = false;
185  if ((firstSensitivitiesAvailable_ ) and
186  (d.first_sensitivities_name == sensitivitiesName_))
187  firstApplySensitivities_ = true;
188  secondApplySensitivities_ = false;
189  if ((secondSensitivitiesAvailable_ ) and
190  (d.second_sensitivities_name == sensitivitiesName_))
191  secondApplySensitivities_ = true;
192 
193  // First try refactored ReadOnly container.
194  RCP<GED> ged;
195  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
196  if (d.gedc->containsDataObject(globalDataKey_ + post))
197  {
198  ged = d.gedc->getDataObject(globalDataKey_ + post);
199  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
200  }
201 
202  // Otherwise, try the old path.
203  else if (d.gedc->containsDataObject(globalDataKey_))
204  {
205  ged = d.gedc->getDataObject(globalDataKey_);
206 
207  // Try to extract the linear object container.
208  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged);
209  auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
210 
211  // Handle the LOCPair case.
212  {
213  auto locPair = rcp_dynamic_cast<LPGED>(ged);
214  if (not locPair.is_null())
215  {
216  RCP<LOC> loc = locPair->getGhostedLOC();
217  epetraContainer = rcp_dynamic_cast<ELOC>(loc);
218  } // end if (not locPair.is_null())
219  } // end of the LOCPair case
220  if ((xEvRoGed_.is_null()) and (not epetraContainer.is_null()))
221  {
222  if (useTimeDerivativeSolutionVector_)
223  x_ = epetraContainer->get_dxdt();
224  else // if (not useTimeDerivativeSolutionVector_)
225  x_ = epetraContainer->get_x();
226  } // end if ((xEvRoGed_.is_null()) and (not epetraContainer.is_null()))
227  } // end if we're doing things the new or old way
228 
229  // Ensure that we actually have something.
230  TEUCHOS_TEST_FOR_EXCEPTION((x_.is_null()) and (xEvRoGed_.is_null()),
231  logic_error, "GatherSolution_Epetra_Hessian::preEvaluate(): Unable to " \
232  "find solution vector.")
233 
234  // Don't try to extract dx if it's not required.
235  if (not secondApplySensitivities_)
236  return;
237 
238  // Now parse the second derivative direction.
239  if (d.gedc->containsDataObject(sensitivities2ndPrefix_ + globalDataKey_))
240  {
241  ged = d.gedc->getDataObject(sensitivities2ndPrefix_ + globalDataKey_);
242  dxEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
243 
244  // Ensure that we actually have something.
245  TEUCHOS_TEST_FOR_EXCEPTION(dxEvRoGed_.is_null(), logic_error, "Cannot " \
246  "find sensitivity vector associated with \"" + sensitivities2ndPrefix_ +
247  globalDataKey_ + "\" and \"" + post + "\".")
248  } // end of parsing the second derivative direction
249 } // end of preEvaluate() (Hessian Specialization)
250 
252 //
253 // evaluateFields() (Hessian Specialization)
254 //
256 template<typename TRAITS, typename LO, typename GO>
257 void
259 evaluateFields(
260  typename TRAITS::EvalData workset)
261 {
262  using PHX::MDField;
263  using std::size_t;
264  using std::string;
265  using std::vector;
266  using Teuchos::ArrayRCP;
267  using Teuchos::ptrFromRef;
268  using Teuchos::RCP;
269  using Teuchos::rcp_dynamic_cast;
270  using Thyra::SpmdVectorBase;
271 
272  // For convenience, pull out some objects from the workset.
273  string blockId(this->wda(workset).block_id);
274  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
275  int numFields(gatherFields_.size()), numCells(localCellIds.size());
276 
277  // Set a sensitivity seed value.
278  double seedValue(0);
279  if (firstApplySensitivities_)
280  {
281  if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
282  seedValue = workset.alpha;
283  else if (gatherSeedIndex_ < 0)
284  seedValue = workset.beta;
285  else if (not useTimeDerivativeSolutionVector_)
286  seedValue = workset.gather_seeds[gatherSeedIndex_];
287  else
288  TEUCHOS_ASSERT(false);
289  } // end if (firstApplySensitivities_)
290 
291  // Turn off sensitivies. This may be faster if we don't expand the term, but
292  // I suspect not, because anywhere it is used the full complement of
293  // sensitivies will be needed anyway.
294  if (not firstApplySensitivities_)
295  seedValue = 0.0;
296 
297  // Interface worksets handle DOFs from two element blocks. The derivative
298  // offset for the other element block must be shifted by the derivative side
299  // of my element block.
300  int dos(0);
301  if (this->wda.getDetailsIndex() == 1)
302  {
303  // Get the DOF count for my element block.
304  dos = globalIndexer_->getElementBlockGIDCount(workset.details(0).block_id);
305  } // end if (this->wda.getDetailsIndex() == 1)
306 
307  if (x_.is_null())
308  {
309  // Loop over the fields to be gathered.
310  for (int fieldIndex(0); fieldIndex < numFields; ++fieldIndex)
311  {
312  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldIndex];
313  int fieldNum(fieldIds_[fieldIndex]);
314  const vector<int>& elmtOffset =
315  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
316  int numBases(elmtOffset.size());
317 
318  // Gather operation for each cell in the workset.
319  for (int cell(0); cell < numCells; ++cell)
320  {
321  size_t cellLocalId(localCellIds[cell]);
322  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
323 
324  // Loop over the basis functions and fill the fields.
325  for (int basis(0); basis < numBases; ++basis)
326  {
327  int offset(elmtOffset[basis]), lid(LIDs[offset]);
328  field(cell, basis) = (*xEvRoGed_)[lid];
329  } // end loop over the basis functions
330  } // end loop over localCellIds
331  } // end loop over the fields to be gathered
332  }
333  else // if (not x_.is_null())
334  {
335  // Loop over the fields to be gathered.
336  for (int fieldIndex(0); fieldIndex < numFields; ++fieldIndex)
337  {
338  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldIndex];
339  int fieldNum(fieldIds_[fieldIndex]);
340  const vector<int>& elmtOffset =
341  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
342  int numBases(elmtOffset.size());
343 
344  // Gather operation for each cell in the workset.
345  for (int cell(0); cell < numCells; ++cell)
346  {
347  size_t cellLocalId(localCellIds[cell]);
348  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
349 
350  // Loop over the basis functions and fill the fields.
351  for (int basis(0); basis < numBases; ++basis)
352  {
353  int offset(elmtOffset[basis]), lid(LIDs[offset]);
354  field(cell, basis) = (*x_)[lid];
355  } // end loop over the basis functions
356  } // end loop over localCellIds
357  } // end loop over the fields to be gathered
358  } // end if (x_.is_null()) or not
359 
360  // Deal with the first sensitivities.
361  if (firstApplySensitivities_)
362  {
363  // Loop over the fields to be gathered.
364  for (int fieldIndex(0); fieldIndex < numFields; ++fieldIndex)
365  {
366  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldIndex];
367  int fieldNum(fieldIds_[fieldIndex]);
368  const vector<int>& elmtOffset =
369  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
370  int numBases(elmtOffset.size());
371 
372  // Gather operation for each cell in the workset.
373  for (int cell(0); cell < numCells; ++cell)
374  {
375  // Loop over the basis functions and fill the fields.
376  for (int basis(0); basis < numBases; ++basis)
377  {
378  int offset(elmtOffset[basis]);
379  field(cell, basis).fastAccessDx(dos + offset) = seedValue;
380  } // end loop over the basis functions
381  } // end loop over localCellIds
382  } // end loop over the fields to be gathered
383  } // end if (firstApplySensitivities_)
384 
385  // Deal with the second sensitivities.
386  if (secondApplySensitivities_)
387  {
388  // Loop over the fields to be gathered.
389  for (int fieldIndex(0); fieldIndex < numFields; ++fieldIndex)
390  {
391  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldIndex];
392  int fieldNum(fieldIds_[fieldIndex]);
393  const vector<int>& elmtOffset =
394  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
395  int numBases(elmtOffset.size());
396 
397  // Gather operation for each cell in the workset.
398  for (int cell(0); cell < numCells; ++cell)
399  {
400  size_t cellLocalId(localCellIds[cell]);
401  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
402 
403  // Loop over the basis functions and fill the fields.
404  for (int basis(0); basis < numBases; ++basis)
405  {
406  int offset(elmtOffset[basis]), lid(LIDs[offset]);
407  field(cell, basis).val().fastAccessDx(0) = (*dxEvRoGed_)[lid];
408  } // end loop over the basis functions
409  } // end loop over localCellIds
410  } // end loop over the fields to be gathered
411  } // end if (secondApplySensitivities_)
412 } // end of evaluateFields() (Hessian Specialization)
413 
414 #endif // Panzer_BUILD_HESSIAN_SUPPORT
415 
416 #endif // __Panzer_GatherSolution_Epetra_Hessian_impl_hpp__
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
void setParameterList(const Teuchos::ParameterList &pl)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::string getGlobalDataKey() const
Name of the global evaluation data container to use for the source vector (all types) ...
int numFields
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
This class provides a boundary exchange communication mechanism for vectors.
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
bool secondSensitivitiesAvailable()
Are second derivative sensitivies enabled or disabled (Hessian only)
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we&#39;ll contribute, or in which we&#39;ll store, the result of computing this integral...
const std::vector< std::string > & getIndexerNames() const
Description and data layouts associated with a particular basis.
#define TEUCHOS_ASSERT(assertion_test)
std::string getSecondSensitivityDataKeyPrefix()
What prefix to use for the GEDC for second derivative sensitivity direction (Hessian only) ...
std::string getSensitivitiesName() const
The name of the sensitivities. Enables sensitivities at &quot;preEvaluate&quot; time (Jacobian and Hessian) ...