Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NavierStokesAdvection.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NavierStokesAdvection.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: Evaluation of the Navier Stokes advective term
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <IncNavierStokesSolver/AdvectionTerms/NavierStokesAdvection.h>
37 
38 namespace Nektar
39 {
40  string NavierStokesAdvection::className = GetAdvectionTermFactory().RegisterCreatorFunction("Convective", NavierStokesAdvection::create);
41  string NavierStokesAdvection::className2 = GetAdvectionTermFactory().RegisterCreatorFunction("NonConservative", NavierStokesAdvection::create);
42 
43  /**
44  * Constructor. Creates ...
45  *
46  * \param
47  * \param
48  */
49 
50  NavierStokesAdvection::NavierStokesAdvection(
51  const LibUtilities::SessionReaderSharedPtr& pSession,
52  const SpatialDomains::MeshGraphSharedPtr& pGraph):
53  AdvectionTerm(pSession, pGraph)
54 
55  {
56 
57  }
58 
59  NavierStokesAdvection::~NavierStokesAdvection()
60  {
61  }
62 
63  //Advection function
64 
65 
66  //Evaluation of the advective terms
67  void NavierStokesAdvection::v_ComputeAdvectionTerm(
68  Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
69  const Array<OneD, Array<OneD, NekDouble> > &pV,
70  const Array<OneD, const NekDouble> &pU,
71  Array<OneD, NekDouble> &pOutarray,
72  int pVelocityComponent,
73  NekDouble m_time,
74  Array<OneD, NekDouble> &pWk)
75  {
76  // use dimension of Velocity vector to dictate dimension of operation
77  int ndim = pV.num_elements();
78  Array<OneD, Array<OneD, NekDouble> > AdvVel (pV.num_elements());
79  Array<OneD, NekDouble> Outarray;
80 
81 
82  int nPointsTot = pFields[0]->GetNpoints();
83  Array<OneD, NekDouble> grad0,grad1,grad2,wkSp;
84 
85  NekDouble OneDptscale = 1.5; // factor to rescale 1d points in dealiasing
86 
87  if(m_specHP_dealiasing)
88  {
89  // Get number of points to dealias a quadratic non-linearity
90  nPointsTot = pFields[0]->Get1DScaledTotPoints(OneDptscale);
91  }
92 
93  grad0 = Array<OneD, NekDouble> (nPointsTot);
94 
95  // interpolate Advection velocity
96  int nadv = pV.num_elements();
97  if(m_specHP_dealiasing) // interpolate advection field to higher space.
98  {
99  AdvVel[0] = Array<OneD, NekDouble> (nPointsTot*(nadv+1));
100  for(int i = 0; i < nadv; ++i)
101  {
102  if(i)
103  {
104  AdvVel[i] = AdvVel[i-1]+nPointsTot;
105  }
106  // interpolate infield to 3/2 dimension
107  pFields[0]->PhysInterp1DScaled(OneDptscale,pV[i],AdvVel[i]);
108  }
109 
110  Outarray = AdvVel[nadv-1] + nPointsTot;
111  }
112  else
113  {
114  for(int i = 0; i < nadv; ++i)
115  {
116  AdvVel[i] = pV[i];
117  }
118 
119  Outarray = pOutarray;
120  }
121 
122  wkSp = Array<OneD, NekDouble> (nPointsTot);
123 
124 
125  // Evaluate V\cdot Grad(u)
126  switch(ndim)
127  {
128  case 1:
129  pFields[0]->PhysDeriv(pU,grad0);
130  Vmath::Vmul(nPointsTot,grad0,1,pV[0],1,pOutarray,1);
131  break;
132  case 2:
133  {
134  grad1 = Array<OneD, NekDouble> (nPointsTot);
135  pFields[0]->PhysDeriv(pU,grad0,grad1);
136 
137  if(m_specHP_dealiasing) // interpolate gradient field
138  {
139  pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
140  Vmath::Vcopy(nPointsTot,wkSp,1,grad0,1);
141  pFields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
142  Vmath::Vcopy(nPointsTot,wkSp,1,grad1,1);
143  }
144 
145  Vmath::Vmul (nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
146  Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,Outarray,1,Outarray,1);
147 
148  if(m_specHP_dealiasing) // Galerkin project solution back to origianl space
149  {
150  pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,pOutarray);
151  }
152 
153  }
154  break;
155  case 3:
156  grad1 = Array<OneD, NekDouble> (pFields[0]->GetNpoints());
157  grad2 = Array<OneD, NekDouble> (pFields[0]->GetNpoints());
158 
159  if(pFields[0]->GetWaveSpace() == false && m_homogen_dealiasing == true )
160  {
161  ASSERTL0(m_specHP_dealiasing == false,"Spectral/hp element dealaising is not set up for this option");
162 
163  pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
164 
165  pFields[0]->DealiasedProd(pV[0],grad0,grad0,m_CoeffState);
166  pFields[0]->DealiasedProd(pV[1],grad1,grad1,m_CoeffState);
167  pFields[0]->DealiasedProd(pV[2],grad2,grad2,m_CoeffState);
168  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,pOutarray,1);
169  Vmath::Vadd(nPointsTot,grad2,1,pOutarray,1,pOutarray,1);
170  }
171  else if(pFields[0]->GetWaveSpace() == true && m_homogen_dealiasing == false)
172  {
173  // take d/dx, d/dy gradients in physical Fourier space
174  pFields[0]->PhysDeriv(pV[pVelocityComponent],grad0,grad1);
175 
176  // Take d/dz derivative using wave space field
177  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],pU,
178  pOutarray);
179  pFields[0]->HomogeneousBwdTrans(pOutarray,grad2);
180 
181  if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
182  {
183  pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
184  Vmath::Vmul(nPointsTot,wkSp,1,AdvVel[0],1,Outarray,1);
185  }
186  else
187  {
188  Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
189  }
190 
191  if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
192  {
193  pFields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
194  Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[1],1,Outarray,1,
195  Outarray,1);
196  }
197  else
198  {
199  Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,Outarray,1,
200  Outarray,1);
201  }
202 
203  if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
204  {
205  pFields[0]->PhysInterp1DScaled(OneDptscale,grad2,wkSp);
206  Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[2],1,Outarray,1,Outarray,1);
207  pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,grad2);
208  pFields[0]->HomogeneousFwdTrans(grad2,pOutarray);
209  }
210  else
211  {
212  Vmath::Vvtvp(nPointsTot,grad2,1,AdvVel[2],1,Outarray,1,grad0,1);
213  pFields[0]->HomogeneousFwdTrans(grad0,pOutarray);
214  }
215  }
216  else if(pFields[0]->GetWaveSpace() == false && m_homogen_dealiasing == false)
217  {
218 
219  pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
220 
221  if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
222  {
223  pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
224  Vmath::Vmul(nPointsTot,wkSp,1,AdvVel[0],1,Outarray,1);
225  }
226  else
227  {
228  Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
229  }
230 
231 
232  if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
233  {
234  pFields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
235  Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[1],1,Outarray,1,
236  Outarray,1);
237  }
238  else
239  {
240  Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,Outarray,1,
241  Outarray,1);
242  }
243 
244  if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
245  {
246  pFields[0]->PhysInterp1DScaled(OneDptscale,grad2,wkSp);
247  Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[2],1,Outarray,1,Outarray,1);
248  pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,pOutarray);
249  }
250  else
251  {
252  Vmath::Vvtvp(nPointsTot,grad2,1,AdvVel[2],1,Outarray,1,pOutarray,1);
253  }
254  }
255  else if(pFields[0]->GetWaveSpace() == true && m_homogen_dealiasing == true)
256  {
257  ASSERTL0(m_specHP_dealiasing == false,"Spectral/hp element dealaising is not set up for this option");
258 
259  pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
260 
261  pFields[0]->HomogeneousBwdTrans(grad0, pOutarray);
262  pFields[0]->DealiasedProd(pV[0], pOutarray, grad0,
263  m_CoeffState);
264 
265  pFields[0]->HomogeneousBwdTrans(grad1,pOutarray);
266  pFields[0]->DealiasedProd(pV[1], pOutarray, grad1,
267  m_CoeffState);
268 
269  pFields[0]->HomogeneousBwdTrans(grad2,pOutarray);
270  pFields[0]->DealiasedProd(pV[2], pOutarray, grad2,
271  m_CoeffState);
272 
273  Vmath::Vadd(nPointsTot, grad0, 1, grad1, 1, grad0, 1);
274  Vmath::Vadd(nPointsTot, grad0, 1, grad2, 1, grad0, 1);
275 
276  pFields[0]->HomogeneousFwdTrans(grad0,pOutarray);
277  }
278  else
279  {
280  ASSERTL0(false, "Advection term calculation not implented or "
281  "possible with the current problem set up");
282  }
283  break;
284  default:
285  ASSERTL0(false,"dimension unknown");
286  }
287  }
288 
289 } //end of namespace
290 
Generated on Sun Mar 15 2015 00:11:33 for Nektar++ by   doxygen 1.8.1.2