Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessCombineAvg.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessCombineAvg.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Combines two fld files containing average fields.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <string>
38 using namespace std;
39 
40 #include "ProcessCombineAvg.h"
41 
44 
45 namespace Nektar
46 {
47 namespace FieldUtils
48 {
49 
50 ModuleKey ProcessCombineAvg::className =
52  ModuleKey(eProcessModule, "combineAvg"),
53  ProcessCombineAvg::create,
54  "combine two fields containing averages (and possibly Reynolds "
55  "stresses). Must specify fromfld.");
56 
57 ProcessCombineAvg::ProcessCombineAvg(FieldSharedPtr f) : ProcessModule(f)
58 {
59  m_config["fromfld"] =
60  ConfigOption(false, "NotSet", "Fld file form which to add field");
61 
62  ASSERTL0(m_config["fromfld"].as<string>().compare("NotSet") != 0,
63  "Need to specify fromfld=file.fld ");
64 }
65 
67 {
68 }
69 
70 void ProcessCombineAvg::Process(po::variables_map &vm)
71 {
72  if (m_f->m_verbose)
73  {
74  if (m_f->m_comm->TreatAsRankZero())
75  {
76  cout << "ProcessCombineAvg: Combining new fld into input avg fld..."
77  << endl;
78  }
79  }
80 
81  ASSERTL0(m_f->m_exp.size() != 0, "No input expansion defined");
82 
83  int nfields = m_f->m_fielddef[0]->m_fields.size();
84  int nq = m_f->m_exp[0]->GetTotPoints();
85  int expdim = m_f->m_graph->GetMeshDimension();
86  int spacedim = expdim;
87  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
88  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
89  {
90  spacedim += m_f->m_fielddef[0]->m_numHomogeneousDir;
91  }
92 
93  // Allocate storage for new field and correction (for Reynolds stress)
94  Array<OneD, Array<OneD, NekDouble> > fromPhys(nfields);
95  Array<OneD, Array<OneD, NekDouble> > correction(nfields);
96  for (int j = 0; j < nfields; ++j)
97  {
98  fromPhys[j] = Array<OneD, NekDouble>(nq, 0.0);
99  correction[j] = Array<OneD, NekDouble>(nq, 0.0);
100  }
101 
102  string fromfld = m_config["fromfld"].as<string>();
103  FieldSharedPtr fromField = boost::shared_ptr<Field>(new Field());
104  LibUtilities::FieldMetaDataMap fromFieldMetaDataMap;
105 
106  // Set up ElementGIDs in case of parallel processing
107  Array<OneD, int> ElementGIDs(m_f->m_exp[0]->GetExpSize());
108  for (int i = 0; i < m_f->m_exp[0]->GetExpSize(); ++i)
109  {
110  ElementGIDs[i] = m_f->m_exp[0]->GetExp(i)->GetGeom()->GetGlobalID();
111  }
112  // Import fromfld file
113  m_f->FieldIOForFile(fromfld)->Import(
114  fromfld, fromField->m_fielddef, fromField->m_data, fromFieldMetaDataMap,
115  ElementGIDs);
116  ASSERTL0(fromField->m_fielddef[0]->m_fields.size() == nfields,
117  "Mismatch in number of fields");
118  // Extract data to fromPhys
119  for (int j = 0; j < nfields; ++j)
120  {
121  ASSERTL0(fromField->m_fielddef[0]->m_fields[j] ==
122  m_f->m_fielddef[0]->m_fields[j],
123  "Field names do not match.");
124 
125  // load new field (overwrite m_f->m_exp coeffs for now)
126  for (int i = 0; i < fromField->m_data.size(); ++i)
127  {
128  m_f->m_exp[j]->ExtractDataToCoeffs(
129  fromField->m_fielddef[i], fromField->m_data[i],
130  fromField->m_fielddef[i]->m_fields[j],
131  m_f->m_exp[j]->UpdateCoeffs());
132  }
133  m_f->m_exp[j]->BwdTrans(m_f->m_exp[j]->GetCoeffs(), fromPhys[j]);
134  }
135 
136  // Load number of samples in each file
137  ASSERTL0(m_f->m_fieldMetaDataMap.count("NumberOfFieldDumps") != 0,
138  "Missing NumberOfFieldDumps metadata.");
139  ASSERTL0(fromFieldMetaDataMap.count("NumberOfFieldDumps") != 0,
140  "Missing NumberOfFieldDumps metadata.");
141  string s_num;
142  s_num = m_f->m_fieldMetaDataMap["NumberOfFieldDumps"];
143  int na = atoi(s_num.c_str());
144  s_num = fromFieldMetaDataMap["NumberOfFieldDumps"];
145  int nb = atoi(s_num.c_str());
146 
147  // Look for Reynolds stresses
148  int stress = -1;
149  for (int j = 0; j < nfields; ++j)
150  {
151  if (m_f->m_fielddef[0]->m_fields[j] == "uu")
152  {
153  stress = j;
154  break;
155  }
156  }
157 
158  // Calculate correction for Reynolds stresses
159  if (stress != -1)
160  {
161  Array<OneD, NekDouble> tmp(nq, 0.0);
162  int n = stress;
163  // Follow same numbering as FilterReynoldsStresses
164  for (int i = 0; i < spacedim; ++i)
165  {
166  for (int j = i; j < spacedim; ++j, ++n)
167  {
168  // correction is zero for averages and
169  // = (\bar{x_a}-\bar{x_b})*(\bar{y_a}-\bar{y_b})*na*nb/N
170  // for Reynolds stresses
171  NekDouble fac = ((NekDouble)(na * nb)) / ((NekDouble)(na + nb));
172  Vmath::Vsub(nq, m_f->m_exp[i]->GetPhys(), 1, fromPhys[i], 1,
173  correction[n], 1);
174  Vmath::Vsub(nq, m_f->m_exp[j]->GetPhys(), 1, fromPhys[j], 1,
175  tmp, 1);
176  Vmath::Vmul(nq, correction[n], 1, tmp, 1, correction[n], 1);
177  Vmath::Smul(nq, fac, correction[n], 1, correction[n], 1);
178  }
179  }
180  }
181  // Combine fields
182  for (int j = 0; j < nfields; ++j)
183  {
184  // The new value is: (x_a*na + x_b*nb + correction)/N
185  Vmath::Smul(nq, 1.0 * na, m_f->m_exp[j]->GetPhys(), 1,
186  m_f->m_exp[j]->UpdatePhys(), 1);
187  Vmath::Svtvp(nq, 1.0 * nb, fromPhys[j], 1, m_f->m_exp[j]->GetPhys(), 1,
188  m_f->m_exp[j]->UpdatePhys(), 1);
189  Vmath::Vadd(nq, m_f->m_exp[j]->GetPhys(), 1, correction[j], 1,
190  m_f->m_exp[j]->UpdatePhys(), 1);
191  Vmath::Smul(nq, 1.0 / (na + nb), m_f->m_exp[j]->GetPhys(), 1,
192  m_f->m_exp[j]->UpdatePhys(), 1);
193 
194  m_f->m_exp[j]->FwdTrans_IterPerExp(m_f->m_exp[j]->GetPhys(),
195  m_f->m_exp[j]->UpdateCoeffs());
196  }
197 
198  // Update metadata
199  m_f->m_fieldMetaDataMap["NumberOfFieldDumps"] =
200  boost::lexical_cast<std::string>(na + nb);
201  NekDouble t0 = -1;
202  NekDouble finTime = -1;
203  if (m_f->m_fieldMetaDataMap.count("InitialTime"))
204  {
205  string s_t = m_f->m_fieldMetaDataMap["InitialTime"];
206  NekDouble t = atof(s_t.c_str());
207 
208  t0 = t;
209  }
210  if (fromFieldMetaDataMap.count("InitialTime"))
211  {
212  string s_t = fromFieldMetaDataMap["InitialTime"];
213  NekDouble t = atof(s_t.c_str());
214 
215  if (t0 == -1)
216  {
217  t0 = t;
218  }
219  else
220  {
221  t0 = std::min(t0, t);
222  }
223  }
224  if (m_f->m_fieldMetaDataMap.count("FinalTime"))
225  {
226  string s_t = m_f->m_fieldMetaDataMap["FinalTime"];
227  NekDouble t = atof(s_t.c_str());
228 
229  finTime = std::max(t0, t);
230  }
231  if (fromFieldMetaDataMap.count("FinalTime"))
232  {
233  string s_t = fromFieldMetaDataMap["FinalTime"];
234  NekDouble t = atof(s_t.c_str());
235 
236  finTime = std::max(t0, t);
237  }
238  if (t0 != -1)
239  {
240  m_f->m_fieldMetaDataMap["InitialTime"] =
241  boost::lexical_cast<std::string>(t0);
242  }
243  if (finTime != -1)
244  {
245  m_f->m_fieldMetaDataMap["FinalTime"] =
246  boost::lexical_cast<std::string>(finTime);
247  }
248 
249  // Update field def and data
250  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
251  m_f->m_exp[0]->GetFieldDefinitions();
252  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
253  for (int i = 0; i < nfields; ++i)
254  {
255  for (int j = 0; j < FieldDef.size(); ++j)
256  {
257  FieldDef[j]->m_fields.push_back(m_f->m_fielddef[0]->m_fields[i]);
258  m_f->m_exp[i]->AppendFieldData(FieldDef[j], FieldData[j]);
259  }
260  }
261 
262  m_f->m_fielddef = FieldDef;
263  m_f->m_data = FieldData;
264 }
265 }
266 }
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Represents a command-line configuration option.
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
STL namespace.
pair< ModuleType, string > ModuleKey
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:54
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:740
double NekDouble
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:343
virtual void Process(po::variables_map &vm)
Write mesh to output file.
Abstract base class for processing modules.
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215