Panzer Version of the Day
Loading...
Searching...
No Matches
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
55template<typename EvalT,typename TRAITS>
57ReorderADValues_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// **********************************************************************
80template<typename EvalT,typename TRAITS>
82ReorderADValues_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// **********************************************************************
108template<typename EvalT,typename TRAITS>
110evaluateFields(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
124template<typename TRAITS>
126ReorderADValues_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
156template<typename TRAITS>
158ReorderADValues_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// **********************************************************************
202template<typename TRAITS>
204evaluateFields(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// **********************************************************************
313template<typename TRAITS>
315buildSrcToDestMap(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// **********************************************************************
343template<typename TRAITS>
345buildSrcToDestMap(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
Reorders the ad values of a specified field to match a different unique global indexer.
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)