Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_ReorderADValues_Evaluator_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_ReorderADValues_Evaluator_impl_hpp__
44 #define __Panzer_ReorderADValues_Evaluator_impl_hpp__
45 
46 
47 #include "Teuchos_RCP.hpp"
48 #include "Teuchos_Assert.hpp"
49 #include "Teuchos_FancyOStream.hpp"
50 
51 #include "Panzer_GlobalIndexer.hpp"
52 
53 #include "Phalanx_DataLayout.hpp"
54 
55 template<typename EvalT,typename TRAITS>
56 panzer::ReorderADValues_Evaluator<EvalT, TRAITS>::
57 ReorderADValues_Evaluator(const std::string & outPrefix,
58  const std::vector<std::string> & inFieldNames,
59  const std::vector<Teuchos::RCP<PHX::DataLayout> > & fieldLayouts,
60  const std::string & /* elementBlock */,
61  const GlobalIndexer & /* indexerSrc */,
62  const GlobalIndexer & /* indexerDest */)
63 {
64  TEUCHOS_ASSERT(inFieldNames.size()==fieldLayouts.size());
65 
66  // build the vector of fields that this is dependent on
67  for (std::size_t eq = 0; eq < inFieldNames.size(); ++eq) {
68  inFields_.push_back(PHX::MDField<const ScalarT>(inFieldNames[eq],fieldLayouts[eq]));
69  outFields_.push_back(PHX::MDField<ScalarT>(outPrefix+inFieldNames[eq],fieldLayouts[eq]));
70 
71  // tell the field manager that we depend on this field
72  this->addDependentField(inFields_[eq]);
73  this->addEvaluatedField(outFields_[eq]);
74  }
75 
76  this->setName(outPrefix+" Reorder AD Values");
77 }
78 
79 // **********************************************************************
80 template<typename EvalT,typename TRAITS>
81 panzer::ReorderADValues_Evaluator<EvalT, TRAITS>::
82 ReorderADValues_Evaluator(const std::string & outPrefix,
83  const std::vector<std::string> & inFieldNames,
84  const std::vector<std::string> & inDOFs,
85  const std::vector<std::string> & outDOFs,
86  const std::vector<Teuchos::RCP<PHX::DataLayout> > & fieldLayouts,
87  const std::string & /* elementBlock */,
88  const GlobalIndexer & /* indexerSrc */,
89  const GlobalIndexer & /* indexerDest */)
90 {
91  TEUCHOS_ASSERT(inFieldNames.size()==fieldLayouts.size());
92  TEUCHOS_ASSERT(inDOFs.size()==outDOFs.size());
93 
94  // build the vector of fields that this is dependent on
95  for (std::size_t eq = 0; eq < inFieldNames.size(); ++eq) {
96  inFields_.push_back(PHX::MDField<const ScalarT>(inFieldNames[eq],fieldLayouts[eq]));
97  outFields_.push_back(PHX::MDField<ScalarT>(outPrefix+inFieldNames[eq],fieldLayouts[eq]));
98 
99  // tell the field manager that we depend on this field
100  this->addDependentField(inFields_[eq]);
101  this->addEvaluatedField(outFields_[eq]);
102  }
103 
104  this->setName("Reorder AD Values");
105 }
106 
107 // **********************************************************************
108 template<typename EvalT,typename TRAITS>
109 void panzer::ReorderADValues_Evaluator<EvalT, TRAITS>::
110 evaluateFields(typename TRAITS::EvalData /* workset */)
111 {
112  // just copy fields if there is no AD data
113  //for(std::size_t i = 0; i < inFields_.size(); ++i)
114  // for(typename PHX::MDField<ScalarT>::size_type j = 0; j < inFields_[i].size(); ++j)
115  // outFields_[i][j] = inFields_[i][j];
116  for(std::size_t i = 0; i < inFields_.size(); ++i)
117  outFields_[i].deep_copy(inFields_[i]);
118 }
119 
120 // **********************************************************************
121 // Jacobian
122 // **********************************************************************
123 
124 template<typename TRAITS>
125 panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
126 ReorderADValues_Evaluator(const std::string & outPrefix,
127  const std::vector<std::string> & inFieldNames,
128  const std::vector<Teuchos::RCP<PHX::DataLayout> > & fieldLayouts,
129  const std::string & elementBlock,
130  const GlobalIndexer & indexerSrc,
131  const GlobalIndexer & indexerDest)
132 {
133  TEUCHOS_ASSERT(inFieldNames.size()==fieldLayouts.size());
134 
135  // build the vector of fields that this is dependent on
136  inFields_.resize(inFieldNames.size());
137  outFields_.resize(inFieldNames.size());
138  for (std::size_t eq = 0; eq < inFieldNames.size(); ++eq) {
139  inFields_[eq] = PHX::MDField<const ScalarT>(inFieldNames[eq],fieldLayouts[eq]);
140  outFields_[eq] = PHX::MDField<ScalarT>(outPrefix+inFieldNames[eq],fieldLayouts[eq]);
141 
142  // tell the field manager that we depend on this field
143  this->addDependentField(inFields_[eq]);
144  this->addEvaluatedField(outFields_[eq]);
145  }
146 
147  buildSrcToDestMap(elementBlock,
148  indexerSrc,
149  indexerDest);
150 
151  this->setName(outPrefix+" Reorder AD Values");
152 }
153 
154 // **********************************************************************
155 
156 template<typename TRAITS>
157 panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
158 ReorderADValues_Evaluator(const std::string & outPrefix,
159  const std::vector<std::string> & inFieldNames,
160  const std::vector<std::string> & inDOFs,
161  const std::vector<std::string> & outDOFs,
162  const std::vector<Teuchos::RCP<PHX::DataLayout> > & fieldLayouts,
163  const std::string & elementBlock,
164  const GlobalIndexer & indexerSrc,
165  const GlobalIndexer & indexerDest)
166 {
167  TEUCHOS_ASSERT(inFieldNames.size()==fieldLayouts.size());
168  TEUCHOS_ASSERT(inDOFs.size()==outDOFs.size());
169 
170  // build the vector of fields that this is dependent on
171  std::map<int,int> fieldNumberMaps;
172  for (std::size_t eq = 0; eq < inFieldNames.size(); ++eq) {
173  inFields_.push_back(PHX::MDField<const ScalarT>(inFieldNames[eq],fieldLayouts[eq]));
174  outFields_.push_back(PHX::MDField<ScalarT>(outPrefix+inFieldNames[eq],fieldLayouts[eq]));
175 
176  // tell the field manager that we depend on this field
177  this->addDependentField(inFields_[eq]);
178  this->addEvaluatedField(outFields_[eq]);
179  // Don't share so we can avoid zeroing out off blck Jacobian entries
180  this->addUnsharedField(outFields_[eq].fieldTag().clone());
181  }
182 
183  // build a int-int map that associates fields
184  for(std::size_t i=0;i<inDOFs.size();i++) {
185  int srcFieldNum = indexerSrc.getFieldNum(inDOFs[i]);
186  int dstFieldNum = indexerDest.getFieldNum(outDOFs[i]);
187  TEUCHOS_ASSERT(srcFieldNum>=0);
188  TEUCHOS_ASSERT(dstFieldNum>=0);
189 
190  fieldNumberMaps[srcFieldNum] = dstFieldNum;
191  }
192 
193  buildSrcToDestMap(elementBlock,
194  fieldNumberMaps,
195  indexerSrc,
196  indexerDest);
197 
198  this->setName("Reorder AD Values");
199 }
200 
201 // **********************************************************************
202 template<typename TRAITS>
203 void panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
204 evaluateFields(typename TRAITS::EvalData /* workset */)
205 {
206  // for AD data do a reordering
207 
208  // TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"ERROR: panzer::ReorderADValues_Evaluator: This is currently broken for the Kokkkos Transition! Contact Drekar team to fix!");
209 
210  for(std::size_t fieldIndex = 0; fieldIndex < inFields_.size(); ++fieldIndex) {
211 
212  const PHX::MDField<const ScalarT>& inField = inFields_[fieldIndex];
213  const PHX::MDField<ScalarT>& outField = outFields_[fieldIndex];
214 
215  if(inField.size()>0) {
216 
217  switch (inFields_[fieldIndex].rank()) {
218  case (1):
219  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i) {
220  outField(i).val() = inField(i).val();
221  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
222  outField(i).fastAccessDx(dx) = inField(i).fastAccessDx(dstFromSrcMap_[dx]);
223  }
224  break;
225  case (2):
226  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
227  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j) {
228  outField(i,j).val() = inField(i,j).val();
229  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
230  outField(i,j).fastAccessDx(dx) = inField(i,j).fastAccessDx(dstFromSrcMap_[dx]);
231  }
232  break;
233  case (3):
234  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
235  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
236  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k) {
237  outField(i,j,k).val() = inField(i,j,k).val();
238  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
239  outField(i,j,k).fastAccessDx(dx) = inField(i,j,k).fastAccessDx(dstFromSrcMap_[dx]);
240  }
241  break;
242  case (4):
243  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
244  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
245  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
246  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l) {
247  outField(i,j,k,l).val() = inField(i,j,k,l).val();
248  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
249  outField(i,j,k,l).fastAccessDx(dx) = inField(i,j,k,l).fastAccessDx(dstFromSrcMap_[dx]);
250  }
251  break;
252  case (5):
253  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
254  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
255  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
256  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l)
257  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.extent(4); ++m) {
258  outField(i,j,k,l,m).val() = inField(i,j,k,l,m).val();
259  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
260  outField(i,j,k,l,m).fastAccessDx(dx) = inField(i,j,k,l,m).fastAccessDx(dstFromSrcMap_[dx]);
261  }
262  break;
263  case (6):
264  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
265  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
266  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
267  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l)
268  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.extent(4); ++m)
269  for (typename PHX::MDField<ScalarT>::size_type n = 0; n < inField.extent(5); ++n) {
270  outField(i,j,k,l,m,n).val() = inField(i,j,k,l,m,n).val();
271  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
272  outField(i,j,k,l,m,n).fastAccessDx(dx) = inField(i,j,k,l,m,n).fastAccessDx(dstFromSrcMap_[dx]);
273  }
274  break;
275  }
276 
277  }
278 
279  }
280 
281 //Irina TOFIX
282 /*
283  for(std::size_t i = 0; i < inFields_.size(); ++i) {
284 
285  for(typename PHX::MDField<ScalarT>::size_type j = 0; j < inFields_[i].size(); ++j) {
286  // allocated scalar fields
287  outFields_[i][j] = ScalarT(dstFromSrcMap_.size(), inFields_[i][j].val());
288 
289  ScalarT & outField = outFields_[i][j];
290  const ScalarT & inField = inFields_[i][j];
291 
292  // the jacobian must be initialized, otherwise its just a value copy
293  if(inField.size()>0) {
294  // loop over the sensitivity indices: all DOFs on a cell
295  outField.resize(dstFromSrcMap_.size());
296 
297  // copy jacobian entries correctly reordered
298  for(std::size_t k=0;k<dstFromSrcMap_.size();k++)
299  outField.fastAccessDx(k) = inField.fastAccessDx(dstFromSrcMap_[k]);
300  }
301 
302  outField.val() = inField.val();
303  }
304  }
305 */
306 }
307 
308 // **********************************************************************
309 template<typename TRAITS>
310 void panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
311 buildSrcToDestMap(const std::string & elementBlock,
312  const GlobalIndexer & indexerSrc,
313  const GlobalIndexer & indexerDest)
314 {
315  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
316  out.setOutputToRootOnly(0);
317 
318  TEUCHOS_ASSERT(indexerSrc.getComm()!=Teuchos::null);
319  TEUCHOS_ASSERT(indexerDest.getComm()!=Teuchos::null);
320 
321  const std::vector<int> & dstFieldsNum = indexerDest.getBlockFieldNumbers(elementBlock);
322 
323  // build a map between destination field numbers and source field numbers
324  std::map<int,int> fieldNumberMaps;
325  for(std::size_t i=0;i<dstFieldsNum.size();i++) {
326  std::string fieldName = indexerDest.getFieldString(dstFieldsNum[i]);
327 
328  int srcFieldNum = indexerSrc.getFieldNum(fieldName);
329  if(srcFieldNum>=0)
330  fieldNumberMaps[srcFieldNum] = dstFieldsNum[i];
331  else
332  out << "Warning: Reorder AD Values can't find field \"" << fieldName << "\"" << std::endl;
333  }
334 
335  buildSrcToDestMap(elementBlock,fieldNumberMaps,indexerSrc,indexerDest);
336 }
337 
338 // **********************************************************************
339 template<typename TRAITS>
340 void panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
341 buildSrcToDestMap(const std::string & elementBlock,
342  const std::map<int,int> & fieldNumberMaps,
343  const GlobalIndexer & indexerSrc,
344  const GlobalIndexer & indexerDest)
345 {
346  int maxDest = -1;
347  std::map<int,int> offsetMap; // map from source to destination offsets
348  for(std::map<int,int>::const_iterator itr=fieldNumberMaps.begin();
349  itr!=fieldNumberMaps.end();++itr) {
350  int srcField = itr->first;
351  int dstField = itr->second;
352 
353  const std::vector<int> & srcOffsets = indexerSrc.getGIDFieldOffsets(elementBlock,srcField);
354  const std::vector<int> & dstOffsets = indexerDest.getGIDFieldOffsets(elementBlock,dstField);
355 
356  // field should be the same size
357  TEUCHOS_ASSERT(srcOffsets.size()==dstOffsets.size());
358  for(std::size_t i=0;i<srcOffsets.size();i++) {
359  offsetMap[srcOffsets[i]] = dstOffsets[i];
360 
361  // provides a size for allocating an array below: we will be able
362  // to index into dstFromSrcMap_ in a simple way
363  maxDest = dstOffsets[i]>maxDest ? dstOffsets[i] : maxDest;
364  }
365  }
366 
367  // Build map
368  TEUCHOS_ASSERT(maxDest>0);
369  dstFromSrcMap_ = std::vector<int>(maxDest+1,-1);
370  for(std::map<int,int>::const_iterator itr=offsetMap.begin();
371  itr!=offsetMap.end();++itr) {
372  dstFromSrcMap_[itr->second] = itr->first;
373  }
374 }
375 
376 // **********************************************************************
377 
378 #endif
379 
panzer::GlobalIndexer::getFieldString
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
panzer::GlobalIndexer::getComm
virtual Teuchos::RCP< Teuchos::Comm< int > > getComm() const =0
panzer::ReorderADValues_Evaluator::outFields_
std::vector< PHX::MDField< ScalarT > > outFields_
Definition: Panzer_ReorderADValues_Evaluator_decl.hpp:99
panzer::ReorderADValues_Evaluator::evaluateFields
void evaluateFields(typename TRAITS::EvalData d)
Definition: Panzer_ReorderADValues_Evaluator_impl.hpp:110
Teuchos::basic_FancyOStream
panzer::ReorderADValues_Evaluator::inFields_
std::vector< PHX::MDField< const ScalarT > > inFields_
Definition: Panzer_ReorderADValues_Evaluator_decl.hpp:96
panzer::GlobalIndexer::getFieldNum
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
Teuchos::basic_FancyOStream::setOutputToRootOnly
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
panzer::ReorderADValues_Evaluator
Reorders the ad values of a specified field to match a different unique global indexer.
Definition: Panzer_ReorderADValues_Evaluator_decl.hpp:70
Panzer_GlobalIndexer.hpp
panzer::GlobalIndexer::getGIDFieldOffsets
virtual const std::vector< int > & getGIDFieldOffsets(const std::string &blockId, int fieldNum) const =0
Use the field pattern so that you can find a particular field in the GIDs array.
panzer::GlobalIndexer::getBlockFieldNumbers
virtual const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const =0
panzer::ReorderADValues_Evaluator::ReorderADValues_Evaluator
ReorderADValues_Evaluator(const std::string &outPrefix, const std::vector< std::string > &inFieldNames, const std::vector< Teuchos::RCP< PHX::DataLayout > > &fieldLayouts, const std::string &elementBlock, const GlobalIndexer &indexerSrc, const GlobalIndexer &indexerDest)
Definition: Panzer_ReorderADValues_Evaluator_impl.hpp:57
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Teuchos_RCP.hpp
Generated by   doxygen 1.8.5