Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_GlobalStatistics_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_GLOBAL_STATISTICS_IMPL_HPP
44#define PANZER_GLOBAL_STATISTICS_IMPL_HPP
45
46#include "Intrepid2_FunctionSpaceTools.hpp"
50#include "Panzer_GlobalData.hpp"
52#include "Phalanx_DataLayout_MDALayout.hpp"
53#include "Teuchos_ScalarTraits.hpp"
54#include "Teuchos_CommHelpers.hpp"
55#include <iomanip>
56
57namespace panzer {
58
59//**********************************************************************
60template<typename EvalT, typename Traits>
63 const Teuchos::ParameterList& p)
64{
65 comm = p.get< Teuchos::RCP<const Teuchos::Comm<int> > >("Comm");
66
67 global_data = p.get<Teuchos::RCP<panzer::GlobalData> >("Global Data");
68
69 // Expects a string that is a Colon separated list of field names to compute statistics on.
70 // for example the string "UX:UY:UZ:PRESSURE" would be separated into a vector with
71 // four fields, "UX", "UY", "UZ", and "PRESSURE".
72 std::string names_string = p.get<std::string>("Names");
73 std::vector<std::string> names;
74 panzer::StringTokenizer(names, names_string);
75
76 Teuchos::RCP<panzer::IntegrationRule> ir = p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR");
77
78 field_values.clear();
79 for (typename std::vector<std::string>::const_iterator name = names.begin(); name != names.end(); ++name)
80 field_values.push_back(PHX::MDField<const ScalarT,Cell,IP>(*name, ir->dl_scalar));
81
82 Teuchos::RCP<PHX::MDALayout<Cell> > cell_dl = Teuchos::rcp(new PHX::MDALayout<Cell>(ir->dl_scalar->extent(0)));
83 volumes = PHX::MDField<ScalarT,Cell>("Cell Volumes",cell_dl);
84
85 tmp = PHX::MDField<ScalarT,Cell>("GlobalStatistics:tmp:"+names_string,cell_dl);
86 ones = PHX::MDField<ScalarT,Cell,IP>("GlobalStatistics:ones:"+names_string,ir->dl_scalar);
87
88 this->addEvaluatedField(volumes);
89 this->addEvaluatedField(tmp);
90 this->addEvaluatedField(ones);
91 for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::const_iterator field = field_values.begin();
92 field != field_values.end(); ++field) {
93 this->addDependentField(*field);
94 }
95
96 averages.resize(field_values.size());
97 maxs.resize(field_values.size());
98 mins.resize(field_values.size());
99 global_maxs.resize(field_values.size());
100 global_mins.resize(field_values.size());
101 global_averages.resize(field_values.size());
102
103 ir_order = ir->cubature_degree;
104
105 std::string n = "GlobalStatistics: " + names_string;
106 this->setName(n);
107}
108
109//**********************************************************************
110template<typename EvalT, typename Traits>
111void
114 typename Traits::SetupData sd,
116{
117 ir_index = panzer::getIntegrationRuleIndex(ir_order,(*sd.worksets_)[0], this->wda);
118 auto l_ones = ones.get_static_view();
119 Kokkos::parallel_for("GlobalStatistics", l_ones.extent(0), KOKKOS_LAMBDA(int cell) {
120 for (std::size_t ip = 0; ip < l_ones.extent(1); ++ip)
121 l_ones(cell,ip) = 1.0;
122 });
123}
124
125//**********************************************************************
126template<typename EvalT, typename Traits>
127void
130 typename Traits::EvalData workset)
131{
132 if (workset.num_cells == 0)
133 return;
134
135 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::integrate(volumes.get_view(),
136 ones.get_view(),
137 (this->wda(workset).int_rules[ir_index])->weighted_measure.get_view());
138 auto volumes_h = Kokkos::create_mirror_view(as_view(volumes));
139 Kokkos::deep_copy(volumes_h, as_view(volumes));
140
141 for (index_t cell = 0; cell < workset.num_cells; ++cell)
142 total_volume += volumes_h(cell);
143
144 typename std::vector<PHX::MDField<ScalarT,Cell,IP> >::size_type field_index = 0;
145 for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::iterator field = field_values.begin();
146 field != field_values.end(); ++field,++field_index) {
147
148 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::integrate(tmp.get_view(),
149 field->get_view(),
150 (this->wda(workset).int_rules[ir_index])->weighted_measure.get_view());
151 auto tmp_h = Kokkos::create_mirror_view(tmp.get_static_view());
152 auto field_h = Kokkos::create_mirror_view( field->get_static_view());
153 Kokkos::deep_copy(tmp_h, tmp.get_static_view());
154 Kokkos::deep_copy(field_h, field->get_static_view());
155
156
157 for (index_t cell = 0; cell < workset.num_cells; ++cell) {
158 averages[field_index] += tmp_h(cell);
159
160 for (typename PHX::MDField<ScalarT,Cell,IP>::size_type ip = 0; ip < (field->extent(1)); ++ip) {
161 maxs[field_index] = std::max( field_h(cell,ip), maxs[field_index]);
162 mins[field_index] = std::min( field_h(cell,ip), mins[field_index]);
163 }
164 }
165
166 }
167}
168
169//**********************************************************************
170template<typename EvalT, typename Traits>
171void
174 typename Traits::PreEvalData /* data */)
175{
176 total_volume = Teuchos::ScalarTraits<ScalarT>::zero();
177
178 for (typename std::vector<ScalarT>::iterator field = averages.begin(); field != averages.end(); ++field)
179 *field = Teuchos::ScalarTraits<ScalarT>::zero();
180
181 for (typename std::vector<ScalarT>::iterator field = maxs.begin(); field != maxs.end(); ++field)
182 *field = Teuchos::ScalarTraits<ScalarT>::rmin();
183
184 for (typename std::vector<ScalarT>::iterator field = mins.begin(); field != mins.end(); ++field)
185 *field = Teuchos::ScalarTraits<ScalarT>::rmax();
186}
187
188//**********************************************************************
189template<typename EvalT, typename Traits>
190void
193 typename Traits::PostEvalData /* data */)
194{
195 this->postprocess(*(global_data->os));
196}
197
198//**********************************************************************
199template<typename EvalT, typename TRAITS>
201{
202 // throw unless specialized for residual evaluations
203 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"SHOULD NEVER BE CALLED!");
204}
205
206//**********************************************************************
207template<>
209{
210 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, static_cast<int>(1), &total_volume, &global_total_volume);
211 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, static_cast<int>(averages.size()), &averages[0], &global_averages[0]);
212 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, static_cast<int>(maxs.size()), &maxs[0], &global_maxs[0]);
213 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MIN, static_cast<int>(mins.size()), &mins[0], &global_mins[0]);
214
215 for (std::vector<ScalarT>::size_type i = 0; i < field_values.size(); ++i)
216 global_averages[i] /= global_total_volume;
217
218 if (comm->getRank() == 0) {
219
220 panzer::ios_all_saver saver(os);
221
222 std::size_t precision = 8;
223 os << std::scientific << std::showpoint << std::setprecision(precision) << std::left;
224
225 std::size_t name_width = 0;
226 for (std::vector<ScalarT>::size_type i = 0; i < field_values.size(); ++i)
227 name_width = std::max(name_width,field_values[i].fieldTag().name().size());
228
229 std::size_t value_width = precision + 7;
230
231 os << std::setw(name_width) << "Field"
232 << " " << std::setw(value_width) << "Average"
233 << " " << std::setw(value_width) << "Maximum (@IP)"
234 << " " << std::setw(value_width) << "Minimum (@IP)"
235 << std::endl;
236
237 for (std::vector<ScalarT>::size_type i = 0; i < field_values.size(); ++i) {
238 os << std::setw(name_width) << field_values[i].fieldTag().name()
239 << " " << std::setw(value_width) << global_averages[i]
240 << " " << std::setw(value_width) << global_maxs[i]
241 << " " << std::setw(value_width) << global_mins[i] << std::endl;
242 }
243
244 }
245
246}
247
248//**********************************************************************
249template<typename EvalT, typename TRAITS>
251{
252 return tmp.fieldTag();
253}
254
255
256} // namespace panzer
257
258
259
260#endif
261
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
GlobalStatistics(const Teuchos::ParameterList &p)
void postEvaluate(typename Traits::PostEvalData d)
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
const PHX::FieldTag & getRequiredFieldTag()
void preEvaluate(typename Traits::PreEvalData d)
void evaluateFields(typename Traits::EvalData d)
int num_cells
DEPRECATED - use: numCells()
void StringTokenizer(std::vector< std::string > &tokens, const std::string &str, const std::string delimiters, bool trim)
Tokenize a string, put tokens in a vector.
std::vector< int >::size_type getIntegrationRuleIndex(int ir_degree, const panzer::Workset &workset, WorksetDetailsAccessor &wda)
Teuchos::RCP< const std::vector< panzer::Workset > > worksets_