Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessC0Projection.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessC0Projection.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: Computes C0 projection.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <string>
38 using namespace std;
39 
40 #include "ProcessC0Projection.h"
41 
44 
45 namespace Nektar
46 {
47 namespace FieldUtils
48 {
49 
50 ModuleKey ProcessC0Projection::className =
52  ModuleKey(eProcessModule, "C0Projection"),
53  ProcessC0Projection::create,
54  "Computes C0 projection.");
55 
56 ProcessC0Projection::ProcessC0Projection(FieldSharedPtr f) : ProcessModule(f)
57 {
58  m_config["fields"] = ConfigOption(false, "All", "Start field to project");
59  m_config["localtoglobalmap"] = ConfigOption(
60  true, "0", "Just perform a local to global mapping and back");
61  m_config["usexmlbcs"] = ConfigOption(
62  true, "0", "Use boundary conditions given in xml file. Requires all "
63  "projected fields to be defined in xml file");
64  m_config["helmsmoothing"] = ConfigOption(
65  false, "Not Set", "Use a Helmholtz smoother to remove high frequency "
66  "components above specified L");
67 
68  f->m_declareExpansionAsContField = true;
69 }
70 
72 {
73 }
74 
75 void ProcessC0Projection::Process(po::variables_map &vm)
76 {
77  if (m_f->m_verbose)
78  {
79  if (m_f->m_comm->TreatAsRankZero())
80  {
81  cout << "ProcessC0Projection: Projecting field into C0 space..."
82  << endl;
83  }
84  }
85 
86  // ensure not using diagonal preconditioner since tends not to converge fo
87  // mass matrix
88  if (m_f->m_graph->GetMeshDimension() == 3)
89  {
90  if (boost::iequals(m_f->m_session->GetSolverInfo("GLOBALSYSSOLN"),
91  "IterativeStaticCond"))
92  {
93  if (boost::iequals(m_f->m_session->GetSolverInfo("PRECONDITIONER"),
94  "Diagonal"))
95  {
96  m_f->m_session->SetSolverInfo("PRECONDITIONER",
97  "LowEnergyBlock");
98  }
99  if (boost::iequals(m_f->m_session->GetSolverInfo("PRECONDITIONER"),
100  "FullLinearSpaceWithDiagonal"))
101  {
102  m_f->m_session->SetSolverInfo(
103  "PRECONDITIONER", "FullLinearSpaceWithLowEnergyBlock");
104  }
105 
106  if (m_f->m_verbose)
107  {
108  if (m_f->m_comm->GetRank() == 0)
109  {
110  cout << "Resetting diagonal precondition to low energy "
111  "block "
112  << endl;
113  }
114  }
115  }
116  }
117  bool JustPerformLocToGloMap = m_config["localtoglobalmap"].as<bool>();
118  bool HelmSmoother =
119  (boost::iequals(m_config["helmsmoothing"].as<string>(), "Not Set"))
120  ? false
121  : true;
122  int nfields = m_f->m_exp.size();
123  Array<OneD, MultiRegions::ExpListSharedPtr> C0ProjectExp(nfields);
124  if (m_config["usexmlbcs"].as<bool>())
125  {
126  for (int i = 0; i < nfields; ++i)
127  {
128  C0ProjectExp[i] = m_f->m_exp[i];
129  }
130  }
131  else
132  {
133  // generate a C0 expansion field with no boundary conditions.
134  bool savedef = m_f->m_declareExpansionAsContField;
135  bool savedef2 = m_f->m_requireBoundaryExpansion;
136  m_f->m_declareExpansionAsContField = true;
137  m_f->m_requireBoundaryExpansion = false;
138  C0ProjectExp[0] = m_f->AppendExpList(
139  m_f->m_fielddef[0]->m_numHomogeneousDir, "DefaultVar", true);
140  m_f->m_declareExpansionAsContField = savedef;
141  m_f->m_requireBoundaryExpansion = savedef2;
142  for (int i = 1; i < nfields; ++i)
143  {
144  C0ProjectExp[i] = C0ProjectExp[0];
145  }
146  }
147 
148  string fields = m_config["fields"].as<string>();
149  vector<unsigned int> processFields;
150 
151  if (fields.compare("All") == 0)
152  {
153  for (int i = 0; i < nfields; ++i)
154  {
155  processFields.push_back(i);
156  }
157  }
158  else
159  {
160  ASSERTL0(
161  ParseUtils::GenerateOrderedVector(fields.c_str(), processFields),
162  "Failed to interpret field string in C0Projection");
163  }
164 
165  for (int i = 0; i < processFields.size(); ++i)
166  {
167  ASSERTL0(processFields[i] < nfields,
168  "Attempt to process field that is larger than then number of "
169  "fields available");
170 
171  if (m_f->m_verbose)
172  {
173  if (m_f->m_comm->GetRank() == 0)
174  {
175  cout << "\t Processing field: " << processFields[i] << endl;
176  }
177  }
178 
179  if (JustPerformLocToGloMap)
180  {
181  int ncoeffs = m_f->m_exp[0]->GetNcoeffs();
182  Vmath::Vcopy(ncoeffs, m_f->m_exp[processFields[i]]->GetCoeffs(), 1,
183  C0ProjectExp[processFields[i]]->UpdateCoeffs(), 1);
184  C0ProjectExp[processFields[i]]->LocalToGlobal();
185  C0ProjectExp[processFields[i]]->GlobalToLocal();
186  Vmath::Vcopy(ncoeffs, C0ProjectExp[processFields[i]]->GetCoeffs(),
187  1, m_f->m_exp[processFields[i]]->UpdateCoeffs(), 1);
188  }
189  else
190  {
191  if (HelmSmoother)
192  {
193  int dim = m_f->m_graph->GetSpaceDimension();
194  int npoints = m_f->m_exp[0]->GetNpoints();
195  NekDouble lambda = m_config["helmsmoothing"].as<NekDouble>();
196  lambda = 2 * M_PI / lambda;
197  lambda = lambda * lambda;
198 
199  if (m_f->m_verbose)
200  {
201  cout << "Setting up Helmholtz smoother with lambda = "
202  << lambda << endl;
203  }
204 
206  Array<OneD, NekDouble> forcing(npoints);
207  factors[StdRegions::eFactorLambda] = -lambda;
208 
209  Array<OneD, Array<OneD, NekDouble> > Velocity(dim);
210  for (int j = 0; j < dim; ++j)
211  {
212  Velocity[j] = Array<OneD, NekDouble>(npoints, 0.0);
213  }
214 
215  C0ProjectExp[processFields[i]]->BwdTrans(
216  m_f->m_exp[processFields[i]]->GetCoeffs(),
217  m_f->m_exp[processFields[i]]->UpdatePhys());
218 
219  Vmath::Smul(npoints, -lambda,
220  m_f->m_exp[processFields[i]]->GetPhys(), 1, forcing,
221  1);
222 
223  // Note we are using the
224  // LinearAdvectionDiffusionReaction solver here
225  // instead of HelmSolve since lambda is negative and
226  // so matrices are not positive definite. Ideally
227  // should allow for negative lambda coefficient in
228  // HelmSolve
229  C0ProjectExp[processFields[i]]
230  ->LinearAdvectionDiffusionReactionSolve(
231  Velocity, forcing,
232  m_f->m_exp[processFields[i]]->UpdateCoeffs(), -lambda);
233  }
234  else
235  {
236  C0ProjectExp[processFields[i]]->BwdTrans(
237  m_f->m_exp[processFields[i]]->GetCoeffs(),
238  m_f->m_exp[processFields[i]]->UpdatePhys());
239  C0ProjectExp[processFields[i]]->FwdTrans(
240  m_f->m_exp[processFields[i]]->GetPhys(),
241  m_f->m_exp[processFields[i]]->UpdateCoeffs());
242  }
243  }
244  }
245 
246  // reset FieldDef in case of serial input and parallel output
247  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
248  m_f->m_exp[0]->GetFieldDefinitions();
249  // reset up FieldData with new values before projecting
250  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
251 
252  for (int i = 0; i < nfields; ++i)
253  {
254  for (int j = 0; j < FieldDef.size(); ++j)
255  {
256  FieldDef[j]->m_fields.push_back(m_f->m_fielddef[0]->m_fields[i]);
257  m_f->m_exp[i]->AppendFieldData(FieldDef[j], FieldData[j]);
258  }
259  }
260 
261  m_f->m_fielddef = FieldDef;
262  m_f->m_data = FieldData;
263 }
264 }
265 }
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:97
Represents a command-line configuration option.
STL namespace.
pair< ModuleType, string > ModuleKey
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
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
Abstract base class for processing modules.
virtual void Process(po::variables_map &vm)
Write mesh to output file.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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