Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 using namespace std;
40 
41 namespace Nektar
42 {
43  string UnsteadyAdvection::className = SolverUtils::GetEquationSystemFactory().
44  RegisterCreatorFunction("UnsteadyAdvection",
45  UnsteadyAdvection::create,
46  "Unsteady Advection equation.");
47 
48  UnsteadyAdvection::UnsteadyAdvection(
50  : UnsteadySystem(pSession),
51  AdvectionSystem(pSession)
52  {
53  m_planeNumber = 0;
54  }
55 
56  /**
57  * @brief Initialisation object for the unsteady linear advection equation.
58  */
60  {
61  // Call to the initialisation object of UnsteadySystem
63 
64  m_session->LoadParameter("wavefreq", m_waveFreq, 0.0);
65  // Read the advection velocities from session file
66 
67  std::vector<std::string> vel;
68  vel.push_back("Vx");
69  vel.push_back("Vy");
70  vel.push_back("Vz");
71 
72  // Resize the advection velocities vector to dimension of the problem
73  vel.resize(m_spacedim);
74 
75  // Store in the global variable m_velocity the advection velocities
77  EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
78 
79  // Type of advection class to be used
80  switch(m_projectionType)
81  {
82  // Continuous field
84  {
85  string advName;
86  m_session->LoadSolverInfo(
87  "AdvectionType", advName, "NonConservative");
89  GetAdvectionFactory().CreateInstance(advName, advName);
91  {
92  m_advObject->SetFluxVector(
94  }
95  else
96  {
97  m_advObject->SetFluxVector(
99  }
100  break;
101  }
102  // Discontinuous field
104  {
105  // Do not forwards transform initial condition
106  m_homoInitialFwd = false;
107 
108  // Define the normal velocity fields
109  if (m_fields[0]->GetTrace())
110  {
112  }
113 
114  string advName;
115  string riemName;
116  m_session->LoadSolverInfo(
117  "AdvectionType", advName, "WeakDG");
119  GetAdvectionFactory().CreateInstance(advName, advName);
121  {
122  m_advObject->SetFluxVector(
124  }
125  else
126  {
127  m_advObject->SetFluxVector(
129  }
130  m_session->LoadSolverInfo(
131  "UpwindType", riemName, "Upwind");
134  m_riemannSolver->SetScalar(
136 
137  m_advObject->SetRiemannSolver(m_riemannSolver);
138  m_advObject->InitObject(m_session, m_fields);
139  break;
140  }
141  default:
142  {
143  ASSERTL0(false, "Unsupported projection type.");
144  break;
145  }
146  }
147 
148  // If explicit it computes RHS and PROJECTION for the time integration
150  {
153  }
154  // Otherwise it gives an error (no implicit integration)
155  else
156  {
157  ASSERTL0(false, "Implicit unsteady Advection not set up.");
158  }
159  }
160 
161  /**
162  * @brief Unsteady linear advection equation destructor.
163  */
165  {
166  }
167 
168  /**
169  * @brief Get the normal velocity for the linear advection equation.
170  */
172  {
173  // Number of trace (interface) points
174  int i;
175  int nTracePts = GetTraceNpoints();
176 
177  // Auxiliary variable to compute the normal velocity
178  Array<OneD, NekDouble> tmp(nTracePts);
179 
180  // Reset the normal velocity
181  Vmath::Zero(nTracePts, m_traceVn, 1);
182 
183  for (i = 0; i < m_velocity.num_elements(); ++i)
184  {
185  m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
186 
187  Vmath::Vvtvp(nTracePts,
188  m_traceNormals[i], 1,
189  tmp, 1,
190  m_traceVn, 1,
191  m_traceVn, 1);
192  }
193 
194  return m_traceVn;
195  }
196 
197  /**
198  * @brief Compute the right-hand side for the linear advection equation.
199  *
200  * @param inarray Given fields.
201  * @param outarray Calculated solution.
202  * @param time Time.
203  */
205  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
206  Array<OneD, Array<OneD, NekDouble> >&outarray,
207  const NekDouble time)
208  {
209  // Counter variable
210  int i;
211 
212  // Number of fields (variables of the problem)
213  int nVariables = inarray.num_elements();
214 
215  // Number of solution points
216  int nSolutionPts = GetNpoints();
217 
218  // RHS computation using the new advection base class
219  m_advObject->Advect(nVariables, m_fields, m_velocity, inarray,
220  outarray, time);
221 
222  // Negate the RHS
223  for (i = 0; i < nVariables; ++i)
224  {
225  Vmath::Neg(nSolutionPts, outarray[i], 1);
226  }
227  }
228 
229  /**
230  * @brief Compute the projection for the linear advection equation.
231  *
232  * @param inarray Given fields.
233  * @param outarray Calculated solution.
234  * @param time Time.
235  */
237  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
238  Array<OneD, Array<OneD, NekDouble> >&outarray,
239  const NekDouble time)
240  {
241  // Counter variable
242  int i;
243 
244  // Number of fields (variables of the problem)
245  int nVariables = inarray.num_elements();
246 
247  // Set the boundary conditions
248  SetBoundaryConditions(time);
249 
250  // Switch on the projection type (Discontinuous or Continuous)
251  switch(m_projectionType)
252  {
253  // Discontinuous projection
255  {
256  // Number of quadrature points
257  int nQuadraturePts = GetNpoints();
258 
259  // Just copy over array
260  for(i = 0; i < nVariables; ++i)
261  {
262  Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
263  }
264  break;
265  }
266 
267  // Continuous projection
270  {
271  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs(),0.0);
272  for(i = 0; i < nVariables; ++i)
273  {
274  m_fields[i]->FwdTrans(inarray[i], coeffs);
275  m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
276  }
277  break;
278  }
279 
280  default:
281  ASSERTL0(false,"Unknown projection scheme");
282  break;
283  }
284  }
285 
286  /**
287  * @brief Return the flux vector for the linear advection equation.
288  *
289  * @param i Component of the flux vector to calculate.
290  * @param physfield Fields.
291  * @param flux Resulting flux.
292  */
294  const Array<OneD, Array<OneD, NekDouble> > &physfield,
296  {
297  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
298  "Dimension of flux array and velocity array do not match");
299 
300  int i , j;
301  int nq = physfield[0].num_elements();
302 
303  for (i = 0; i < flux.num_elements(); ++i)
304  {
305  for (j = 0; j < flux[0].num_elements(); ++j)
306  {
307  Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
308  flux[i][j], 1);
309  }
310  }
311  }
312 
313  /**
314  * @brief Return the flux vector for the linear advection equation using
315  * the dealiasing technique.
316  *
317  * @param i Component of the flux vector to calculate.
318  * @param physfield Fields.
319  * @param flux Resulting flux.
320  */
322  const Array<OneD, Array<OneD, NekDouble> > &physfield,
324  {
325  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
326  "Dimension of flux array and velocity array do not match");
327 
328  int i, j;
329  int nq = physfield[0].num_elements();
330  int nVariables = physfield.num_elements();
331 
332  // Factor to rescale 1d points in dealiasing
333  NekDouble OneDptscale = 2;
334 
336  advVel_plane(m_velocity.num_elements());
337 
338  // Get number of points to dealias a cubic non-linearity
339  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
340 
341  // Initialisation of higher-space variables
342  Array<OneD, Array<OneD, NekDouble> >physfieldInterp(nVariables);
344  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >fluxInterp(nVariables);
345 
346  // Interpolation to higher space of physfield
347  for (i = 0; i < nVariables; ++i)
348  {
349  physfieldInterp[i] = Array<OneD, NekDouble>(nq);
350  fluxInterp[i] = Array<OneD, Array<OneD, NekDouble> >(m_expdim);
351  for (j = 0; j < m_expdim; ++j)
352  {
353  fluxInterp[i][j] = Array<OneD, NekDouble>(nq);
354  }
355 
356  m_fields[0]->PhysInterp1DScaled(
357  OneDptscale, physfield[i], physfieldInterp[i]);
358  }
359 
360  // Interpolation to higher space of velocity
361  for (j = 0; j < m_expdim; ++j)
362  {
363  velocityInterp[j] = Array<OneD, NekDouble>(nq);
364 
365  m_fields[0]->PhysInterp1DScaled(
366  OneDptscale, m_velocity[j], velocityInterp[j]);
367  }
368 
369  // Evaluation of flux vector in the higher space
370  for (i = 0; i < flux.num_elements(); ++i)
371  {
372  for (j = 0; j < flux[0].num_elements(); ++j)
373  {
374  Vmath::Vmul(nq, physfieldInterp[i], 1, velocityInterp[j], 1,
375  fluxInterp[i][j], 1);
376  }
377  }
378 
379  // Galerkin project solution back to original space
380  for (i = 0; i < nVariables; ++i)
381  {
382  for (j = 0; j < m_spacedim; ++j)
383  {
384  m_fields[0]->PhysGalerkinProjection1DScaled(
385  OneDptscale, fluxInterp[i][j], flux[i][j]);
386  }
387  }
388  }
389 
391  {
393  }
394 }
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
virtual void v_InitObject()
Initialise the object.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point.
int m_expdim
Expansion dimension.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
STL namespace.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
Base class for unsteady solvers.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
RiemannSolverFactory & GetRiemannSolverFactory()
int m_spacedim
Spatial dimension (>= expansion dim).
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
EquationSystemFactory & GetEquationSystemFactory()
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
void GetFluxVectorDeAlias(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point using dealiasing.
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, NekDouble > m_traceVn
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
virtual ~UnsteadyAdvection()
Destructor.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
A base class for PDEs which include an advection component.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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