Nektar++
UnsteadyInviscidBurgers.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: UnsteadyInviscidBurgers.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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Unsteady inviscid Burgers solve routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38namespace Nektar
39{
40
43 "UnsteadyInviscidBurgers", UnsteadyInviscidBurgers::create,
44 "Inviscid Burgers equation");
45
49 : UnsteadySystem(pSession, pGraph), AdvectionSystem(pSession, pGraph)
50{
51}
52
53/**
54 * @brief Initialisation object for the inviscid Burgers equation.
55 */
57{
58 AdvectionSystem::v_InitObject(DeclareFields);
59
60 // Type of advection class to be used
61 switch (m_projectionType)
62 {
63 // Discontinuous field
65 {
66 // Do not forwards transform initial condition
67 m_homoInitialFwd = false;
68
69 // Define the normal velocity fields
70 if (m_fields[0]->GetTrace())
71 {
73 }
74
75 // Advection term
76 std::string advName;
77 std::string riemName;
78 m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
80 advName, advName);
82 this);
83 m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
86 riemName, m_session);
87 m_riemannSolver->SetScalar(
89 m_advObject->SetRiemannSolver(m_riemannSolver);
90 m_advObject->InitObject(m_session, m_fields);
91 break;
92 }
93 // Continuous field
95 {
96 std::string advName;
97 m_session->LoadSolverInfo("AdvectionType", advName,
98 "NonConservative");
100 advName, advName);
102 this);
103 break;
104 }
105 default:
106 {
107 ASSERTL0(false, "Unsupported projection type.");
108 break;
109 }
110 }
111
112 // Forcing terms
113 m_forcing = SolverUtils::Forcing::Load(m_session, shared_from_this(),
114 m_fields, m_fields.size());
115
117 {
120 }
121 else
122 {
123 ASSERTL0(false, "Implicit unsteady Advection not set up.");
124 }
125}
126
127/**
128 * @brief Compute the right-hand side for the inviscid Burgers equation.
129 *
130 * @param inarray Given fields.
131 * @param outarray Calculated solution.
132 * @param time Time.
133 */
135 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
136 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
137{
138 // Number of fields (variables of the problem)
139 int nVariables = inarray.size();
140
141 // Number of solution points
142 int nSolutionPts = GetNpoints();
143
145 timer.Start();
146 // RHS computation using the new advection base class
147 m_advObject->Advect(nVariables, m_fields, inarray, inarray, outarray, time);
148 timer.Stop();
149 // Elapsed time
150 timer.AccumulateRegion("Advect");
151
152 // Negate the RHS
153 for (int i = 0; i < nVariables; ++i)
154 {
155 Vmath::Neg(nSolutionPts, outarray[i], 1);
156 }
157
158 // Add forcing terms
159 for (auto &x : m_forcing)
160 {
161 // set up non-linear terms
162 x->Apply(m_fields, inarray, outarray, time);
163 }
164}
165
166/**
167 * @brief Compute the projection for the inviscid Burgers equation.
168 *
169 * @param inarray Given fields.
170 * @param outarray Calculated solution.
171 * @param time Time.
172 */
174 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
175 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
176{
177 // Number of variables of the problem
178 int nVariables = inarray.size();
179
180 // Set the boundary conditions
182
183 // Switch on the projection type (Discontinuous or Continuous)
184 switch (m_projectionType)
185 {
186 // Discontinuous projection
188 {
189 // Just copy over array
190 if (inarray != outarray)
191 {
192 int npoints = GetNpoints();
193
194 for (int i = 0; i < nVariables; ++i)
195 {
196 Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
197 }
198 }
199 break;
200 }
201 // Continuous projection
203 {
205
206 for (int i = 0; i < nVariables; ++i)
207 {
208 m_fields[i]->FwdTrans(inarray[i], coeffs);
209 m_fields[i]->BwdTrans(coeffs, outarray[i]);
210 }
211 break;
212 }
213 default:
214 {
215 ASSERTL0(false, "Unknown projection scheme");
216 break;
217 }
218 }
219}
220
221/**
222 * @brief Get the normal velocity for the inviscid Burgers equation.
223 */
225{
226 // Number of trace (interface) points
227 int nTracePts = GetTraceNpoints();
228
229 // Number of solution points
230 int nSolutionPts = GetNpoints();
231
232 // Auxiliary variables to compute the normal velocity
233 Array<OneD, NekDouble> Fwd(nTracePts);
234 Array<OneD, NekDouble> Bwd(nTracePts);
235 Array<OneD, NekDouble> physfield(nSolutionPts);
236
237 // Reset the normal velocity
238 Vmath::Zero(nTracePts, m_traceVn, 1);
239
240 for (int i = 0; i < m_spacedim; ++i)
241 {
242 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), physfield);
243 m_fields[i]->GetFwdBwdTracePhys(physfield, Fwd, Bwd, true);
244 Vmath::Vadd(nTracePts, Fwd, 1, Bwd, 1, Fwd, 1);
245 Vmath::Smul(nTracePts, 0.5, Fwd, 1, Fwd, 1);
246 Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Fwd, 1, m_traceVn, 1,
247 m_traceVn, 1);
248 }
249 Vmath::Smul(nTracePts, 0.5, m_traceVn, 1, m_traceVn, 1);
250
251 return m_traceVn;
252}
253
254/**
255 * @brief Return the flux vector for the inviscid Burgers equation.
256 *
257 * @param physfield Fields.
258 * @param flux Resulting flux.
259 */
261 const Array<OneD, Array<OneD, NekDouble>> &physfield,
263{
264 const int nq = GetNpoints();
265
266 for (int i = 0; i < flux.size(); ++i)
267 {
268 for (int j = 0; j < flux[0].size(); ++j)
269 {
270 Vmath::Vmul(nq, physfield[i], 1, physfield[i], 1, flux[i][j], 1);
271 Vmath::Smul(nq, 0.5, flux[i][j], 1, flux[i][j], 1);
272 }
273 }
274}
275
277{
279}
280
281} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:70
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:76
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
UnsteadyInviscidBurgers(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
static std::string className
Name of class.
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
EquationSystemFactory & GetEquationSystemFactory()
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
double NekDouble
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.hpp:72
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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.hpp:366
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825