Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterDirichletResidual_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_ScatterDirichletResidual_BlockedEpetra_Hessian_impl_hpp__
44#define __Panzer_ScatterDirichletResidual_BlockedEpetra_Hessian_impl_hpp__
45
46// only do this if required by the user
47#ifdef Panzer_BUILD_HESSIAN_SUPPORT
48
49// the includes for this file come in as a result of the includes in the main
50// blocked Epetra scatter dirichlet residual file
51
52
53namespace panzer {
54
55// **************************************************************
56// Hessian Specialization
57// **************************************************************
58template<typename TRAITS,typename LO,typename GO>
60ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer<LO,int> > > & rIndexers,
61 const std::vector<Teuchos::RCP<const GlobalIndexer<LO,int> > > & cIndexers,
62 const Teuchos::ParameterList& p,
63 bool /* useDiscreteAdjoint */)
64 : rowIndexers_(rIndexers)
65 , colIndexers_(cIndexers)
66 , globalDataKey_("Residual Scatter Container")
67{
68 std::string scatterName = p.get<std::string>("Scatter Name");
69 scatterHolder_ =
70 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
71
72 // get names to be evaluated
73 const std::vector<std::string>& names =
74 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
75
76 // grab map from evaluated names to field names
77 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
78
79 Teuchos::RCP<PHX::DataLayout> dl =
80 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
81
82 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
83 local_side_id_ = p.get<int>("Local Side ID");
84
85 // build the vector of fields that this is dependent on
86 scatterFields_.resize(names.size());
87 for (std::size_t eq = 0; eq < names.size(); ++eq) {
88 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
89
90 // tell the field manager that we depend on this field
91 this->addDependentField(scatterFields_[eq]);
92 }
93
94 checkApplyBC_ = p.get<bool>("Check Apply BC");
95 if (checkApplyBC_) {
96 applyBC_.resize(names.size());
97 for (std::size_t eq = 0; eq < names.size(); ++eq) {
98 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
99 this->addDependentField(applyBC_[eq]);
100 }
101 }
102
103 // this is what this evaluator provides
104 this->addEvaluatedField(*scatterHolder_);
105
106 if (p.isType<std::string>("Global Data Key"))
107 globalDataKey_ = p.get<std::string>("Global Data Key");
108
109 if(colIndexers_.size()==0)
110 colIndexers_ = rowIndexers_;
111
112 this->setName(scatterName+" Scatter Dirichlet Residual : BlockedEpetra (Hessian)");
113}
114
115template<typename TRAITS,typename LO,typename GO>
116void
118postRegistrationSetup(typename TRAITS::SetupData /* d */,
120{
121 indexerIds_.resize(scatterFields_.size());
122 subFieldIds_.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
129 indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
130 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
131 }
132
133 // get the number of nodes (Should be renamed basis)
134 num_nodes = scatterFields_[0].extent(1);
135 num_eq = scatterFields_.size();
136}
137
138template<typename TRAITS,typename LO,typename GO>
139void
141preEvaluate(typename TRAITS::PreEvalData d)
142{
144
145 using Teuchos::rcp_dynamic_cast;
146
147 // extract dirichlet counter from container
148 Teuchos::RCP<const BLOC> blockContainer
149 = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
150
151 dirichletCounter_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
152 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
153
154 // extract linear object container
155 blockContainer = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_),true);
156 TEUCHOS_ASSERT(!Teuchos::is_null(blockContainer));
157
158 Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockContainer->get_A());
159}
160
161template<typename TRAITS,typename LO,typename GO>
162void
164evaluateFields(typename TRAITS::EvalData workset)
165{
166 using Teuchos::RCP;
167 using Teuchos::ArrayRCP;
168 using Teuchos::ptrFromRef;
169 using Teuchos::rcp_dynamic_cast;
170
171 using Thyra::SpmdVectorBase;
172
173 // for convenience pull out some objects from workset
174 std::string blockId = this->wda(workset).block_id;
175 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
176
177 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
178
179 std::vector<int> blockOffsets;
180 computeBlockOffsets(blockId,colIndexers_,blockOffsets);
181
182 std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
183
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 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
190 int rowIndexer = indexerIds_[fieldIndex];
191 int subFieldNum = subFieldIds_[fieldIndex];
192
193 // loop over each field to be scattered
194 Teuchos::ArrayRCP<double> local_dc;
195 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
196 ->getNonconstLocalData(ptrFromRef(local_dc));
197
198 auto subRowIndexer = rowIndexers_[rowIndexer];
199 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
200 const std::vector<int> & subElmtOffset = subIndicePair.first;
201 const std::vector<int> & subBasisIdMap = subIndicePair.second;
202
203 // scatter operation for each cell in workset
204 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
205 std::size_t cellLocalId = localCellIds[worksetCellIndex];
206
207 auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
208
209 // loop over basis functions
210 for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
211 int offset = subElmtOffset[basis];
212 int lid = rLIDs[offset];
213 if(lid<0) // not on this processor
214 continue;
215
216 int basisId = subBasisIdMap[basis];
217
218 if (checkApplyBC_)
219 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
220 continue;
221
222 // zero out matrix row
223 for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
224 int start = blockOffsets[colIndexer];
225 int end = blockOffsets[colIndexer+1];
226
227 if(end-start<=0)
228 continue;
229
230 // check hash table for jacobian sub block
231 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
232 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
233 // if you didn't find one before, add it to the hash table
234 if(subJac==Teuchos::null) {
235 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
236
237 // block operator is null, don't do anything (it is excluded)
238 if(Teuchos::is_null(tOp))
239 continue;
240
241 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
242 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
243 jacEpetraBlocks[blockIndex] = subJac;
244 }
245
246 int numEntries = 0;
247 int * rowIndices = 0;
248 double * rowValues = 0;
249
250 subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
251
252 for(int i=0;i<numEntries;i++)
253 rowValues[i] = 0.0;
254 }
255
256 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
257
258 local_dc[lid] = 1.0; // mark row as dirichlet
259
260 // loop over the sensitivity indices: all DOFs on a cell
261 std::vector<double> jacRow(scatterField.size(),0.0);
262
263 for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
264 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
265
266 for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
267 int start = blockOffsets[colIndexer];
268 int end = blockOffsets[colIndexer+1];
269
270 if(end-start<=0)
271 continue;
272
273 auto subColIndexer = colIndexers_[colIndexer];
274 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
275
276 TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
277
278 // check hash table for jacobian sub block
279 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
280 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
281
282 // if you didn't find one before, add it to the hash table
283 if(subJac==Teuchos::null) {
284 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
285
286 // block operator is null, don't do anything (it is excluded)
287 if(Teuchos::is_null(tOp))
288 continue;
289
290 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
291 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
292 jacEpetraBlocks[blockIndex] = subJac;
293 }
294
295 // Sum Jacobian
296 int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
297 if(err!=0) {
298 std::stringstream ss;
299 ss << "Failed inserting row: " << " (" << lid << "): ";
300 for(int i=0;i<end-start;i++)
301 ss << cLIDs[i] << " ";
302 ss << std::endl;
303 ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
304
305 ss << "scatter field = ";
306 scatterFields_[fieldIndex].print(ss);
307 ss << std::endl;
308
309 TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
310 }
311
312 }
313 }
314 }
315 }
316}
317
318}
319
320// **************************************************************
321#endif
322
323#endif
Pushes residual values into the residual vector for a Newton-based solve.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
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)