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  AdvectionSystem(pSession)
50  {
51  m_planeNumber = 0;
52  }
53 
54  /**
55  * @brief Initialisation object for the unsteady linear advection equation.
56  */
58  {
59  // Call to the initialisation object of UnsteadySystem
60  AdvectionSystem::v_InitObject();
61 
62  m_session->LoadParameter("wavefreq", m_waveFreq, 0.0);
63  // Read the advection velocities from session file
64 
65  std::vector<std::string> vel;
66  vel.push_back("Vx");
67  vel.push_back("Vy");
68  vel.push_back("Vz");
69 
70  // Resize the advection velocities vector to dimension of the problem
71  vel.resize(m_spacedim);
72 
73  // Store in the global variable m_velocity the advection velocities
75  EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
76 
77  // Type of advection class to be used
78  switch(m_projectionType)
79  {
80  // Continuous field
82  {
83  string advName;
84  m_session->LoadSolverInfo(
85  "AdvectionType", advName, "NonConservative");
87  GetAdvectionFactory().CreateInstance(advName, advName);
89  {
90  m_advObject->SetFluxVector(
92  }
93  else
94  {
95  m_advObject->SetFluxVector(
97  }
98  break;
99  }
100  // Discontinuous field
102  {
103  // Do not forwards transform initial condition
104  m_homoInitialFwd = false;
105 
106  // Define the normal velocity fields
107  if (m_fields[0]->GetTrace())
108  {
110  }
111 
112  string advName;
113  string riemName;
114  m_session->LoadSolverInfo(
115  "AdvectionType", advName, "WeakDG");
117  GetAdvectionFactory().CreateInstance(advName, advName);
119  {
120  m_advObject->SetFluxVector(
122  }
123  else
124  {
125  m_advObject->SetFluxVector(
127  }
128  m_session->LoadSolverInfo(
129  "UpwindType", riemName, "Upwind");
132  m_riemannSolver->SetScalar(
134 
135  m_advObject->SetRiemannSolver(m_riemannSolver);
136  m_advObject->InitObject(m_session, m_fields);
137  break;
138  }
139  default:
140  {
141  ASSERTL0(false, "Unsupported projection type.");
142  break;
143  }
144  }
145 
146  // If explicit it computes RHS and PROJECTION for the time integration
148  {
151  }
152  // Otherwise it gives an error (no implicit integration)
153  else
154  {
155  ASSERTL0(false, "Implicit unsteady Advection not set up.");
156  }
157  }
158 
159  /**
160  * @brief Unsteady linear advection equation destructor.
161  */
163  {
164  }
165 
166  /**
167  * @brief Get the normal velocity for the linear advection equation.
168  */
170  {
171  // Number of trace (interface) points
172  int i;
173  int nTracePts = GetTraceNpoints();
174 
175  // Auxiliary variable to compute the normal velocity
176  Array<OneD, NekDouble> tmp(nTracePts);
177 
178  // Reset the normal velocity
179  Vmath::Zero(nTracePts, m_traceVn, 1);
180 
181  for (i = 0; i < m_velocity.num_elements(); ++i)
182  {
183  m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
184 
185  Vmath::Vvtvp(nTracePts,
186  m_traceNormals[i], 1,
187  tmp, 1,
188  m_traceVn, 1,
189  m_traceVn, 1);
190  }
191 
192  return m_traceVn;
193  }
194 
195  /**
196  * @brief Compute the right-hand side for the linear advection equation.
197  *
198  * @param inarray Given fields.
199  * @param outarray Calculated solution.
200  * @param time Time.
201  */
203  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
204  Array<OneD, Array<OneD, NekDouble> >&outarray,
205  const NekDouble time)
206  {
207  // Counter variable
208  int i;
209 
210  // Number of fields (variables of the problem)
211  int nVariables = inarray.num_elements();
212 
213  // Number of solution points
214  int nSolutionPts = GetNpoints();
215 
216  // RHS computation using the new advection base class
217  m_advObject->Advect(nVariables, m_fields, m_velocity, inarray,
218  outarray, time);
219 
220  // Negate the RHS
221  for (i = 0; i < nVariables; ++i)
222  {
223  Vmath::Neg(nSolutionPts, outarray[i], 1);
224  }
225  }
226 
227  /**
228  * @brief Compute the projection for the linear advection equation.
229  *
230  * @param inarray Given fields.
231  * @param outarray Calculated solution.
232  * @param time Time.
233  */
235  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
236  Array<OneD, Array<OneD, NekDouble> >&outarray,
237  const NekDouble time)
238  {
239  // Counter variable
240  int i;
241 
242  // Number of fields (variables of the problem)
243  int nVariables = inarray.num_elements();
244 
245  // Set the boundary conditions
246  SetBoundaryConditions(time);
247 
248  // Switch on the projection type (Discontinuous or Continuous)
249  switch(m_projectionType)
250  {
251  // Discontinuous projection
253  {
254  // Number of quadrature points
255  int nQuadraturePts = GetNpoints();
256 
257  // Just copy over array
258  for(i = 0; i < nVariables; ++i)
259  {
260  Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
261  }
262  break;
263  }
264 
265  // Continuous projection
268  {
269  Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs(),0.0);
270  for(i = 0; i < nVariables; ++i)
271  {
272  m_fields[i]->FwdTrans(inarray[i], coeffs);
273  m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
274  }
275  break;
276  }
277 
278  default:
279  ASSERTL0(false,"Unknown projection scheme");
280  break;
281  }
282  }
283 
284  /**
285  * @brief Return the flux vector for the linear advection equation.
286  *
287  * @param i Component of the flux vector to calculate.
288  * @param physfield Fields.
289  * @param flux Resulting flux.
290  */
292  const Array<OneD, Array<OneD, NekDouble> > &physfield,
294  {
295  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
296  "Dimension of flux array and velocity array do not match");
297 
298  int i , j;
299  int nq = physfield[0].num_elements();
300 
301  for (i = 0; i < flux.num_elements(); ++i)
302  {
303  for (j = 0; j < flux[0].num_elements(); ++j)
304  {
305  Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
306  flux[i][j], 1);
307  }
308  }
309  }
310 
311  /**
312  * @brief Return the flux vector for the linear advection equation using
313  * the dealiasing technique.
314  *
315  * @param i Component of the flux vector to calculate.
316  * @param physfield Fields.
317  * @param flux Resulting flux.
318  */
320  const Array<OneD, Array<OneD, NekDouble> > &physfield,
322  {
323  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
324  "Dimension of flux array and velocity array do not match");
325 
326  int i, j;
327  int nq = physfield[0].num_elements();
328  int nVariables = physfield.num_elements();
329 
330  // Factor to rescale 1d points in dealiasing
331  NekDouble OneDptscale = 2;
332 
334  advVel_plane(m_velocity.num_elements());
335 
336  // Get number of points to dealias a cubic non-linearity
337  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
338 
339  // Initialisation of higher-space variables
340  Array<OneD, Array<OneD, NekDouble> >physfieldInterp(nVariables);
342  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >fluxInterp(nVariables);
343 
344  // Interpolation to higher space of physfield
345  for (i = 0; i < nVariables; ++i)
346  {
347  physfieldInterp[i] = Array<OneD, NekDouble>(nq);
348  fluxInterp[i] = Array<OneD, Array<OneD, NekDouble> >(m_expdim);
349  for (j = 0; j < m_expdim; ++j)
350  {
351  fluxInterp[i][j] = Array<OneD, NekDouble>(nq);
352  }
353 
354  m_fields[0]->PhysInterp1DScaled(
355  OneDptscale, physfield[i], physfieldInterp[i]);
356  }
357 
358  // Interpolation to higher space of velocity
359  for (j = 0; j < m_expdim; ++j)
360  {
361  velocityInterp[j] = Array<OneD, NekDouble>(nq);
362 
363  m_fields[0]->PhysInterp1DScaled(
364  OneDptscale, m_velocity[j], velocityInterp[j]);
365  }
366 
367  // Evaluation of flux vector in the higher space
368  for (i = 0; i < flux.num_elements(); ++i)
369  {
370  for (j = 0; j < flux[0].num_elements(); ++j)
371  {
372  Vmath::Vmul(nq, physfieldInterp[i], 1, velocityInterp[j], 1,
373  fluxInterp[i][j], 1);
374  }
375  }
376 
377  // Galerkin project solution back to original space
378  for (i = 0; i < nVariables; ++i)
379  {
380  for (j = 0; j < m_spacedim; ++j)
381  {
382  m_fields[0]->PhysGalerkinProjection1DScaled(
383  OneDptscale, fluxInterp[i][j], flux[i][j]);
384  }
385  }
386  }
387 
389  {
390  AdvectionSystem::v_GenerateSummary(s);
391  }
392 }
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:161
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:428
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
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)
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.
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).
UnsteadyAdvection(const LibUtilities::SessionReaderSharedPtr &pSession)
Session reader.
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:382
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.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
static std::string className
Name of class.
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