Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_GatherSolution_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_GatherSolution_BlockedEpetra_Hessian_impl_hpp__
44 #define __Panzer_GatherSolution_BlockedEpetra_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 // Panzer
61 #include "Panzer_PureBasis.hpp"
62 #include "Panzer_GlobalIndexer.hpp"
63 
64 // Phalanx
65 #include "Phalanx_DataLayout.hpp"
66 
67 // Teuchos
68 #include "Teuchos_Assert.hpp"
69 #include "Teuchos_FancyOStream.hpp"
70 
71 // Thyra
72 #include "Thyra_ProductVectorBase.hpp"
73 #include "Thyra_SpmdVectorBase.hpp"
74 
76 //
77 // Initializing Constructor (Hessian Specialization)
78 //
80 template<typename TRAITS, typename LO, typename GO>
83  const std::vector<Teuchos::RCP<const GlobalIndexer<LO, int>>>&
84  indexers,
85  const Teuchos::ParameterList& p)
86  :
87  indexers_(indexers)
88 {
89  using panzer::PureBasis;
90  using PHX::MDField;
91  using PHX::print;
92  using std::size_t;
93  using std::string;
94  using std::vector;
95  using Teuchos::RCP;
97  input.setParameterList(p);
98  const vector<string>& names = input.getDofNames();
99  RCP<const PureBasis> basis = input.getBasis();
100  indexerNames_ = input.getIndexerNames();
101  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
102  globalDataKey_ = input.getGlobalDataKey();
103  gatherSeedIndex_ = input.getGatherSeedIndex();
104  sensitivitiesName_ = input.getSensitivitiesName();
105  firstSensitivitiesAvailable_ = input.firstSensitivitiesAvailable();
106  secondSensitivitiesAvailable_ = input.secondSensitivitiesAvailable();
107  sensitivities2ndPrefix_ = input.getSecondSensitivityDataKeyPrefix();
108 
109  // Allocate the fields.
110  int numFields(names.size());
111  gatherFields_.resize(numFields);
112  for (int fd(0); fd < numFields; ++fd)
113  {
114  gatherFields_[fd] =
115  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
116  this->addEvaluatedField(gatherFields_[fd]);
117  } // end loop over names
118 
119  // Figure out what the first active name is.
120  string firstName("<none>"), n("GatherSolution (BlockedEpetra");
121  if (numFields > 0)
122  firstName = names[0];
123  if (not firstSensitivitiesAvailable_)
124  n += ", No First Sensitivities";
125  n += "): " + firstName + " (" + print<EvalT>() + ")";
126  this->setName(n);
127 } // end of Initializing Constructor (Hessian Specialization)
128 
130 //
131 // postRegistrationSetup() (Hessian Specialization)
132 //
134 template<typename TRAITS, typename LO, typename GO>
135 void
138  typename TRAITS::SetupData /* d */,
139  PHX::FieldManager<TRAITS>& /* fm */)
140 {
141  using std::size_t;
142  using std::string;
143  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
144  int numFields(gatherFields_.size());
145  indexerIds_.resize(numFields);
146  subFieldIds_.resize(numFields);
147  for (int fd(0); fd < numFields; ++fd)
148  {
149  // Get the field ID from the DOF manager.
150  const string& fieldName(indexerNames_[fd]);
151  indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
152  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
153  TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
154  } // end loop over gatherFields_
155  indexerNames_.clear();
156 } // end of postRegistrationSetup() (Hessian Specialization)
157 
159 //
160 // preEvaluate() (Hessian Specialization)
161 //
163 template<typename TRAITS, typename LO, typename GO>
164 void
167  typename TRAITS::PreEvalData d)
168 {
169  using std::endl;
170  using std::logic_error;
171  using std::string;
172  using Teuchos::RCP;
173  using Teuchos::rcp_dynamic_cast;
174  using Teuchos::typeName;
176  using BELOC = BlockedEpetraLinearObjContainer;
178 
179  // Manage sensitivities.
180  firstApplySensitivities_ = false;
181  if ((firstSensitivitiesAvailable_ ) and
182  (d.first_sensitivities_name == sensitivitiesName_))
183  firstApplySensitivities_ = true;
184  secondApplySensitivities_ = false;
185  if ((secondSensitivitiesAvailable_ ) and
186  (d.second_sensitivities_name == sensitivitiesName_))
187  secondApplySensitivities_ = true;
188 
189  // First try the refactored ReadOnly container.
191  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
192  if (d.gedc->containsDataObject(globalDataKey_ + post))
193  {
194  ged = d.gedc->getDataObject(globalDataKey_ + post);
195  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
196  }
197 
198  // Otherwise, try the old path.
199  else if (d.gedc->containsDataObject(globalDataKey_))
200  {
201  ged = d.gedc->getDataObject(globalDataKey_);
202 
203  // Try to extract the linear object container.
204  auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
205  auto blockedContainer = rcp_dynamic_cast<const BELOC>(ged);
206  if (not roGed.is_null())
207  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
208  else if (not blockedContainer.is_null())
209  {
210  if (useTimeDerivativeSolutionVector_)
211  x_ = rcp_dynamic_cast<ProductVectorBase<double>>
212  (blockedContainer->get_dxdt());
213  else // if (not useTimeDerivativeSolutionVector_)
214  x_ = rcp_dynamic_cast<ProductVectorBase<double>>
215  (blockedContainer->get_x());
216  } // end if roGed or blockedContainer is non-null
217  } // end if we're doing things the new or old way
218 
219  // Ensure that we actually have something.
220  TEUCHOS_ASSERT((not x_.is_null()) or (not xBvRoGed_.is_null()))
221 
222  // Don't try to extract dx if it's not required.
223  if (not secondApplySensitivities_)
224  return;
225 
226  // Now parse the second derivative direction.
227  if (d.gedc->containsDataObject(sensitivities2ndPrefix_ + globalDataKey_))
228  {
229  ged = d.gedc->getDataObject(sensitivities2ndPrefix_ + globalDataKey_);
230  dxBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
231  } // end if (d.gedc->containsDataObject(...))
232 
233  // Ensure that we actually have something.
234  TEUCHOS_TEST_FOR_EXCEPTION(dxBvRoGed_.is_null(), logic_error,
235  "Cannot find sensitivity vector associated with \"" +
236  sensitivities2ndPrefix_ + globalDataKey_ + "\" and \"" + post + "\".");
237 } // end of preEvaluate() (Hessian Specialization)
238 
240 //
241 // evaluateFields() (Hessian Specialization)
242 //
244 template<typename TRAITS, typename LO, typename GO>
245 void
248  typename TRAITS::EvalData workset)
249 {
250  using PHX::MDField;
251  using std::size_t;
252  using std::string;
253  using std::vector;
254  using Teuchos::ArrayRCP;
255  using Teuchos::ptrFromRef;
256  using Teuchos::RCP;
257  using Teuchos::rcp_dynamic_cast;
259  using Thyra::SpmdVectorBase;
260  using Thyra::VectorBase;
261 
262  // For convenience, pull out some objects from the workset.
263  string blockId(this->wda(workset).block_id);
264  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
265  int numFields(gatherFields_.size()), numCells(localCellIds.size());
266  double seedValue(0);
267  if (firstApplySensitivities_)
268  {
269  if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
270  seedValue = workset.alpha;
271  else if (gatherSeedIndex_ < 0)
272  seedValue = workset.beta;
273  else if (not useTimeDerivativeSolutionVector_)
274  seedValue = workset.gather_seeds[gatherSeedIndex_];
275  else
276  TEUCHOS_ASSERT(false);
277  } // end if (firstApplySensitivities_)
278 
279  // Turn off sensitivies. This may be faster if we don't expand the term, but
280  // I suspect not, because anywhere it is used the full complement of
281  // sensitivies will be needed anyway.
282  if (not firstApplySensitivities_)
283  seedValue = 0.0;
284 
285  vector<int> blockOffsets;
286  computeBlockOffsets(blockId, indexers_, blockOffsets);
287  if (x_.is_null())
288  {
289  // Loop over the fields to be gathered.
290  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
291  {
292  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
293  int indexerId(indexerIds_[fieldInd]),
294  subFieldNum(subFieldIds_[fieldInd]);
295 
296  // Grab the local data for inputing.
297  auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
298  auto subRowIndexer = indexers_[indexerId];
299  const vector<int>& elmtOffset =
300  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
301  int numBases(elmtOffset.size());
302 
303  // Gather operation for each cell in the workset.
304  for (int cell(0); cell < numCells; ++cell)
305  {
306  size_t cellLocalId(localCellIds[cell]);
307  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
308 
309  // Loop over the basis functions and fill the fields.
310  for (int basis(0); basis < numBases; ++basis)
311  {
312  int offset(elmtOffset[basis]), lid(LIDs[offset]);
313  field(cell, basis) = (*xEvRoGed)[lid];
314  } // end loop over the basis functions
315  } // end loop over localCellIds
316  } // end loop over the fields to be gathered
317  }
318  else // if (not x_.is_null())
319  {
320  // Loop over the fields to be gathered.
321  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
322  {
323  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
324  int indexerId(indexerIds_[fieldInd]),
325  subFieldNum(subFieldIds_[fieldInd]);
326 
327  // Grab the local data for inputing.
329  rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
330  getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
331  auto subRowIndexer = indexers_[indexerId];
332  const vector<int>& elmtOffset =
333  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
334  int numBases(elmtOffset.size());
335 
336  // Gather operation for each cell in the workset.
337  for (int cell(0); cell < numCells; ++cell)
338  {
339  size_t cellLocalId(localCellIds[cell]);
340  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
341 
342  // Loop over the basis functions and fill the fields.
343  for (int basis(0); basis < numBases; ++basis)
344  {
345  int offset(elmtOffset[basis]), lid(LIDs[offset]);
346  field(cell, basis) = x[lid];
347  } // end loop over the basis functions
348  } // end loop over localCellIds
349  } // end loop over the fields to be gathered
350  } // end if (x_.is_null()) or not
351 
352  // Deal with the first sensitivities.
353  if (firstApplySensitivities_)
354  {
355  // Loop over the fields to be gathered.
356  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
357  {
358  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
359  int indexerId(indexerIds_[fieldInd]),
360  subFieldNum(subFieldIds_[fieldInd]);
361  auto subRowIndexer = indexers_[indexerId];
362  const vector<int>& elmtOffset =
363  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
364  int startBlkOffset(blockOffsets[indexerId]),
365  numBases(elmtOffset.size());
366 
367  // Gather operation for each cell in the workset.
368  for (int cell(0); cell < numCells; ++cell)
369  {
370  // Loop over the basis functions and fill the fields.
371  for (int basis(0); basis < numBases; ++basis)
372  {
373  int offset(elmtOffset[basis]);
374  field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
375  } // end loop over the basis functions
376  } // end loop over localCellIds
377  } // end loop over the fields to be gathered
378  } // end if (firstApplySensitivities_)
379 
380  // Deal with the second sensitivities.
381  if (secondApplySensitivities_)
382  {
383  // Loop over the fields to be gathered.
384  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
385  {
386  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
387  int indexerId(indexerIds_[fieldInd]),
388  subFieldNum(subFieldIds_[fieldInd]);
389 
390  // Grab the local data for inputing.
391  auto dxEvRoGed = dxBvRoGed_->getGEDBlock(indexerId);
392  auto subRowIndexer = indexers_[indexerId];
393  const vector<int>& elmtOffset =
394  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
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 = subRowIndexer->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_BlockedEpetra_Hessian_impl_hpp__
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
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) ...
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer >> &ugis, std::vector< int > &blockOffsets)
int numFields
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer >> &ugis)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Post-Registration Setup.
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
Constructor.
bool secondSensitivitiesAvailable()
Are second derivative sensitivies enabled or disabled (Hessian only)
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
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
void evaluateFields(typename TRAITS::EvalData d)
Evaluate Fields.
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) ...
std::string typeName(const T &t)