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_UniqueGlobalIndexer.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 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>
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 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>
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 UniqueGlobalIndexerBase & indexerSrc,
131  const UniqueGlobalIndexerBase & 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 UniqueGlobalIndexerBase & indexerSrc,
165  const UniqueGlobalIndexerBase & 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  }
180 
181  // build a int-int map that associates fields
182  for(std::size_t i=0;i<inDOFs.size();i++) {
183  int srcFieldNum = indexerSrc.getFieldNum(inDOFs[i]);
184  int dstFieldNum = indexerDest.getFieldNum(outDOFs[i]);
185  TEUCHOS_ASSERT(srcFieldNum>=0);
186  TEUCHOS_ASSERT(dstFieldNum>=0);
187 
188  fieldNumberMaps[srcFieldNum] = dstFieldNum;
189  }
190 
191  buildSrcToDestMap(elementBlock,
192  fieldNumberMaps,
193  indexerSrc,
194  indexerDest);
195 
196  this->setName("Reorder AD Values");
197 }
198 
199 // **********************************************************************
200 template<typename TRAITS>
201 void panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
202 evaluateFields(typename TRAITS::EvalData /* workset */)
203 {
204  // for AD data do a reordering
205 
206  // 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!");
207 
208  for(std::size_t fieldIndex = 0; fieldIndex < inFields_.size(); ++fieldIndex) {
209 
210  const PHX::MDField<const ScalarT>& inField = inFields_[fieldIndex];
211  const PHX::MDField<ScalarT>& outField = outFields_[fieldIndex];
212 
213  if(inField.size()>0) {
214 
215  switch (inFields_[fieldIndex].rank()) {
216  case (1):
217  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i) {
218  outField(i).val() = inField(i).val();
219  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
220  outField(i).fastAccessDx(dx) = inField(i).fastAccessDx(dstFromSrcMap_[dx]);
221  }
222  break;
223  case (2):
224  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
225  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j) {
226  outField(i,j).val() = inField(i,j).val();
227  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
228  outField(i,j).fastAccessDx(dx) = inField(i,j).fastAccessDx(dstFromSrcMap_[dx]);
229  }
230  break;
231  case (3):
232  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
233  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
234  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k) {
235  outField(i,j,k).val() = inField(i,j,k).val();
236  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
237  outField(i,j,k).fastAccessDx(dx) = inField(i,j,k).fastAccessDx(dstFromSrcMap_[dx]);
238  }
239  break;
240  case (4):
241  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
242  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
243  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
244  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l) {
245  outField(i,j,k,l).val() = inField(i,j,k,l).val();
246  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
247  outField(i,j,k,l).fastAccessDx(dx) = inField(i,j,k,l).fastAccessDx(dstFromSrcMap_[dx]);
248  }
249  break;
250  case (5):
251  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
252  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
253  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
254  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l)
255  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.extent(4); ++m) {
256  outField(i,j,k,l,m).val() = inField(i,j,k,l,m).val();
257  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
258  outField(i,j,k,l,m).fastAccessDx(dx) = inField(i,j,k,l,m).fastAccessDx(dstFromSrcMap_[dx]);
259  }
260  break;
261  case (6):
262  for (typename PHX::MDField<ScalarT>::size_type i = 0; i < inField.extent(0); ++i)
263  for (typename PHX::MDField<ScalarT>::size_type j = 0; j < inField.extent(1); ++j)
264  for (typename PHX::MDField<ScalarT>::size_type k = 0; k < inField.extent(2); ++k)
265  for (typename PHX::MDField<ScalarT>::size_type l = 0; l < inField.extent(3); ++l)
266  for (typename PHX::MDField<ScalarT>::size_type m = 0; m < inField.extent(4); ++m)
267  for (typename PHX::MDField<ScalarT>::size_type n = 0; n < inField.extent(5); ++n) {
268  outField(i,j,k,l,m,n).val() = inField(i,j,k,l,m,n).val();
269  for (typename PHX::MDField<ScalarT>::size_type dx = 0; dx < Teuchos::as<typename PHX::MDField<ScalarT>::size_type>(dstFromSrcMap_.size()); ++dx)
270  outField(i,j,k,l,m,n).fastAccessDx(dx) = inField(i,j,k,l,m,n).fastAccessDx(dstFromSrcMap_[dx]);
271  }
272  break;
273  }
274 
275  }
276 
277  }
278 
279 //Irina TOFIX
280 /*
281  for(std::size_t i = 0; i < inFields_.size(); ++i) {
282 
283  for(typename PHX::MDField<ScalarT>::size_type j = 0; j < inFields_[i].size(); ++j) {
284  // allocated scalar fields
285  outFields_[i][j] = ScalarT(dstFromSrcMap_.size(), inFields_[i][j].val());
286 
287  ScalarT & outField = outFields_[i][j];
288  const ScalarT & inField = inFields_[i][j];
289 
290  // the jacobian must be initialized, otherwise its just a value copy
291  if(inField.size()>0) {
292  // loop over the sensitivity indices: all DOFs on a cell
293  outField.resize(dstFromSrcMap_.size());
294 
295  // copy jacobian entries correctly reordered
296  for(std::size_t k=0;k<dstFromSrcMap_.size();k++)
297  outField.fastAccessDx(k) = inField.fastAccessDx(dstFromSrcMap_[k]);
298  }
299 
300  outField.val() = inField.val();
301  }
302  }
303 */
304 }
305 
306 // **********************************************************************
307 template<typename TRAITS>
308 void panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
309 buildSrcToDestMap(const std::string & elementBlock,
310  const UniqueGlobalIndexerBase & indexerSrc,
311  const UniqueGlobalIndexerBase & indexerDest)
312 {
313  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
314  out.setOutputToRootOnly(0);
315 
316  TEUCHOS_ASSERT(indexerSrc.getComm()!=Teuchos::null);
317  TEUCHOS_ASSERT(indexerDest.getComm()!=Teuchos::null);
318 
319  const std::vector<int> & dstFieldsNum = indexerDest.getBlockFieldNumbers(elementBlock);
320 
321  // build a map between destination field numbers and source field numbers
322  std::map<int,int> fieldNumberMaps;
323  for(std::size_t i=0;i<dstFieldsNum.size();i++) {
324  std::string fieldName = indexerDest.getFieldString(dstFieldsNum[i]);
325 
326  int srcFieldNum = indexerSrc.getFieldNum(fieldName);
327  if(srcFieldNum>=0)
328  fieldNumberMaps[srcFieldNum] = dstFieldsNum[i];
329  else
330  out << "Warning: Reorder AD Values can't find field \"" << fieldName << "\"" << std::endl;
331  }
332 
333  buildSrcToDestMap(elementBlock,fieldNumberMaps,indexerSrc,indexerDest);
334 }
335 
336 // **********************************************************************
337 template<typename TRAITS>
338 void panzer::ReorderADValues_Evaluator<typename TRAITS::Jacobian, TRAITS>::
339 buildSrcToDestMap(const std::string & elementBlock,
340  const std::map<int,int> & fieldNumberMaps,
341  const UniqueGlobalIndexerBase & indexerSrc,
342  const UniqueGlobalIndexerBase & indexerDest)
343 {
344  int maxDest = -1;
345  std::map<int,int> offsetMap; // map from source to destination offsets
346  for(std::map<int,int>::const_iterator itr=fieldNumberMaps.begin();
347  itr!=fieldNumberMaps.end();++itr) {
348  int srcField = itr->first;
349  int dstField = itr->second;
350 
351  const std::vector<int> & srcOffsets = indexerSrc.getGIDFieldOffsets(elementBlock,srcField);
352  const std::vector<int> & dstOffsets = indexerDest.getGIDFieldOffsets(elementBlock,dstField);
353 
354  // field should be the same size
355  TEUCHOS_ASSERT(srcOffsets.size()==dstOffsets.size());
356  for(std::size_t i=0;i<srcOffsets.size();i++) {
357  offsetMap[srcOffsets[i]] = dstOffsets[i];
358 
359  // provides a size for allocating an array below: we will be able
360  // to index into dstFromSrcMap_ in a simple way
361  maxDest = dstOffsets[i]>maxDest ? dstOffsets[i] : maxDest;
362  }
363  }
364 
365  // Build map
366  TEUCHOS_ASSERT(maxDest>0);
367  dstFromSrcMap_ = std::vector<int>(maxDest+1,-1);
368  for(std::map<int,int>::const_iterator itr=offsetMap.begin();
369  itr!=offsetMap.end();++itr) {
370  dstFromSrcMap_[itr->second] = itr->first;
371  }
372 }
373 
374 // **********************************************************************
375 
376 #endif
377 
panzer::UniqueGlobalIndexerBase::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::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 UniqueGlobalIndexerBase &indexerSrc, const UniqueGlobalIndexerBase &indexerDest)
Definition: Panzer_ReorderADValues_Evaluator_impl.hpp:57
panzer::UniqueGlobalIndexerBase::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
Teuchos::basic_FancyOStream::setOutputToRootOnly
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
panzer::UniqueGlobalIndexerBase::getBlockFieldNumbers
virtual const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const =0
Panzer_UniqueGlobalIndexer.hpp
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::UniqueGlobalIndexerBase::getFieldString
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
panzer::UniqueGlobalIndexerBase::getFieldNum
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
Teuchos_RCP.hpp
Generated by   doxygen 1.8.5