Nektar++
ForcingWavyness.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ForcingWavyness.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: Implementation of a forcing term to simulate a wavy body
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <MultiRegions/ExpList.h>
38 
39 namespace Nektar
40 {
41 namespace SolverUtils
42 {
43 
44 std::string ForcingWavyness::className =
46  ForcingWavyness::create, "Wavyness Forcing");
47 
48 /**
49  *
50  */
53  : Forcing(pSession)
54 {
55 }
56 
57 
58 /**
59  *
60  */
63  const unsigned int &pNumForcingFields,
64  const TiXmlElement *pForce)
65 {
66  // Just 3D homogenous 1D problems can use this techinque
67  ASSERTL0(pFields[0]->GetExpType() == MultiRegions::e3DH1D,
68  "Wavyness implemented just for 3D Homogenous 1D expansions.");
69 
70  m_session->LoadParameter("Kinvis", m_kinvis);
71  // forcing size (it must be 3)
72  m_NumVariable = pNumForcingFields;
75  int phystot = pFields[0]->GetTotPoints();
76 
77  // Allocation of forcing and geometry memory
78  for (int i = 0; i < m_NumVariable; ++i)
79  {
80  m_Forcing[i] = Array<OneD, NekDouble>(phystot, 0.0);
81  }
82  for (int i = 0; i < m_wavyGeom.num_elements(); i++)
83  {
84  m_wavyGeom[i] = Array<OneD, NekDouble>(phystot, 0.0);
85  }
86 
87  const TiXmlElement* funcNameElmt;
88  funcNameElmt = pForce->FirstChildElement("WAVYNESS");
89  ASSERTL0(funcNameElmt, "Requires WAVYNESS tag, specifying function "
90  "name which prescribes wavy function.");
91 
92  string funcName = funcNameElmt->GetText();
93  ASSERTL0(m_session->DefinesFunction(funcName),
94  "Function '" + funcName + "' not defined.");
95 
96  std::string s_FieldStr = m_session->GetVariable(0);
97  ASSERTL0(m_session->DefinesFunction(funcName, s_FieldStr),
98  "Variable '" + s_FieldStr + "' not defined.");
99 
100  EvaluateFunction(pFields, m_session, s_FieldStr, m_wavyGeom[0],
101  funcName);
102 
103  // Calculate derivatives of transformation
104  for (int i = 1; i < 4; i++)
105  {
106  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
107  m_wavyGeom[i - 1], m_wavyGeom[i]);
108  }
109 
110  Vmath::Vmul(phystot, m_wavyGeom[1], 1, m_wavyGeom[1], 1,
111  m_wavyGeom[4], 1);
112  Vmath::Vmul(phystot, m_wavyGeom[1], 1, m_wavyGeom[2], 1,
113  m_wavyGeom[5], 1);
114 }
115 
116 
117 /**
118  *
119  */
122  const Array<OneD, Array<OneD, NekDouble> > &inarray,
123  Array<OneD, Array<OneD, NekDouble> > &outarray,
124  const NekDouble &time)
125 {
126  //calcualte the forcing components Ax,Ay,Az and put them in m_Forcing
127  CalculateForcing(fields);
128 
129  // Apply forcing terms
130  for (int i = 0; i < m_NumVariable; i++)
131  {
132  Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1, m_Forcing[i], 1,
133  outarray[i], 1);
134  }
135 }
136 
137 
138 /**
139  *
140  */
143 {
144  int np = fields[0]->GetNpoints();
145 
146  Array<OneD, NekDouble> U, V, W, P, tmp1, tmp2, tmp3, Wz, Wzz, Px;
147 
148  U = Array<OneD, NekDouble>(np);
149  V = Array<OneD, NekDouble>(np);
150  W = Array<OneD, NekDouble>(np);
151  P = Array<OneD, NekDouble>(np);
152 
153  tmp1 = Array<OneD, NekDouble>(np);
154  tmp2 = Array<OneD, NekDouble>(np);
155  tmp3 = Array<OneD, NekDouble>(np);
156 
157  Wz = Array<OneD, NekDouble>(np);
158  Wzz = Array<OneD, NekDouble>(np);
159  Px = Array<OneD, NekDouble>(np);
160 
161  fields[0]->HomogeneousBwdTrans(fields[0]->GetPhys(), U);
162  fields[0]->HomogeneousBwdTrans(fields[1]->GetPhys(), V);
163  fields[0]->HomogeneousBwdTrans(fields[2]->GetPhys(), W);
164  fields[0]->HomogeneousBwdTrans(fields[3]->GetPhys(), P);
165  //-------------------------------------------------------------------------
166  // Ax calculation
167  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], P, Px); // Px
168  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
169  fields[3]->GetPhys(), tmp2); // Pz
170  // Pz in physical space
171  fields[0]->HomogeneousBwdTrans(tmp2, tmp3);
172  // Pz * Xz
173  Vmath::Vmul(np, tmp3, 1, m_wavyGeom[1], 1, tmp3, 1);
174  // Px * Xz^2
175  Vmath::Vmul(np, Px, 1, m_wavyGeom[4], 1, tmp1, 1);
176  // W^2
177  Vmath::Vmul(np, W, 1, W, 1, tmp2, 1);
178  // Xzz * W^2
179  Vmath::Vmul(np, tmp2, 1, m_wavyGeom[2], 1, tmp2, 1);
180  // Pz * Xz - Px * Xz^2
181  Vmath::Vsub(np, tmp3, 1, tmp1, 1, m_Forcing[0], 1);
182  // A0 = Pz * Xz - Px * Xz^2 - Xzz * W^2
183  Vmath::Vsub(np, m_Forcing[0], 1, tmp2, 1, m_Forcing[0], 1);
184 
185  // here part to be multiplied by 1/Re we use P to store it, since we dont't
186  // need it anymore
187  // W * Xzzz
188  Vmath::Vmul(np, W, 1, m_wavyGeom[3], 1, P, 1);
189  // P = W * Xzzz + Xz * Xzz
190  Vmath::Vadd(np, P, 1, m_wavyGeom[5], 1, P, 1);
191  // Wz
192  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
193  fields[2]->GetPhys(), tmp1);
194  // Wzz
195  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2], tmp1, tmp2);
196  //Wz in physical space
197  fields[0]->HomogeneousBwdTrans(tmp1, Wz);
198  //Wzz in physical space
199  fields[0]->HomogeneousBwdTrans(tmp2, Wzz);
200  // Wzz * Xz
201  Vmath::Vmul(np, Wzz, 1, m_wavyGeom[1], 1, tmp1, 1);
202  // P = W * Xzzz + Xz * Xzz - Wzz * Xz
203  Vmath::Vsub(np, P, 1, tmp1, 1, P, 1);
204  // Ux
205  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], U, tmp1);
206  // Ux * Xzz
207  Vmath::Vmul(np, tmp1, 1, m_wavyGeom[2], 1, tmp2, 1);
208  // P = W * Xzzz + Xz * Xzz - Wzz * Xz - Ux * Xzz
209  Vmath::Vsub(np, P, 1, tmp2, 1, P, 1);
210  // Uxx
211  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp1, tmp2);
212  // Uxx * Xz^2
213  Vmath::Vmul(np, tmp2, 1, m_wavyGeom[4], 1, tmp2, 1);
214  // P = W * Xzzz + Xz * Xzz - Wzz * Xz - Ux * Xzz + Uxx * Xz^2
215  Vmath::Vadd(np, P, 1, tmp2, 1, P, 1);
216  // Uz
217  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
218  fields[0]->GetPhys(), tmp1);
219  // Uzx
220  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp1, tmp2);
221  // Uzx in physical space
222  fields[0]->HomogeneousBwdTrans(tmp2, tmp1);
223  // Uzx * Xz
224  Vmath::Vmul(np, tmp1, 1, m_wavyGeom[1], 1, tmp1, 1);
225  // 2 * Uzx * Xz
226  Vmath::Smul(np, 2.0, tmp1, 1, tmp2, 1);
227  // P = W * Xzzz + Xz * Xzz - Wzz * Xz - Ux * Xzz + Uxx * Xz^2 - 2 * Uzx * Xz
228  Vmath::Vsub(np, P, 1, tmp2, 1, P, 1);
229  // *1/Re
230  Vmath::Smul(np, m_kinvis, P, 1, tmp1, 1);
231  Vmath::Vadd(np, tmp1, 1, m_Forcing[0], 1, tmp2, 1);
232  // back to Fourier Space
233  fields[0]->HomogeneousFwdTrans(tmp2, m_Forcing[0]);
234 
235  //-------------------------------------------------------------------------
236  // Ay calucaltion
237 
238  // Vx
239  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], V, tmp1);
240  // Vxx
241  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp1, tmp2);
242  // Vx * Xzz
243  Vmath::Vmul(np, tmp1, 1, m_wavyGeom[2], 1, tmp1, 1);
244  // Vxx * Xz^2
245  Vmath::Vmul(np, tmp2, 1, m_wavyGeom[4], 1, tmp2, 1);
246  // Vxx * Xz^2 - Vx* Xzz
247  Vmath::Vsub(np, tmp2, 1, tmp1, 1, m_Forcing[1], 1);
248  //Vz
249  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
250  fields[1]->GetPhys(), tmp3);
251  //Vzx
252  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp3, tmp2);
253  // Vzx physical space
254  fields[0]->HomogeneousBwdTrans(tmp2, tmp1);
255  // Vzx * Xz
256  Vmath::Vmul(np, tmp1, 1, m_wavyGeom[1], 1, tmp2, 1);
257  // 2 * Vzx * Xz
258  Vmath::Smul(np, 2.0, tmp2, 1, tmp1, 1);
259  // Vxx * Xz^2 - Vx* Xzz - Vxz * Xz
260  Vmath::Vsub(np, m_Forcing[1], 1, tmp1, 1, tmp3, 1);
261  // * 1/Re
262  Vmath::Smul(np, m_kinvis, tmp3, 1, tmp1, 1);
263  // back to Fourier Space
264  fields[0]->HomogeneousFwdTrans(tmp1, m_Forcing[1]);
265 
266  //-------------------------------------------------------------------------
267  // Az calculation
268  // Px * Xz
269  Vmath::Vmul(np, Px, 1, m_wavyGeom[1], 1, P, 1);
270  // Wzz - Xzz
271  Vmath::Vsub(np, Wzz, 1, m_wavyGeom[2], 1, Wzz, 1);
272  // Wx
273  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], W, tmp1);
274  // Wx * Xzz
275  Vmath::Vmul(np, tmp1, 1, m_wavyGeom[2], 1, tmp2, 1);
276  // Wzz - Xzz - Wx * Xzz
277  Vmath::Vsub(np, Wzz, 1, tmp2, 1, Wzz, 1);
278  // Wxx
279  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp1, tmp2);
280  // Wxx * Xz^2
281  Vmath::Vmul(np, tmp2, 1, m_wavyGeom[4], 1, tmp2, 1);
282  // Wzz - Xzz - Wx * Xzz + Wxx * Xz^2
283  Vmath::Vadd(np, Wzz, 1, tmp2, 1, Wzz, 1);
284  // Wzx
285  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], Wz, tmp1);
286  // Wzx * Xz
287  Vmath::Vmul(np, tmp1, 1, m_wavyGeom[1], 1, tmp1, 1);
288  // 2 * Vzx * Xz
289  Vmath::Smul(np, 2.0, tmp1, 1, tmp1, 1);
290  // Wzz - Xzz - Wx * Xzz + Wxx * Xz^2 - 2 * Vzx * Xz
291  Vmath::Vsub(np, Wzz, 1, tmp1, 1, Wzz, 1);
292  // * 1/Re
293  Vmath::Smul(np, m_kinvis, Wzz, 1, tmp1, 1);
294  Vmath::Vadd(np, tmp1, 1, P, 1, tmp2, 1);
295  // back to Fourier Space
296  fields[0]->HomogeneousFwdTrans(tmp2, m_Forcing[2]);
297 }
298 
299 }
300 }
void CalculateForcing(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
static SOLVER_UTILS_EXPORT ForcingSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
Creates an instance of this class.
Array< OneD, Array< OneD, NekDouble > > m_Forcing
Evaluated forcing function.
Definition: Forcing.h:97
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, Array< OneD, NekDouble > > m_wavyGeom
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:42
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, MultiRegions::ExpListSharedPtr > pFields, LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
Definition: Forcing.cpp:140
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
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:199
double NekDouble
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
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:329
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
int m_NumVariable
Number of variables.
Definition: Forcing.h:99
ForcingWavyness(const LibUtilities::SessionReaderSharedPtr &pSession)
virtual SOLVER_UTILS_EXPORT void v_Apply(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:70
virtual SOLVER_UTILS_EXPORT void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
static std::string className
Name of the class.
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:285
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:169
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215