Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UnsteadyAdvection.cpp
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////
2 //
3 // File UnsteadyAdvection.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: Unsteady linear advection solve routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
38 
39 namespace Nektar
40 {
42  RegisterCreatorFunction("UnsteadyAdvection",
44  "Unsteady Advection equation.");
45 
48  : UnsteadySystem(pSession)
49  {
50  m_planeNumber = 0;
51  }
52 
53  /**
54  * @brief Initialisation object for the unsteady linear advection equation.
55  */
57  {
58  // Call to the initialisation object of UnsteadySystem
60 
61  m_session->LoadParameter("wavefreq", m_waveFreq, 0.0);
62  // Read the advection velocities from session file
63 
64  std::vector<std::string> vel;
65  vel.push_back("Vx");
66  vel.push_back("Vy");
67  vel.push_back("Vz");
68 
69  // Resize the advection velocities vector to dimension of the problem
70  vel.resize(m_spacedim);
71 
72  // Store in the global variable m_velocity the advection velocities
73  m_velocity = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
74  EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
75 
76  // Type of advection class to be used
77  switch(m_projectionType)
78  {
79  // Continuous field
81  {
82  string advName;
83  m_session->LoadSolverInfo(
84  "AdvectionType", advName, "NonConservative");
86  GetAdvectionFactory().CreateInstance(advName, advName);
88  {
89  m_advection->SetFluxVector(
91  }
92  else
93  {
94  m_advection->SetFluxVector(
96  }
97  break;
98  }
99  // Discontinuous field
101  {
102  // Do not forwards transform initial condition
103  m_homoInitialFwd = false;
104 
105  // Define the normal velocity fields
106  if (m_fields[0]->GetTrace())
107  {
108  m_traceVn = Array<OneD, NekDouble>(GetTraceNpoints());
109  }
110 
111  string advName;
112  string riemName;
113  m_session->LoadSolverInfo(
114  "AdvectionType", advName, "WeakDG");
116  GetAdvectionFactory().CreateInstance(advName, advName);
118  {
119  m_advection->SetFluxVector(
121  }
122  else
123  {
124  m_advection->SetFluxVector(
126  }
127  m_session->LoadSolverInfo(
128  "UpwindType", riemName, "Upwind");
131  m_riemannSolver->SetScalar(
133 
134  m_advection->SetRiemannSolver(m_riemannSolver);
135  m_advection->InitObject(m_session, m_fields);
136  break;
137  }
138  default:
139  {
140  ASSERTL0(false, "Unsupported projection type.");
141  break;
142  }
143  }
144 
145  // If explicit it computes RHS and PROJECTION for the time integration
147  {
150  }
151  // Otherwise it gives an error (no implicit integration)
152  else
153  {
154  ASSERTL0(false, "Implicit unsteady Advection not set up.");
155  }
156  }
157 
158  /**
159  * @brief Unsteady linear advection equation destructor.
160  */
162  {
163  }
164 
165  /**
166  * @brief Get the normal velocity for the linear advection equation.
167  */
168  Array<OneD, NekDouble> &UnsteadyAdvection::GetNormalVelocity()
169  {
170  // Number of trace (interface) points
171  int i;
172  int nTracePts = GetTraceNpoints();
173 
174  // Auxiliary variable to compute the normal velocity
175  Array<OneD, NekDouble> tmp(nTracePts);
176 
177  // Reset the normal velocity
178  Vmath::Zero(nTracePts, m_traceVn, 1);
179 
180  for (i = 0; i < m_velocity.num_elements(); ++i)
181  {
182  m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
183 
184  Vmath::Vvtvp(nTracePts,
185  m_traceNormals[i], 1,
186  tmp, 1,
187  m_traceVn, 1,
188  m_traceVn, 1);
189  }
190 
191  return m_traceVn;
192  }
193 
194  /**
195  * @brief Compute the right-hand side for the linear advection equation.
196  *
197  * @param inarray Given fields.
198  * @param outarray Calculated solution.
199  * @param time Time.
200  */
202  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
203  Array<OneD, Array<OneD, NekDouble> >&outarray,
204  const NekDouble time)
205  {
206  // Counter variable
207  int i;
208 
209  // Number of fields (variables of the problem)
210  int nVariables = inarray.num_elements();
211 
212  // Number of solution points
213  int nSolutionPts = GetNpoints();
214 
215  // RHS computation using the new advection base class
216  m_advection->Advect(nVariables, m_fields, m_velocity, inarray,
217  outarray);
218 
219  // Negate the RHS
220  for (i = 0; i < nVariables; ++i)
221  {
222  Vmath::Neg(nSolutionPts, outarray[i], 1);
223  }
224  }
225 
226  /**
227  * @brief Compute the projection for the linear advection equation.
228  *
229  * @param inarray Given fields.
230  * @param outarray Calculated solution.
231  * @param time Time.
232  */
234  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
235  Array<OneD, Array<OneD, NekDouble> >&outarray,
236  const NekDouble time)
237  {
238  // Counter variable
239  int i;
240 
241  // Number of fields (variables of the problem)
242  int nVariables = inarray.num_elements();
243 
244  // Set the boundary conditions
245  SetBoundaryConditions(time);
246 
247  // Switch on the projection type (Discontinuous or Continuous)
248  switch(m_projectionType)
249  {
250  // Discontinuous projection
252  {
253  // Number of quadrature points
254  int nQuadraturePts = GetNpoints();
255 
256  // Just copy over array
257  for(i = 0; i < nVariables; ++i)
258  {
259  Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
260  }
261  break;
262  }
263 
264  // Continuous projection
267  {
268  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs(),0.0);
269  for(i = 0; i < nVariables; ++i)
270  {
271  m_fields[i]->FwdTrans(inarray[i], coeffs);
272  m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
273  }
274  break;
275  }
276 
277  default:
278  ASSERTL0(false,"Unknown projection scheme");
279  break;
280  }
281  }
282 
283  /**
284  * @brief Return the flux vector for the linear advection equation.
285  *
286  * @param i Component of the flux vector to calculate.
287  * @param physfield Fields.
288  * @param flux Resulting flux.
289  */
291  const Array<OneD, Array<OneD, NekDouble> > &physfield,
292  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
293  {
294  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
295  "Dimension of flux array and velocity array do not match");
296 
297  int i , j;
298  int nq = physfield[0].num_elements();
299 
300  for (i = 0; i < flux.num_elements(); ++i)
301  {
302  for (j = 0; j < flux[0].num_elements(); ++j)
303  {
304  Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
305  flux[i][j], 1);
306  }
307  }
308  }
309 
310  /**
311  * @brief Return the flux vector for the linear advection equation using
312  * the dealiasing technique.
313  *
314  * @param i Component of the flux vector to calculate.
315  * @param physfield Fields.
316  * @param flux Resulting flux.
317  */
319  const Array<OneD, Array<OneD, NekDouble> > &physfield,
320  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
321  {
322  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
323  "Dimension of flux array and velocity array do not match");
324 
325  int i, j;
326  int nq = physfield[0].num_elements();
327  int nVariables = physfield.num_elements();
328 
329  // Factor to rescale 1d points in dealiasing
330  NekDouble OneDptscale = 2;
331 
332  Array<OneD, Array<OneD, NekDouble> >
333  advVel_plane(m_velocity.num_elements());
334 
335  // Get number of points to dealias a cubic non-linearity
336  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
337 
338  // Initialisation of higher-space variables
339  Array<OneD, Array<OneD, NekDouble> >physfieldInterp(nVariables);
340  Array<OneD, Array<OneD, NekDouble> >velocityInterp(m_expdim);
341  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >fluxInterp(nVariables);
342 
343  // Interpolation to higher space of physfield
344  for (i = 0; i < nVariables; ++i)
345  {
346  physfieldInterp[i] = Array<OneD, NekDouble>(nq);
347  fluxInterp[i] = Array<OneD, Array<OneD, NekDouble> >(m_expdim);
348  for (j = 0; j < m_expdim; ++j)
349  {
350  fluxInterp[i][j] = Array<OneD, NekDouble>(nq);
351  }
352 
353  m_fields[0]->PhysInterp1DScaled(
354  OneDptscale, physfield[i], physfieldInterp[i]);
355  }
356 
357  // Interpolation to higher space of velocity
358  for (j = 0; j < m_expdim; ++j)
359  {
360  velocityInterp[j] = Array<OneD, NekDouble>(nq);
361 
362  m_fields[0]->PhysInterp1DScaled(
363  OneDptscale, m_velocity[j], velocityInterp[j]);
364  }
365 
366  // Evaluation of flux vector in the higher space
367  for (i = 0; i < flux.num_elements(); ++i)
368  {
369  for (j = 0; j < flux[0].num_elements(); ++j)
370  {
371  Vmath::Vmul(nq, physfieldInterp[i], 1, velocityInterp[j], 1,
372  fluxInterp[i][j], 1);
373  }
374  }
375 
376  // Galerkin project solution back to original space
377  for (i = 0; i < nVariables; ++i)
378  {
379  for (j = 0; j < m_spacedim; ++j)
380  {
381  m_fields[0]->PhysGalerkinProjection1DScaled(
382  OneDptscale, fluxInterp[i][j], flux[i][j]);
383  }
384  }
385  }
386 
388  {
390  }
391 }