Panzer Version of the Day
Loading...
Searching...
No Matches
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"
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//
80template<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{
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//
134template<typename TRAITS, typename LO, typename GO>
135void
138 typename TRAITS::SetupData /* d */,
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//
163template<typename TRAITS, typename LO, typename GO>
164void
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;
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.
190 RCP<GlobalEvaluationData> ged;
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//
244template<typename TRAITS, typename LO, typename GO>
245void
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.
328 ArrayRCP<const double> x;
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__
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 encapsulates the needs of a gather operation to do a halo exchange for blocked vectors.
Gathers solution values from the Newton solution vector into the nodal fields of the field manager.
void evaluateFields(typename TRAITS::EvalData d)
Evaluate Fields.
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
Constructor.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Post-Registration Setup.
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.
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer > > &ugis, std::vector< int > &blockOffsets)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer > > &ugis)