Panzer Version of the Day
Loading...
Searching...
No Matches
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
74template<typename TRAITS,typename LO,typename GO>
76ScatterResidual_Epetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
77 const Teuchos::RCP<const panzer::GlobalIndexer> & 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
96 Teuchos::RCP<PHX::DataLayout> dl =
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// **********************************************************************
124template<typename TRAITS,typename LO,typename GO>
126postRegistrationSetup(typename TRAITS::SetupData /* d */,
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// **********************************************************************
139template<typename TRAITS,typename LO,typename GO>
141preEvaluate(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// **********************************************************************
154template<typename TRAITS,typename LO,typename GO>
156evaluateFields(typename TRAITS::EvalData workset)
157{
158 using PHX::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
173 const Teuchos::RCP<const panzer::GlobalIndexer>&
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())),
220 TEUCHOS_ASSERT_EQUALITY(err,0);
221 }
222 } // end rowBasisNum
223 } // end fieldIndex
224 }
225}
226
227// **********************************************************************
228
229#endif
230
231#endif
Pushes residual values into the residual vector for a Newton-based solve.
T * ptrFromStlVector(std::vector< T > &v)