Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterDirichletResidual_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_ScatterDirichletResidual_Epetra_Hessian_impl_hpp__
44#define __Panzer_ScatterDirichletResidual_Epetra_Hessian_impl_hpp__
45
46// only do this if required by the user
47#ifdef Panzer_BUILD_HESSIAN_SUPPORT
48
50
51// the includes for this file come in as a result of the includes in the main
52// Epetra scatter dirichlet residual file
53
54namespace panzer {
55
56// **************************************************************
57// Hessian Specialization
58// **************************************************************
59template<typename TRAITS,typename LO,typename GO>
61ScatterDirichletResidual_Epetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
62 const Teuchos::RCP<const panzer::GlobalIndexer> & cIndexer,
63 const Teuchos::ParameterList& p)
64 : globalIndexer_(indexer)
65 , colGlobalIndexer_(cIndexer)
66 , globalDataKey_("Residual Scatter Container")
67 , preserveDiagonal_(false)
68{
69 std::string scatterName = p.get<std::string>("Scatter Name");
70 scatterHolder_ =
71 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
72
73 // get names to be evaluated
74 const std::vector<std::string>& names =
75 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
76
77 // grab map from evaluated names to field names
78 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
79
80 Teuchos::RCP<PHX::DataLayout> dl =
81 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
82
83 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
84 local_side_id_ = p.get<int>("Local Side ID");
85
86 // build the vector of fields that this is dependent on
87 scatterFields_.resize(names.size());
88 for (std::size_t eq = 0; eq < names.size(); ++eq) {
89 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
90
91 // tell the field manager that we depend on this field
92 this->addDependentField(scatterFields_[eq]);
93 }
94
95 checkApplyBC_ = p.get<bool>("Check Apply BC");
96 if (checkApplyBC_) {
97 applyBC_.resize(names.size());
98 for (std::size_t eq = 0; eq < names.size(); ++eq) {
99 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
100 this->addDependentField(applyBC_[eq]);
101 }
102 }
103
104 // this is what this evaluator provides
105 this->addEvaluatedField(*scatterHolder_);
106
107 if (p.isType<std::string>("Global Data Key"))
108 globalDataKey_ = p.get<std::string>("Global Data Key");
109
110 if (p.isType<bool>("Preserve Diagonal"))
111 preserveDiagonal_ = p.get<bool>("Preserve Diagonal");
112
113 this->setName(scatterName+" Scatter Residual (Hessian)");
114}
115
116template<typename TRAITS,typename LO,typename GO>
117void
119postRegistrationSetup(typename TRAITS::SetupData /* d */,
121{
122 fieldIds_.resize(scatterFields_.size());
123
124 // load required field numbers for fast use
125 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
126 // get field ID from DOF manager
127 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
128 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
129 }
130
131 // get the number of nodes (Should be renamed basis)
132 num_nodes = scatterFields_[0].extent(1);
133 num_eq = scatterFields_.size();
134}
135
136template<typename TRAITS,typename LO,typename GO>
137void
139preEvaluate(typename TRAITS::PreEvalData d)
140{
141 // extract linear object container
142 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
143
144 if(epetraContainer_==Teuchos::null) {
145 // extract linear object container
146 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
147 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc,true);
148
149 dirichletCounter_ = Teuchos::null;
150 }
151 else {
152 // extract dirichlet counter from container
153 Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc->getDataObject("Dirichlet Counter");
154 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(dataContainer,true);
155
156 dirichletCounter_ = epetraContainer->get_f();
157 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
158 }
159}
160
161template<typename TRAITS,typename LO,typename GO>
162void
164evaluateFields(typename TRAITS::EvalData workset)
165{
167 using std::vector;
168
169 // TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
170 // "ScatterDirichletResidual_Epetra<Hessian> is not yet implemented"); // just in case
171
172 PHX::View<const int*, Kokkos::LayoutRight, PHX::Device> cLIDs, rLIDs;
173 std::vector<double> jacRow;
174
175 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
176
177 // for convenience pull out some objects from workset
178 std::string blockId = this->wda(workset).block_id;
179 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
180
181 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
182 TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
183 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
184 // NOTE: A reordering of these loops will likely improve performance
185 // The "getGIDFieldOffsets may be expensive. However the
186 // "getElementGIDs" can be cheaper. However the lookup for LIDs
187 // may be more expensive!
188
189 // scatter operation for each cell in workset
190 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
191 std::size_t cellLocalId = localCellIds[worksetCellIndex];
192
193 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
194 if(useColumnIndexer)
195 cLIDs = colGlobalIndexer_->getElementLIDs(cellLocalId);
196 else
197 cLIDs = rLIDs;
198
199 // loop over each field to be scattered
200 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
201 int fieldNum = fieldIds_[fieldIndex];
202
203 // this call "should" get the right ordering according to the Intrepid2 basis
204 const std::pair<std::vector<int>,std::vector<int> > & indicePair
205 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
206 const std::vector<int> & elmtOffset = indicePair.first;
207 const std::vector<int> & basisIdMap = indicePair.second;
208
209 // loop over basis functions
210 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
211 int offset = elmtOffset[basis];
212 int row = rLIDs[offset];
213 if(row<0) // not on this processor
214 continue;
215
216 int basisId = basisIdMap[basis];
217
218 if (checkApplyBC_)
219 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
220 continue;
221
222 // zero out matrix row
223 {
224 int numEntries = 0;
225 int * rowIndices = 0;
226 double * rowValues = 0;
227
228 Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
229
230 for(int i=0;i<numEntries;i++) {
231 if(preserveDiagonal_) {
232 if(row!=rowIndices[i])
233 rowValues[i] = 0.0;
234 }
235 else
236 rowValues[i] = 0.0;
237 }
238 }
239
240 // int gid = GIDs[offset];
241 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
242
243 if(dirichletCounter_!=Teuchos::null) {
244 // std::cout << "Writing " << row << " " << dirichletCounter_->MyLength() << std::endl;
245 (*dirichletCounter_)[row] = 1.0; // mark row as dirichlet
246 }
247
248 // loop over the sensitivity indices: all DOFs on a cell
249 jacRow.resize(scatterField.size());
250 for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
251 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
252
253 if(!preserveDiagonal_) {
254 int err = Jac->ReplaceMyValues(row, cLIDs.size(),
255 ptrFromStlVector(jacRow), &cLIDs[0]);
256 TEUCHOS_ASSERT(err==0);
257 }
258 }
259 }
260 }
261}
262
263}
264
265// **************************************************************
266#endif
267
268#endif
Pushes residual values into the residual vector for a Newton-based solve.
T * ptrFromStlVector(std::vector< T > &v)