Panzer  Version of the Day
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 
52 
53 #include "Phalanx_DataLayout.hpp"
54 
55 template<typename EvalT,typename 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 UniqueGlobalIndexerBase & indexerSrc,
62  const UniqueGlobalIndexerBase & 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>
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 UniqueGlobalIndexerBase & indexerSrc,
89  const UniqueGlobalIndexerBase & 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>
110 postRegistrationSetup(typename TRAITS::SetupData d,
112 {
113  // load required field numbers for fast use
114  for(std::size_t fd=0;fd<inFields_.size();++fd) {
115  // fill field data object
116  this->utils.setFieldData(inFields_[fd],fm);
117  this->utils.setFieldData(outFields_[fd],fm);
118  }
119 }
120 
121 // **********************************************************************
122 template<typename EvalT,typename TRAITS>
124 evaluateFields(typename TRAITS::EvalData workset)
125 {
126  // just copy fields if there is no AD data
127  //for(std::size_t i = 0; i < inFields_.size(); ++i)
128  // for(typename PHX::MDField<ScalarT>::size_type j = 0; j < inFields_[i].size(); ++j)
129  // outFields_[i][j] = inFields_[i][j];
130  for(std::size_t i = 0; i < inFields_.size(); ++i)
131  outFields_[i].deep_copy(inFields_[i]);
132 }
133 
134 // **********************************************************************
135 // Jacobian
136 // **********************************************************************
137 
138 template<typename TRAITS>
140 ReorderADValues_Evaluator(const std::string & outPrefix,
141  const std::vector<std::string> & inFieldNames,
142  const std::vector<Teuchos::RCP<PHX::DataLayout> > & fieldLayouts,
143  const std::string & elementBlock,
144  const UniqueGlobalIndexerBase & indexerSrc,
145  const UniqueGlobalIndexerBase & indexerDest)
146 {
147  TEUCHOS_ASSERT(inFieldNames.size()==fieldLayouts.size());
148 
149  // build the vector of fields that this is dependent on
150  for (std::size_t eq = 0; eq < inFieldNames.size(); ++eq) {
151  inFields_.push_back(PHX::MDField<const ScalarT>(inFieldNames[eq],fieldLayouts[eq]));
152  outFields_.push_back(PHX::MDField<ScalarT>(outPrefix+inFieldNames[eq],fieldLayouts[eq]));
153 
154  // tell the field manager that we depend on this field
155  this->addDependentField(inFields_[eq]);
156  this->addEvaluatedField(outFields_[eq]);
157  }
158 
159  buildSrcToDestMap(elementBlock,
160  indexerSrc,
161  indexerDest);
162 
163  this->setName(outPrefix+" Reorder AD Values");
164 }
165 
166 // **********************************************************************
167 
168 template<typename TRAITS>
170 ReorderADValues_Evaluator(const std::string & outPrefix,
171  const std::vector<std::string> & inFieldNames,
172  const std::vector<std::string> & inDOFs,
173  const std::vector<std::string> & outDOFs,
174  const std::vector<Teuchos::RCP<PHX::DataLayout> > & fieldLayouts,
175  const std::string & elementBlock,
176  const UniqueGlobalIndexerBase & indexerSrc,
177  const UniqueGlobalIndexerBase & indexerDest)
178 {
179  TEUCHOS_ASSERT(inFieldNames.size()==fieldLayouts.size());
180  TEUCHOS_ASSERT(inDOFs.size()==outDOFs.size());
181 
182  // build the vector of fields that this is dependent on
183  std::map<int,int> fieldNumberMaps;
184  for (std::size_t eq = 0; eq < inFieldNames.size(); ++eq) {
185  inFields_.push_back(PHX::MDField<const ScalarT>(inFieldNames[eq],fieldLayouts[eq]));
186  outFields_.push_back(PHX::MDField<ScalarT>(outPrefix+inFieldNames[eq],fieldLayouts[eq]));
187 
188  // tell the field manager that we depend on this field
189  this->addDependentField(inFields_[eq]);
190  this->addEvaluatedField(outFields_[eq]);
191  }
192 
193  // build a int-int map that associates fields
194  for(std::size_t i=0;i<inDOFs.size();i++) {
195  int srcFieldNum = indexerSrc.getFieldNum(inDOFs[i]);
196  int dstFieldNum = indexerDest.getFieldNum(outDOFs[i]);
197  TEUCHOS_ASSERT(srcFieldNum>=0);
198  TEUCHOS_ASSERT(dstFieldNum>=0);
199 
200  fieldNumberMaps[srcFieldNum] = dstFieldNum;
201  }
202 
203  buildSrcToDestMap(elementBlock,
204  fieldNumberMaps,
205  indexerSrc,
206  indexerDest);
207 
208  this->setName("Reorder AD Values");
209 }
210 
211 // **********************************************************************
212 template<typename TRAITS>
214 postRegistrationSetup(typename TRAITS::SetupData d,
216 {
217  // load required field numbers for fast use
218  for(std::size_t fd=0;fd<inFields_.size();++fd) {
219  // fill field data object
220  this->utils.setFieldData(inFields_[fd],fm);
221  this->utils.setFieldData(outFields_[fd],fm);
222  }
223 }
224 
225 // **********************************************************************
226 template<typename TRAITS>
228 evaluateFields(typename TRAITS::EvalData workset)
229 {
230  // for AD data do a reordering
231 
232  // 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!");
233 
234  for(std::size_t fieldIndex = 0; fieldIndex < inFields_.size(); ++fieldIndex) {
235 
236  const PHX::MDField<const ScalarT>& inField = inFields_[fieldIndex];
237  const PHX::MDField<ScalarT>& outField = outFields_[fieldIndex];
238 
239  if(inField.size()>0) {
240 
241  switch (inFields_[fieldIndex].rank()) {
242  case (1):
243  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.dimension(0); ++i) {
244  outField(i).val() = inField(i).val();
245  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
246  outField(i).fastAccessDx(dx) = inField(i).fastAccessDx(dstFromSrcMap_[dx]);
247  }
248  break;
249  case (2):
250  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.dimension(0); ++i)
251  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.dimension(1); ++j) {
252  outField(i,j).val() = inField(i,j).val();
253  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
254  outField(i,j).fastAccessDx(dx) = inField(i,j).fastAccessDx(dstFromSrcMap_[dx]);
255  }
256  break;
257  case (3):
258  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.dimension(0); ++i)
259  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.dimension(1); ++j)
260  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.dimension(2); ++k) {
261  outField(i,j,k).val() = inField(i,j,k).val();
262  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
263  outField(i,j,k).fastAccessDx(dx) = inField(i,j,k).fastAccessDx(dstFromSrcMap_[dx]);
264  }
265  break;
266  case (4):
267  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.dimension(0); ++i)
268  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.dimension(1); ++j)
269  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.dimension(2); ++k)
270  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.dimension(3); ++l) {
271  outField(i,j,k,l).val() = inField(i,j,k,l).val();
272  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
273  outField(i,j,k,l).fastAccessDx(dx) = inField(i,j,k,l).fastAccessDx(dstFromSrcMap_[dx]);
274  }
275  break;
276  case (5):
277  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.dimension(0); ++i)
278  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.dimension(1); ++j)
279  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.dimension(2); ++k)
280  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.dimension(3); ++l)
281  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.dimension(4); ++m) {
282  outField(i,j,k,l,m).val() = inField(i,j,k,l,m).val();
283  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
284  outField(i,j,k,l,m).fastAccessDx(dx) = inField(i,j,k,l,m).fastAccessDx(dstFromSrcMap_[dx]);
285  }
286  break;
287  case (6):
288  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.dimension(0); ++i)
289  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.dimension(1); ++j)
290  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.dimension(2); ++k)
291  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.dimension(3); ++l)
292  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.dimension(4); ++m)
293  for (typename PHX::MDField<ScalarT>::size_type n = 0; n < inField.dimension(5); ++n) {
294  outField(i,j,k,l,m,n).val() = inField(i,j,k,l,m,n).val();
295  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
296  outField(i,j,k,l,m,n).fastAccessDx(dx) = inField(i,j,k,l,m,n).fastAccessDx(dstFromSrcMap_[dx]);
297  }
298  break;
299  }
300 
301  }
302 
303  }
304 
305 //Irina TOFIX
306 /*
307  for(std::size_t i = 0; i < inFields_.size(); ++i) {
308 
309  for(typename PHX::MDField<ScalarT>::size_type j = 0; j < inFields_[i].size(); ++j) {
310  // allocated scalar fields
311  outFields_[i][j] = ScalarT(dstFromSrcMap_.size(), inFields_[i][j].val());
312 
313  ScalarT & outField = outFields_[i][j];
314  const ScalarT & inField = inFields_[i][j];
315 
316  // the jacobian must be initialized, otherwise its just a value copy
317  if(inField.size()>0) {
318  // loop over the sensitivity indices: all DOFs on a cell
319  outField.resize(dstFromSrcMap_.size());
320 
321  // copy jacobian entries correctly reordered
322  for(std::size_t k=0;k<dstFromSrcMap_.size();k++)
323  outField.fastAccessDx(k) = inField.fastAccessDx(dstFromSrcMap_[k]);
324  }
325 
326  outField.val() = inField.val();
327  }
328  }
329 */
330 }
331 
332 // **********************************************************************
333 template<typename TRAITS>
335 buildSrcToDestMap(const std::string & elementBlock,
336  const UniqueGlobalIndexerBase & indexerSrc,
337  const UniqueGlobalIndexerBase & indexerDest)
338 {
339  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
340  out.setOutputToRootOnly(0);
341 
342  TEUCHOS_ASSERT(indexerSrc.getComm()!=Teuchos::null);
343  TEUCHOS_ASSERT(indexerDest.getComm()!=Teuchos::null);
344 
345  const std::vector<int> & dstFieldsNum = indexerDest.getBlockFieldNumbers(elementBlock);
346 
347  // build a map between destination field numbers and source field numbers
348  std::map<int,int> fieldNumberMaps;
349  for(std::size_t i=0;i<dstFieldsNum.size();i++) {
350  std::string fieldName = indexerDest.getFieldString(dstFieldsNum[i]);
351 
352  int srcFieldNum = indexerSrc.getFieldNum(fieldName);
353  if(srcFieldNum>=0)
354  fieldNumberMaps[srcFieldNum] = dstFieldsNum[i];
355  else
356  out << "Warning: Reorder AD Values can't find field \"" << fieldName << "\"" << std::endl;
357  }
358 
359  buildSrcToDestMap(elementBlock,fieldNumberMaps,indexerSrc,indexerDest);
360 }
361 
362 // **********************************************************************
363 template<typename TRAITS>
365 buildSrcToDestMap(const std::string & elementBlock,
366  const std::map<int,int> & fieldNumberMaps,
367  const UniqueGlobalIndexerBase & indexerSrc,
368  const UniqueGlobalIndexerBase & indexerDest)
369 {
370  int maxDest = -1;
371  std::map<int,int> offsetMap; // map from source to destination offsets
372  for(std::map<int,int>::const_iterator itr=fieldNumberMaps.begin();
373  itr!=fieldNumberMaps.end();++itr) {
374  int srcField = itr->first;
375  int dstField = itr->second;
376 
377  const std::vector<int> & srcOffsets = indexerSrc.getGIDFieldOffsets(elementBlock,srcField);
378  const std::vector<int> & dstOffsets = indexerDest.getGIDFieldOffsets(elementBlock,dstField);
379 
380  // field should be the same size
381  TEUCHOS_ASSERT(srcOffsets.size()==dstOffsets.size());
382  for(std::size_t i=0;i<srcOffsets.size();i++) {
383  offsetMap[srcOffsets[i]] = dstOffsets[i];
384 
385  // provides a size for allocating an array below: we will be able
386  // to index into dstFromSrcMap_ in a simple way
387  maxDest = dstOffsets[i]>maxDest ? dstOffsets[i] : maxDest;
388  }
389  }
390 
391  // Build map
392  TEUCHOS_ASSERT(maxDest>0);
393  dstFromSrcMap_ = std::vector<int>(maxDest+1,-1);
394  for(std::map<int,int>::const_iterator itr=offsetMap.begin();
395  itr!=offsetMap.end();++itr) {
396  dstFromSrcMap_[itr->second] = itr->first;
397  }
398 }
399 
400 // **********************************************************************
401 
402 #endif
403 
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.
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 UniqueGlobalIndexerBase &indexerSrc, const UniqueGlobalIndexerBase &indexerDest)
PHX::MDField< ScalarT, Cell, BASIS > inField
PHX::MDField< ScalarT, Cell > outField
PHX::MDField< ScalarT > vector
virtual Teuchos::RCP< Teuchos::Comm< int > > getComm() const =0
std::vector< PHX::MDField< ScalarT > > outFields_
std::vector< PHX::MDField< const ScalarT > > inFields_
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
virtual const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const =0
Reorders the ad values of a specified field to match a different unique global indexer.
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
#define TEUCHOS_ASSERT(assertion_test)
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.