Nektar++
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 namespace Nektar
39 {
42  "UnsteadyInviscidBurger",
44  "Inviscid Burger equation");
45 
48  : UnsteadySystem(pSession),
49  AdvectionSystem(pSession)
50  {
51  }
52 
53  /**
54  * @brief Initialisation object for the inviscid Burger equation.
55  */
57  {
58  // Call to the initialisation object of UnsteadySystem
59  AdvectionSystem::v_InitObject();
60 
61  // Define the normal velocity fields
62  if (m_fields[0]->GetTrace())
63  {
65  }
66 
67  // Type of advection class to be used
68  switch(m_projectionType)
69  {
70  // Continuous field
72  {
73  string advName;
74  m_session->LoadSolverInfo("AdvectionType", advName, "NonConservative");
76  m_advObject->SetFluxVector (&UnsteadyInviscidBurger::GetFluxVector, this);
77  break;
78  }
79  // Discontinuous field
81  {
82  string advName;
83  string riemName;
84 
85  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
87  m_advObject->SetFluxVector (&UnsteadyInviscidBurger::GetFluxVector, this);
88 
89  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
92 
93  m_advObject->SetRiemannSolver(m_riemannSolver);
94  m_advObject->InitObject (m_session, m_fields);
95  break;
96  }
97  default:
98  {
99  ASSERTL0(false, "Unsupported projection type.");
100  break;
101  }
102  }
103 
104  // If explicit it computes RHS and PROJECTION for the time integration
106  {
109  }
110  // Otherwise it gives an error because (no implicit integration)
111  else
112  {
113  ASSERTL0(false, "Implicit unsteady Advection not set up.");
114  }
115  }
116 
117  /**
118  * @brief Inviscid Burger equation destructor.
119  */
121  {
122  }
123 
124  /**
125  * @brief Get the normal velocity for the inviscid Burger equation.
126  */
128  {
129  // Counter variable
130  int i;
131 
132  // Number of trace (interface) points
133  int nTracePts = GetTraceNpoints();
134 
135  // Number of solution points
136  int nSolutionPts = GetNpoints();
137 
138  // Number of fields (variables of the problem)
139  int nVariables = m_fields.num_elements();
140 
141  // Auxiliary variables to compute the normal velocity
142  Array<OneD, NekDouble> Fwd (nTracePts);
143  Array<OneD, Array<OneD, NekDouble> > physfield (nVariables);
144 
145  // Reset the normal velocity
146  Vmath::Zero(nTracePts, m_traceVn, 1);
147 
148  // The TimeIntegration Class does not update the physical values of the
149  // solution. It is thus necessary to transform back the coefficient into
150  // the physical space and save them in physfield to compute the normal
151  // advection velocity properly. However it remains a critical point.
152  for(i = 0; i < nVariables; ++i)
153  {
154  physfield[i] = Array<OneD, NekDouble>(nSolutionPts);
155  m_fields[i]->BwdTrans_IterPerExp(m_fields[i]->GetCoeffs(),
156  physfield[i]);
157  }
158 
159  /// Extract the physical values at the trace space
160  m_fields[0]->ExtractTracePhys(physfield[0], Fwd);
161 
162  /// Compute the normal velocity
163  for (i = 0; i < m_spacedim; ++i)
164  {
165  Vmath::Vvtvp(nTracePts,
166  m_traceNormals[i], 1,
167  Fwd, 1,
168  m_traceVn, 1,
169  m_traceVn, 1);
170 
171  Vmath::Smul(nTracePts, 0.5, m_traceVn, 1, m_traceVn, 1);
172  }
173  return m_traceVn;
174  }
175 
176  /**
177  * @brief Compute the right-hand side for the inviscid Burger equation.
178  *
179  * @param inarray Given fields.
180  * @param outarray Calculated solution.
181  * @param time Time.
182  */
184  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
185  Array<OneD, Array<OneD, NekDouble> >&outarray,
186  const NekDouble time)
187  {
188  // Counter variable
189  int i;
190 
191  // Number of fields (variables of the problem)
192  int nVariables = inarray.num_elements();
193 
194  // Number of solution points
195  int nSolutionPts = GetNpoints();
196 
197  // !Useless variable for WeakDG and FR!
199 
200  // RHS computation using the new advection base class
201  m_advObject->Advect(nVariables, m_fields, advVel, inarray,
202  outarray, time);
203 
204  // Negate the RHS
205  for (i = 0; i < nVariables; ++i)
206  {
207  Vmath::Neg(nSolutionPts, outarray[i], 1);
208  }
209  }
210 
211  /**
212  * @brief Compute the projection for the inviscid Burger equation.
213  *
214  * @param inarray Given fields.
215  * @param outarray Calculated solution.
216  * @param time Time.
217  */
219  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
220  Array<OneD, Array<OneD, NekDouble> >&outarray,
221  const NekDouble time)
222  {
223  // Counter variable
224  int i;
225 
226  // Number of variables of the problem
227  int nVariables = inarray.num_elements();
228 
229  // Set the boundary conditions
230  SetBoundaryConditions(time);
231 
232  // Switch on the projection type (Discontinuous or Continuous)
233  switch(m_projectionType)
234  {
235  // Discontinuous projection
237  {
238  // Number of quadrature points
239  int nQuadraturePts = GetNpoints();
240 
241  // Just copy over array
242  for(i = 0; i < nVariables; ++i)
243  {
244  Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
245  }
246  break;
247  }
248 
249  // Continuous projection
252  {
254 
255  for(i = 0; i < nVariables; ++i)
256  {
257  m_fields[i]->FwdTrans(inarray[i], coeffs);
258  m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
259  }
260  break;
261  }
262  default:
263  ASSERTL0(false, "Unknown projection scheme");
264  break;
265  }
266  }
267 
268  /**
269  * @brief Return the flux vector for the inviscid Burger equation.
270  *
271  * @param i Component of the flux vector to calculate.
272  * @param physfield Fields.
273  * @param flux Resulting flux.
274  */
276  const Array<OneD, Array<OneD, NekDouble> > &physfield,
278  {
279  const int nq = GetNpoints();
280 
281  for (int i = 0; i < flux.num_elements(); ++i)
282  {
283  for (int j = 0; j < flux[0].num_elements(); ++j)
284  {
285  Vmath::Vmul(nq, physfield[i], 1, physfield[i], 1,
286  flux[i][j], 1);
287  Vmath::Smul(nq, 0.5, flux[i][j], 1, flux[i][j], 1);
288  }
289  }
290  }
291 }
292 
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
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.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
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)
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.
UnsteadyInviscidBurger(const LibUtilities::SessionReaderSharedPtr &pSession)
Session reader.
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.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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
static std::string className
Name of class.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215