Panzer Version of the Day
Loading...
Searching...
No Matches
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//
82template<typename TRAITS, typename LO, typename GO>
85 const Teuchos::RCP<const panzer::GlobalIndexer>& indexer,
86 const Teuchos::ParameterList& p)
87 :
88 globalIndexer_(indexer)
89{
91 using PHX::MDField;
92 using PHX::print;
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 + " (" + print<EvalT>() + ") ";
127 this->setName(n);
128} // end of Initializing Constructor (Hessian Specialization)
129
131//
132// postRegistrationSetup() (Hessian Specialization)
133//
135template<typename TRAITS, typename LO, typename GO>
136void
139 typename TRAITS::SetupData /* d */,
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//
167template<typename TRAITS, typename LO, typename GO>
168void
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;
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//
256template<typename TRAITS, typename LO, typename GO>
257void
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__
int numFields
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
This class provides a boundary exchange communication mechanism for vectors.
Gathers solution values from the Newton solution vector into the nodal fields of the field manager.
std::string getGlobalDataKey() const
Name of the global evaluation data container to use for the source vector (all types)
void setParameterList(const Teuchos::ParameterList &pl)
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
std::string getSecondSensitivityDataKeyPrefix()
What prefix to use for the GEDC for second derivative sensitivity direction (Hessian only)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
const std::vector< std::string > & getIndexerNames() const
std::string getSensitivitiesName() const
The name of the sensitivities. Enables sensitivities at "preEvaluate" time (Jacobian and Hessian)
bool secondSensitivitiesAvailable()
Are second derivative sensitivies enabled or disabled (Hessian only)
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
Description and data layouts associated with a particular basis.