Nektar++
UnsteadyViscousBurgers.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: UnsteadyViscousBurgers.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 viscous Burgers solve routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar
38{
39
42 "UnsteadyViscousBurgers", UnsteadyViscousBurgers::create,
43 "Viscous Burgers equation");
44
48 : UnsteadySystem(pSession, pGraph),
49 UnsteadyInviscidBurgers(pSession, pGraph)
50{
51}
52
53/**
54 * @brief Initialisation object for the viscous Burgers equation.
55 */
57{
59
60 m_session->LoadParameter("epsilon", m_epsilon, 1.0);
61 m_session->MatchSolverInfo("SpectralVanishingViscosity", "True",
62 m_useSpecVanVisc, false);
63 m_session->MatchSolverInfo("SpectralVanishingViscosity", "VarDiff",
65
67 {
68 m_useSpecVanVisc = true;
69 }
70
72 {
73 m_session->LoadParameter("SVVCutoffRatio", m_sVVCutoffRatio, 0.75);
74 m_session->LoadParameter("SVVDiffCoeff", m_sVVDiffCoeff, 0.1);
75 }
76
77 // Type of diffusion classes to be used
78 switch (m_projectionType)
79 {
80 // Discontinuous field
82 {
84 {
85 std::string diffName;
86 m_session->LoadSolverInfo("DiffusionType", diffName, "LDG");
88 diffName, diffName);
89 m_diffusion->SetFluxVector(
91 m_diffusion->InitObject(m_session, m_fields);
92 }
93 break;
94 }
95 // Continuous field
97 {
98 std::string advName;
99 m_session->LoadSolverInfo("AdvectionType", advName,
100 "NonConservative");
101 if (advName.compare("WeakDG") == 0)
102 {
103 // Define the normal velocity fields
104 if (m_fields[0]->GetTrace())
105 {
107 }
108
109 std::string riemName;
110 m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
113 riemName, m_session);
114 m_riemannSolver->SetScalar(
116 m_advObject->SetRiemannSolver(m_riemannSolver);
117 m_advObject->InitObject(m_session, m_fields);
118 }
119
120 // In case of Galerkin explicit diffusion gives an error
122 {
123 ASSERTL0(false, "Explicit Galerkin diffusion not set up.");
124 }
125 break;
126 }
127 default:
128 {
129 ASSERTL0(false, "Unsupported projection type.");
130 break;
131 }
132 }
133
136}
137
138/**
139 * @brief Compute the right-hand side for the viscous Burgers equation.
140 *
141 * @param inarray Given fields.
142 * @param outarray Calculated solution.
143 * @param time Time.
144 */
146 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
147 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
148{
149 // Number of fields (variables of the problem)
150 int nVariables = inarray.size();
151
152 // Number of solution points
153 int nSolutionPts = GetNpoints();
154
155 UnsteadyInviscidBurgers::DoOdeRhs(inarray, outarray, time);
156
158 {
159 Array<OneD, Array<OneD, NekDouble>> outarrayDiff(nVariables);
160
161 for (int i = 0; i < nVariables; ++i)
162 {
163 outarrayDiff[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
164 }
165
166 m_diffusion->Diffuse(nVariables, m_fields, inarray, outarrayDiff);
167
168 for (int i = 0; i < nVariables; ++i)
169 {
170 Vmath::Vadd(nSolutionPts, &outarrayDiff[i][0], 1, &outarray[i][0],
171 1, &outarray[i][0], 1);
172 }
173 }
174}
175
176/* @brief Compute the diffusion term implicitly.
177 *
178 * @param inarray Given fields.
179 * @param outarray Calculated solution.
180 * @param time Time.
181 * @param lambda Diffusion coefficient.
182 */
184 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
185 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
186 const NekDouble lambda)
187{
188 int nvariables = inarray.size();
189 int nq = m_fields[0]->GetNpoints();
190
193
195 {
198 }
199
201 {
203 }
204
206 for (int n = 0; n < nvariables; ++n)
207 {
208 F[n] = Array<OneD, NekDouble>(nq);
209 }
210
211 // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
212 // inarray = input: \hat{rhs} -> output: \hat{Y}
213 // outarray = output: nabla^2 \hat{Y}
214 // where \hat = modal coeffs
215 for (int i = 0; i < nvariables; ++i)
216 {
217 // Multiply 1.0/timestep/lambda
218 Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1,
219 F[i], 1);
220 }
221
222 // Setting boundary conditions
224
226 {
227 static int cnt = 0;
228
229 if (cnt % 10 == 0)
230 {
232 for (int i = 0; i < m_fields.size(); ++i)
233 {
234 m_fields[i]->ClearGlobalLinSysManager();
235 vel[i] = m_fields[i]->UpdatePhys();
236 }
238 }
239 ++cnt;
240 }
241
242 for (int i = 0; i < nvariables; ++i)
243 {
244 // Solve a system of equations with Helmholtz solver
245 m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(), factors,
247
248 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
249 }
250}
251
252/**
253 * @brief Return the flux vector for the diffusion part.
254 *
255 */
257 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray,
258 const Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &qfield,
259 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &viscousTensor)
260{
261 unsigned int nDim = qfield.size();
262 unsigned int nConvectiveFields = qfield[0].size();
263 unsigned int nPts = qfield[0][0].size();
264 for (unsigned int j = 0; j < nDim; ++j)
265 {
266 for (unsigned int i = 0; i < nConvectiveFields; ++i)
267 {
268 Vmath::Smul(nPts, m_epsilon, qfield[j][i], 1, viscousTensor[j][i],
269 1);
270 }
271 }
272}
273
275{
278 {
279 std::stringstream ss;
280 ss << "SVV (cut off = " << m_sVVCutoffRatio
281 << ", coeff = " << m_sVVDiffCoeff << ")";
282 SolverUtils::AddSummaryItem(s, "Smoothing", ss.str());
283 }
284}
285
286} // 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 DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
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.
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff(const Array< OneD, Array< OneD, NekDouble > > vel, StdRegions::VarCoeffMap &varCoeffMap)
Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h t...
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
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 SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
void GetFluxVectorDiff(const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
Evaluate the flux at each solution point for the diffusion part.
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, NekDouble time, NekDouble lambda)
Solve implicitly the diffusion term.
static std::string className
Name of class.
SolverUtils::DiffusionSharedPtr m_diffusion
UnsteadyViscousBurgers(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
StdRegions::VarCoeffMap m_varCoeffLap
Variable Coefficient map for the Laplacian which can be activated as part of SVV or otherwise.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:39
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
StdRegions::ConstFactorMap factors
double NekDouble
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