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 
51 #include "Panzer_GlobalIndexer.hpp"
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 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>
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>
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>
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>
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>
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 
213  const auto & inField_v = inFields_[fieldIndex].get_view();
214  const auto & outField_v = outFields_[fieldIndex].get_view();
215  auto inField = Kokkos::create_mirror_view(inField_v);
216  auto outField = Kokkos::create_mirror_view(outField_v);
217  Kokkos::deep_copy(inField, inField_v);
218 
219  if(inField.size()>0) {
220 
221  switch (inFields_[fieldIndex].rank()) {
222  case (1):
223  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i) {
224  outField(i).val() = inField(i).val();
225  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
226  outField(i).fastAccessDx(dx) = inField(i).fastAccessDx(dstFromSrcMap_[dx]);
227  }
228  break;
229  case (2):
230  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
231  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j) {
232  outField(i,j).val() = inField(i,j).val();
233  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
234  outField(i,j).fastAccessDx(dx) = inField(i,j).fastAccessDx(dstFromSrcMap_[dx]);
235  }
236  break;
237  case (3):
238  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
239  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
240  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k) {
241  outField(i,j,k).val() = inField(i,j,k).val();
242  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
243  outField(i,j,k).fastAccessDx(dx) = inField(i,j,k).fastAccessDx(dstFromSrcMap_[dx]);
244  }
245  break;
246  case (4):
247  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
248  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
249  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
250  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l) {
251  outField(i,j,k,l).val() = inField(i,j,k,l).val();
252  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
253  outField(i,j,k,l).fastAccessDx(dx) = inField(i,j,k,l).fastAccessDx(dstFromSrcMap_[dx]);
254  }
255  break;
256  case (5):
257  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
258  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
259  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
260  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l)
261  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.extent(4); ++m) {
262  outField(i,j,k,l,m).val() = inField(i,j,k,l,m).val();
263  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
264  outField(i,j,k,l,m).fastAccessDx(dx) = inField(i,j,k,l,m).fastAccessDx(dstFromSrcMap_[dx]);
265  }
266  break;
267  case (6):
268  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
269  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
270  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
271  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l)
272  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.extent(4); ++m)
273  for (typename PHX::MDField<ScalarT>::size_type n = 0; n < inField.extent(5); ++n) {
274  outField(i,j,k,l,m,n).val() = inField(i,j,k,l,m,n).val();
275  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
276  outField(i,j,k,l,m,n).fastAccessDx(dx) = inField(i,j,k,l,m,n).fastAccessDx(dstFromSrcMap_[dx]);
277  }
278  break;
279  }
280 
281  }
282  Kokkos::deep_copy(outField_v, outField);
283  }
284 
285 //Irina TOFIX
286 /*
287  for(std::size_t i = 0; i < inFields_.size(); ++i) {
288 
289  for(typename PHX::MDField<ScalarT>::size_type j = 0; j < inFields_[i].size(); ++j) {
290  // allocated scalar fields
291  outFields_[i][j] = ScalarT(dstFromSrcMap_.size(), inFields_[i][j].val());
292 
293  ScalarT & outField = outFields_[i][j];
294  const ScalarT & inField = inFields_[i][j];
295 
296  // the jacobian must be initialized, otherwise its just a value copy
297  if(inField.size()>0) {
298  // loop over the sensitivity indices: all DOFs on a cell
299  outField.resize(dstFromSrcMap_.size());
300 
301  // copy jacobian entries correctly reordered
302  for(std::size_t k=0;k<dstFromSrcMap_.size();k++)
303  outField.fastAccessDx(k) = inField.fastAccessDx(dstFromSrcMap_[k]);
304  }
305 
306  outField.val() = inField.val();
307  }
308  }
309 */
310 }
311 
312 // **********************************************************************
313 template<typename TRAITS>
315 buildSrcToDestMap(const std::string & elementBlock,
316  const GlobalIndexer & indexerSrc,
317  const GlobalIndexer & indexerDest)
318 {
319  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
320  out.setOutputToRootOnly(0);
321 
322  TEUCHOS_ASSERT(indexerSrc.getComm()!=Teuchos::null);
323  TEUCHOS_ASSERT(indexerDest.getComm()!=Teuchos::null);
324 
325  const std::vector<int> & dstFieldsNum = indexerDest.getBlockFieldNumbers(elementBlock);
326 
327  // build a map between destination field numbers and source field numbers
328  std::map<int,int> fieldNumberMaps;
329  for(std::size_t i=0;i<dstFieldsNum.size();i++) {
330  std::string fieldName = indexerDest.getFieldString(dstFieldsNum[i]);
331 
332  int srcFieldNum = indexerSrc.getFieldNum(fieldName);
333  if(srcFieldNum>=0)
334  fieldNumberMaps[srcFieldNum] = dstFieldsNum[i];
335  else
336  out << "Warning: Reorder AD Values can't find field \"" << fieldName << "\"" << std::endl;
337  }
338 
339  buildSrcToDestMap(elementBlock,fieldNumberMaps,indexerSrc,indexerDest);
340 }
341 
342 // **********************************************************************
343 template<typename TRAITS>
345 buildSrcToDestMap(const std::string & elementBlock,
346  const std::map<int,int> & fieldNumberMaps,
347  const GlobalIndexer & indexerSrc,
348  const GlobalIndexer & indexerDest)
349 {
350  int maxDest = -1;
351  std::map<int,int> offsetMap; // map from source to destination offsets
352  for(std::map<int,int>::const_iterator itr=fieldNumberMaps.begin();
353  itr!=fieldNumberMaps.end();++itr) {
354  int srcField = itr->first;
355  int dstField = itr->second;
356 
357  const std::vector<int> & srcOffsets = indexerSrc.getGIDFieldOffsets(elementBlock,srcField);
358  const std::vector<int> & dstOffsets = indexerDest.getGIDFieldOffsets(elementBlock,dstField);
359 
360  // field should be the same size
361  TEUCHOS_ASSERT(srcOffsets.size()==dstOffsets.size());
362  for(std::size_t i=0;i<srcOffsets.size();i++) {
363  offsetMap[srcOffsets[i]] = dstOffsets[i];
364 
365  // provides a size for allocating an array below: we will be able
366  // to index into dstFromSrcMap_ in a simple way
367  maxDest = dstOffsets[i]>maxDest ? dstOffsets[i] : maxDest;
368  }
369  }
370 
371  // Build map
372  TEUCHOS_ASSERT(maxDest>0);
373  dstFromSrcMap_ = std::vector<int>(maxDest+1,-1);
374  for(std::map<int,int>::const_iterator itr=offsetMap.begin();
375  itr!=offsetMap.end();++itr) {
376  dstFromSrcMap_[itr->second] = itr->first;
377  }
378 }
379 
380 // **********************************************************************
381 
382 #endif
383 
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
virtual Teuchos::RCP< Teuchos::Comm< int > > getComm() const =0
std::vector< PHX::MDField< ScalarT > > outFields_
std::vector< PHX::MDField< const ScalarT > > inFields_
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
Reorders the ad values of a specified field to match a different unique global indexer.
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.
virtual const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const =0
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)