Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UnsteadyInviscidBurger.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File UnsteadyInviscidBurger.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 inviscid Burger solve routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42  string UnsteadyInviscidBurger::className
44  "UnsteadyInviscidBurger",
45  UnsteadyInviscidBurger::create,
46  "Inviscid Burger equation");
47 
48  UnsteadyInviscidBurger::UnsteadyInviscidBurger(
50  : UnsteadySystem(pSession),
51  AdvectionSystem(pSession)
52  {
53  }
54 
55  /**
56  * @brief Initialisation object for the inviscid Burger equation.
57  */
59  {
60  // Call to the initialisation object of UnsteadySystem
62 
63  // Define the normal velocity fields
64  if (m_fields[0]->GetTrace())
65  {
67  }
68 
69  // Type of advection class to be used
70  switch(m_projectionType)
71  {
72  // Continuous field
74  {
75  string advName;
76  m_session->LoadSolverInfo("AdvectionType", advName, "NonConservative");
78  m_advObject->SetFluxVector (&UnsteadyInviscidBurger::GetFluxVector, this);
79  break;
80  }
81  // Discontinuous field
83  {
84  string advName;
85  string riemName;
86 
87  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
89  m_advObject->SetFluxVector (&UnsteadyInviscidBurger::GetFluxVector, this);
90 
91  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
94 
95  m_advObject->SetRiemannSolver(m_riemannSolver);
96  m_advObject->InitObject (m_session, m_fields);
97  break;
98  }
99  default:
100  {
101  ASSERTL0(false, "Unsupported projection type.");
102  break;
103  }
104  }
105 
106  // If explicit it computes RHS and PROJECTION for the time integration
108  {
111  }
112  // Otherwise it gives an error because (no implicit integration)
113  else
114  {
115  ASSERTL0(false, "Implicit unsteady Advection not set up.");
116  }
117  }
118 
119  /**
120  * @brief Inviscid Burger equation destructor.
121  */
123  {
124  }
125 
126  /**
127  * @brief Get the normal velocity for the inviscid Burger equation.
128  */
130  {
131  // Counter variable
132  int i;
133 
134  // Number of trace (interface) points
135  int nTracePts = GetTraceNpoints();
136 
137  // Number of solution points
138  int nSolutionPts = GetNpoints();
139 
140  // Number of fields (variables of the problem)
141  int nVariables = m_fields.num_elements();
142 
143  // Auxiliary variables to compute the normal velocity
144  Array<OneD, NekDouble> Fwd (nTracePts);
145  Array<OneD, Array<OneD, NekDouble> > physfield (nVariables);
146 
147  // Reset the normal velocity
148  Vmath::Zero(nTracePts, m_traceVn, 1);
149 
150  // The TimeIntegration Class does not update the physical values of the
151  // solution. It is thus necessary to transform back the coefficient into
152  // the physical space and save them in physfield to compute the normal
153  // advection velocity properly. However it remains a critical point.
154  for(i = 0; i < nVariables; ++i)
155  {
156  physfield[i] = Array<OneD, NekDouble>(nSolutionPts);
157  m_fields[i]->BwdTrans_IterPerExp(m_fields[i]->GetCoeffs(),
158  physfield[i]);
159  }
160 
161  /// Extract the physical values at the trace space
162  m_fields[0]->ExtractTracePhys(physfield[0], Fwd);
163 
164  /// Compute the normal velocity
165  for (i = 0; i < m_spacedim; ++i)
166  {
167  Vmath::Vvtvp(nTracePts,
168  m_traceNormals[i], 1,
169  Fwd, 1,
170  m_traceVn, 1,
171  m_traceVn, 1);
172 
173  Vmath::Smul(nTracePts, 0.5, m_traceVn, 1, m_traceVn, 1);
174  }
175  return m_traceVn;
176  }
177 
178  /**
179  * @brief Compute the right-hand side for the inviscid Burger equation.
180  *
181  * @param inarray Given fields.
182  * @param outarray Calculated solution.
183  * @param time Time.
184  */
186  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
187  Array<OneD, Array<OneD, NekDouble> >&outarray,
188  const NekDouble time)
189  {
190  // Counter variable
191  int i;
192 
193  // Number of fields (variables of the problem)
194  int nVariables = inarray.num_elements();
195 
196  // Number of solution points
197  int nSolutionPts = GetNpoints();
198 
199  // !Useless variable for WeakDG and FR!
201 
202  // RHS computation using the new advection base class
203  m_advObject->Advect(nVariables, m_fields, advVel, inarray,
204  outarray, time);
205 
206  // Negate the RHS
207  for (i = 0; i < nVariables; ++i)
208  {
209  Vmath::Neg(nSolutionPts, outarray[i], 1);
210  }
211  }
212 
213  /**
214  * @brief Compute the projection for the inviscid Burger equation.
215  *
216  * @param inarray Given fields.
217  * @param outarray Calculated solution.
218  * @param time Time.
219  */
221  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
222  Array<OneD, Array<OneD, NekDouble> >&outarray,
223  const NekDouble time)
224  {
225  // Counter variable
226  int i;
227 
228  // Number of variables of the problem
229  int nVariables = inarray.num_elements();
230 
231  // Set the boundary conditions
232  SetBoundaryConditions(time);
233 
234  // Switch on the projection type (Discontinuous or Continuous)
235  switch(m_projectionType)
236  {
237  // Discontinuous projection
239  {
240  // Number of quadrature points
241  int nQuadraturePts = GetNpoints();
242 
243  // Just copy over array
244  for(i = 0; i < nVariables; ++i)
245  {
246  Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
247  }
248  break;
249  }
250 
251  // Continuous projection
254  {
256 
257  for(i = 0; i < nVariables; ++i)
258  {
259  m_fields[i]->FwdTrans(inarray[i], coeffs);
260  m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
261  }
262  break;
263  }
264  default:
265  ASSERTL0(false, "Unknown projection scheme");
266  break;
267  }
268  }
269 
270  /**
271  * @brief Return the flux vector for the inviscid Burger equation.
272  *
273  * @param i Component of the flux vector to calculate.
274  * @param physfield Fields.
275  * @param flux Resulting flux.
276  */
278  const Array<OneD, Array<OneD, NekDouble> > &physfield,
280  {
281  const int nq = GetNpoints();
282 
283  for (int i = 0; i < flux.num_elements(); ++i)
284  {
285  for (int j = 0; j < flux[0].num_elements(); ++j)
286  {
287  Vmath::Vmul(nq, physfield[i], 1, physfield[i], 1,
288  flux[i][j], 1);
289  Vmath::Smul(nq, 0.5, flux[i][j], 1, flux[i][j], 1);
290  }
291  }
292  }
293 }
294 
Array< OneD, NekDouble > m_traceVn
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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.
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.
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
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
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
RiemannSolverFactory & GetRiemannSolverFactory()
int m_spacedim
Spatial dimension (>= expansion dim).
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
virtual ~UnsteadyInviscidBurger()
Destructor.
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
virtual void v_InitObject()
Initialise the object.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
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:359
A base class for PDEs which include an advection component.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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