Nektar++
UnsteadyDiffusion.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: UnsteadyDiffusion.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 diffusion solve routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar
38{
39
42 "UnsteadyDiffusion", UnsteadyDiffusion::create);
43
47 : UnsteadySystem(pSession, pGraph)
48{
49}
50
51/**
52 * @brief Initialisation object for the unsteady diffusion problem.
53 */
54void UnsteadyDiffusion::v_InitObject(bool DeclareFields)
55{
56 UnsteadySystem::v_InitObject(DeclareFields);
57
58 // Load diffusion parameter
59 m_session->LoadParameter("epsilon", m_epsilon, 1.0);
60
61 m_session->MatchSolverInfo("SpectralVanishingViscosity", "True",
62 m_useSpecVanVisc, false);
63
65 {
66 m_session->LoadParameter("SVVCutoffRatio", m_sVVCutoffRatio, 0.75);
67 m_session->LoadParameter("SVVDiffCoeff", m_sVVDiffCoeff, 0.1);
68 }
69
70 int npoints = m_fields[0]->GetNpoints();
71
72 if (m_session->DefinesParameter("d00"))
73 {
74 m_d00 = m_session->GetParameter("d00");
76 Array<OneD, NekDouble>(npoints, m_session->GetParameter("d00"));
77 }
78 if (m_session->DefinesParameter("d11"))
79 {
80 m_d11 = m_session->GetParameter("d11");
82 Array<OneD, NekDouble>(npoints, m_session->GetParameter("d11"));
83 }
84 if (m_session->DefinesParameter("d22"))
85 {
86 m_d22 = m_session->GetParameter("d22");
88 Array<OneD, NekDouble>(npoints, m_session->GetParameter("d22"));
89 }
90
91 switch (m_projectionType)
92 {
94 {
95 std::string diffName;
96
97 // Do not forwards transform initial condition
98 m_homoInitialFwd = false;
99
100 m_session->LoadSolverInfo("DiffusionType", diffName, "LDG");
102 diffName, diffName);
103 m_diffusion->SetFluxVector(&UnsteadyDiffusion::GetFluxVector, this);
104 m_diffusion->InitObject(m_session, m_fields);
105 break;
106 }
108 {
109 // In case of Galerkin explicit diffusion gives an error
111 {
112 ASSERTL0(false, "Explicit Galerkin diffusion not set up.");
113 }
114 // In case of Galerkin implicit diffusion: do nothing
115 break;
116 }
117 default:
118 {
119 ASSERTL0(false, "Unknown projection scheme");
120 break;
121 }
122 }
123
127}
128
129/* @brief Compute the right-hand side for the unsteady diffusion problem.
130 *
131 * @param inarray Given fields.
132 * @param outarray Calculated solution.
133 * @param time Time.
134 */
136 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
138 [[maybe_unused]] const NekDouble time)
139{
140 // Number of fields (variables of the problem)
141 int nVariables = inarray.size();
142
143 // RHS computation using the new diffusion base class
144 m_diffusion->Diffuse(nVariables, m_fields, inarray, outarray);
145}
146
147/**
148 * @brief Compute the projection for the unsteady diffusion problem.
149 *
150 * @param inarray Given fields.
151 * @param outarray Calculated solution.
152 * @param time Time.
153 */
155 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
156 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
157{
158 int i;
159 int nvariables = inarray.size();
161
162 switch (m_projectionType)
163 {
165 {
166 // Just copy over array
167 if (inarray != outarray)
168 {
169 int npoints = GetNpoints();
170
171 for (i = 0; i < nvariables; ++i)
172 {
173 Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
174 }
175 }
176 break;
177 }
179 {
181
182 for (i = 0; i < nvariables; ++i)
183 {
184 m_fields[i]->FwdTrans(inarray[i], coeffs);
185 m_fields[i]->BwdTrans(coeffs, outarray[i]);
186 }
187 break;
188 }
189 default:
190 {
191 ASSERTL0(false, "Unknown projection scheme");
192 break;
193 }
194 }
195}
196
197/**
198 * @brief Implicit solution of the unsteady diffusion problem.
199 */
201 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
203 [[maybe_unused]] const NekDouble time, const NekDouble lambda)
204{
206
207 int nvariables = inarray.size();
208 int npoints = m_fields[0]->GetNpoints();
211
213 {
216 }
217
218 // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
219 // inarray = input: \hat{rhs} -> output: \hat{Y}
220 // outarray = output: nabla^2 \hat{Y}
221 // where \hat = modal coeffs
222 for (int i = 0; i < nvariables; ++i)
223 {
224 // Multiply 1.0/timestep/lambda
225 Vmath::Smul(npoints, -factors[StdRegions::eFactorLambda], inarray[i], 1,
226 outarray[i], 1);
227
228 // Solve a system of equations with Helmholtz solver
229 m_fields[i]->HelmSolve(outarray[i], m_fields[i]->UpdateCoeffs(),
231
232 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
233
234 m_fields[i]->SetPhysState(false);
235 }
236}
237
238/**
239 * @brief Return the flux vector for the unsteady diffusion problem.
240 */
242 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray,
243 const Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &qfield,
244 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &viscousTensor)
245{
246 unsigned int nDim = qfield.size();
247 unsigned int nConvectiveFields = qfield[0].size();
248 unsigned int nPts = qfield[0][0].size();
249
250 NekDouble d[3] = {m_d00, m_d11, m_d22};
251
252 for (unsigned int j = 0; j < nDim; ++j)
253 {
254 for (unsigned int i = 0; i < nConvectiveFields; ++i)
255 {
256 Vmath::Smul(nPts, m_epsilon * d[j], qfield[j][i], 1,
257 viscousTensor[j][i], 1);
258 }
259 }
260}
261
263{
266 {
267 std::stringstream ss;
268 ss << "SVV (cut off = " << m_sVVCutoffRatio
269 << ", coeff = " << m_sVVDiffCoeff << ")";
270 AddSummaryItem(s, "Smoothing", ss.str());
271 }
272}
273
274} // 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 DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
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.
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 v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, NekDouble time, NekDouble lambda)
Implicit solution of the unsteady diffusion problem.
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
Return the flux vector for the unsteady diffusion problem.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
StdRegions::VarCoeffMap m_varcoeff
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void v_InitObject(bool DeclareFields=true) override
Initialisation object for the unsteady diffusion problem.
UnsteadyDiffusion(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
static std::string className
Name of class.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection for the unsteady diffusion problem.
SolverUtils::DiffusionSharedPtr m_diffusion
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
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
std::vector< double > d(NPUPPER *NPUPPER)
StdRegions::ConstFactorMap factors
double NekDouble
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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825