Nektar++
UnsteadyAdvection.cpp
Go to the documentation of this file.
1/////////////////////////////////////////////////////////////////////////////
2//
3// File: UnsteadyAdvection.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 linear advection solve routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
39#include <iostream>
40
41namespace Nektar
42{
43
46 "UnsteadyAdvection", UnsteadyAdvection::create,
47 "Unsteady Advection equation.");
48
52 : UnsteadySystem(pSession, pGraph), AdvectionSystem(pSession, pGraph)
53{
54}
55
56/**
57 * @brief Initialisation object for the unsteady linear advection equation.
58 */
59void UnsteadyAdvection::v_InitObject(bool DeclareFields)
60{
61 // Call to the initialisation object of UnsteadySystem
62 AdvectionSystem::v_InitObject(DeclareFields);
63
64 // Read the advection velocities from session file
65 m_session->LoadParameter("wavefreq", m_waveFreq, 0.0);
66
67 // check to see if it is explicity turned off
68 m_session->MatchSolverInfo("GJPStabilisation", "False",
70
71 // if GJPStabilisation set to False bool will be true and
72 // if not false so negate/revese bool
74
75 m_session->LoadParameter("GJPJumpScale", m_GJPJumpScale, 1.0);
76
77 // Define Velocity fields
78 std::vector<std::string> vel;
79 vel.push_back("Vx");
80 vel.push_back("Vy");
81 vel.push_back("Vz");
82
83 // Resize the advection velocities vector to dimension of the problem
84 vel.resize(m_spacedim);
85
86 // Store in the global variable m_velocity the advection velocities
88 GetFunction("AdvectionVelocity")->Evaluate(vel, m_velocity);
89
90 // Type of advection class to be used
91 switch (m_projectionType)
92 {
93 // Discontinuous field
95 {
96 // Do not forwards transform initial condition
97 m_homoInitialFwd = false;
98
99 // Define the normal velocity fields
100 if (m_fields[0]->GetTrace())
101 {
103 }
104
105 std::string advName;
106 std::string riemName;
107 m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
109 advName, advName);
111 {
112 m_advObject->SetFluxVector(
114 }
115 else
116 {
118 this);
119 }
120 m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
123 riemName, m_session);
124 m_riemannSolver->SetScalar(
126 m_advObject->SetRiemannSolver(m_riemannSolver);
127 m_advObject->InitObject(m_session, m_fields);
128 break;
129 }
130 // Continuous field
133 {
134 std::string advName;
135 m_session->LoadSolverInfo("AdvectionType", advName,
136 "NonConservative");
138 advName, advName);
140 {
141 m_advObject->SetFluxVector(
143 }
144 else
145 {
147 this);
148 }
149 break;
150 }
151 default:
152 {
153 ASSERTL0(false, "Unsupported projection type.");
154 break;
155 }
156 }
157
158 // Forcing terms
159 m_forcing = SolverUtils::Forcing::Load(m_session, shared_from_this(),
160 m_fields, m_fields.size());
161
162 // If explicit it computes RHS and PROJECTION for the time integration
164 {
167 }
168 // Otherwise it gives an error (no implicit integration)
169 else
170 {
171 ASSERTL0(false, "Implicit unsteady Advection not set up.");
172 }
173}
174
175/**
176 * @brief Compute the right-hand side for the linear advection equation.
177 *
178 * @param inarray Given fields.
179 * @param outarray Calculated solution.
180 * @param time Time.
181 */
183 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
184 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
185{
186 // Number of fields (variables of the problem)
187 int nVariables = inarray.size();
188
190 if (m_ALESolver)
191 {
192 timer.Start();
193 Array<OneD, Array<OneD, NekDouble>> tmpIn(nVariables);
194 // If ALE we must take Mu coefficient space to u physical space
196 auto advWeakDGObject =
197 std::dynamic_pointer_cast<SolverUtils::AdvectionWeakDG>(
199 advWeakDGObject->AdvectCoeffs(nVariables, m_fields, m_velocity, tmpIn,
200 outarray, time);
201 timer.Stop();
202 }
203 else
204 {
205 timer.Start();
206 m_advObject->Advect(nVariables, m_fields, m_velocity, inarray, outarray,
207 time);
208 timer.Stop();
209 }
210
211 // Elapsed time
212 timer.AccumulateRegion("Advect");
213
214 // Negate the RHS
215 for (int i = 0; i < nVariables; ++i)
216 {
217 Vmath::Neg(outarray[i].size(), outarray[i], 1);
218 }
219
220 // Add forcing terms
221 for (auto &x : m_forcing)
222 {
223 // set up non-linear terms
224 x->Apply(m_fields, inarray, outarray, time);
225 }
226}
227
228/**
229 * @brief Compute the projection for the linear advection equation.
230 *
231 * @param inarray Given fields.
232 * @param outarray Calculated solution.
233 * @param time Time.
234 */
236 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
237 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
238{
239 // Number of fields (variables of the problem)
240 int nVariables = inarray.size();
241
242 // Perform ALE movement
243 if (m_ALESolver)
244 {
246 }
247
248 // Set the boundary conditions
250
251 // Switch on the projection type (Discontinuous or Continuous)
252 switch (m_projectionType)
253 {
254 // Discontinuous projection
256 {
257 // Just copy over array
258 if (inarray != outarray)
259 {
260 int npoints = GetNpoints();
261
262 for (int i = 0; i < nVariables; ++i)
263 {
264 Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
265 }
266 }
267 break;
268 }
269 // Continuous projection
272 {
273 int ncoeffs = m_fields[0]->GetNcoeffs();
274 Array<OneD, NekDouble> coeffs(ncoeffs, 0.0);
276 {
279
280 Array<OneD, NekDouble> wsp(ncoeffs);
281
282 for (int i = 0; i < nVariables; ++i)
283 {
285 std::dynamic_pointer_cast<MultiRegions::ContField>(
286 m_fields[i]);
287
288 m_fields[i]->IProductWRTBase(inarray[i], wsp);
289
291 cfield->GetGJPForcing();
292
295
296 if (GJPData->IsSemiImplicit())
297 {
298 mtype = StdRegions::eMassGJP;
299 }
300
301 // to set up forcing need initial guess in
302 // physical space
304
305 GJPData->Apply(inarray[i], wsp, NullNekDouble1DArray,
306 scale);
307
308 // Solve the system
310 mtype, cfield->GetLocalToGlobalMap(), factors);
311
312 cfield->GlobalSolve(key, wsp, coeffs, NullNekDouble1DArray);
313
314 m_fields[i]->BwdTrans(coeffs, outarray[i]);
315 }
316 }
317 else
318 {
319 for (int i = 0; i < nVariables; ++i)
320 {
321 m_fields[i]->FwdTrans(inarray[i], coeffs);
322 m_fields[i]->BwdTrans(coeffs, outarray[i]);
323 }
324 }
325 break;
326 }
327 default:
328 ASSERTL0(false, "Unknown projection scheme");
329 break;
330 }
331}
332
333/**
334 * @brief Get the normal velocity for the linear advection equation.
335 */
337{
339 return m_traceVn;
340}
341
343 const Array<OneD, const Array<OneD, NekDouble>> &velfield)
344{
345 // Number of trace (interface) points
346 int nTracePts = GetTraceNpoints();
347 int nPts = m_velocity[0].size();
348
349 // Auxiliary variable to compute the normal velocity
350 Array<OneD, NekDouble> tmp(nPts), tmp2(nTracePts);
351
352 // Reset the normal velocity
353 Vmath::Zero(nTracePts, m_traceVn, 1);
354
355 for (int i = 0; i < velfield.size(); ++i)
356 {
357 // velocity - grid velocity for ALE before getting trace velocity
358 Vmath::Vsub(nPts, velfield[i], 1, m_gridVelocity[i], 1, tmp, 1);
359
360 m_fields[0]->ExtractTracePhys(tmp, tmp2);
361
362 Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, tmp2, 1, m_traceVn, 1,
363 m_traceVn, 1);
364 }
365
366 return m_traceVn;
367}
368
369/**
370 * @brief Return the flux vector for the linear advection equation.
371 *
372 * @param physfield Fields.
373 * @param flux Resulting flux.
374 */
376 const Array<OneD, Array<OneD, NekDouble>> &physfield,
378{
379 ASSERTL1(flux[0].size() == m_velocity.size(),
380 "Dimension of flux array and velocity array do not match");
381
382 const int nq = m_fields[0]->GetNpoints();
383
384 for (int i = 0; i < flux.size(); ++i)
385 {
386 for (int j = 0; j < flux[0].size(); ++j)
387 {
388 for (int k = 0; k < nq; ++k)
389 {
390 // If ALE we need to take off the grid velocity
391 flux[i][j][k] =
392 physfield[i][k] * (m_velocity[j][k] - m_gridVelocity[j][k]);
393 }
394 }
395 }
396}
397
398/**
399 * @brief Return the flux vector for the linear advection equation using
400 * the dealiasing technique.
401 *
402 * @param physfield Fields.
403 * @param flux Resulting flux.
404 */
406 const Array<OneD, Array<OneD, NekDouble>> &physfield,
408{
409 ASSERTL1(flux[0].size() == m_velocity.size(),
410 "Dimension of flux array and velocity array do not match");
411
412 int nq = physfield[0].size();
413 int nVariables = physfield.size();
414
415 // Factor to rescale 1d points in dealiasing
416 NekDouble OneDptscale = 2;
417
419
420 // Get number of points to dealias a cubic non-linearity
421 nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
422
423 // Initialisation of higher-space variables
424 Array<OneD, Array<OneD, NekDouble>> physfieldInterp(nVariables);
426 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> fluxInterp(nVariables);
427
428 // Interpolation to higher space of physfield
429 for (int i = 0; i < nVariables; ++i)
430 {
431 physfieldInterp[i] = Array<OneD, NekDouble>(nq);
433 for (int j = 0; j < m_expdim; ++j)
434 {
435 fluxInterp[i][j] = Array<OneD, NekDouble>(nq);
436 }
437
438 m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
439 physfieldInterp[i]);
440 }
441
442 // Interpolation to higher space of velocity
443 for (int j = 0; j < m_expdim; ++j)
444 {
445 velocityInterp[j] = Array<OneD, NekDouble>(nq);
446
447 m_fields[0]->PhysInterp1DScaled(OneDptscale, m_velocity[j],
448 velocityInterp[j]);
449 }
450
451 // Evaluation of flux vector in the higher space
452 for (int i = 0; i < flux.size(); ++i)
453 {
454 for (int j = 0; j < flux[0].size(); ++j)
455 {
456 Vmath::Vmul(nq, physfieldInterp[i], 1, velocityInterp[j], 1,
457 fluxInterp[i][j], 1);
458 }
459 }
460
461 // Galerkin project solution back to original space
462 for (int i = 0; i < nVariables; ++i)
463 {
464 for (int j = 0; j < m_spacedim; ++j)
465 {
466 m_fields[0]->PhysGalerkinProjection1DScaled(
467 OneDptscale, fluxInterp[i][j], flux[i][j]);
468 }
469 }
470}
471
473{
476 {
478 s, "GJP Stab. Impl. ",
479 m_session->GetSolverInfo("GJPStabilisation"));
480 SolverUtils::AddSummaryItem(s, "GJP Stab. JumpScale", m_GJPJumpScale);
481 }
482}
483
484bool UnsteadyAdvection::v_PreIntegrate([[maybe_unused]] int step)
485{
486 return false;
487}
488
490 std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
491 std::vector<std::string> &variables)
492{
493 bool extraFields;
494 m_session->MatchSolverInfo("OutputExtraFields", "True", extraFields, true);
495
496 if (extraFields && m_ALESolver)
497 {
498 ExtraFldOutputGridVelocity(fieldcoeffs, variables);
499 }
500}
501
504{
506 {
507 m_spaceDim = spaceDim;
508 m_fieldsALE = fields;
509
510 // Initialise grid velocities as 0s
513 for (int i = 0; i < spaceDim; ++i)
514 {
515 m_gridVelocity[i] =
516 Array<OneD, NekDouble>(fields[0]->GetTotPoints(), 0.0);
518 fields[0]->GetTrace()->GetTotPoints(), 0.0);
519 }
520 }
521 ALEHelper::InitObject(spaceDim, fields);
522}
523
524} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
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
Array< OneD, MultiRegions::ExpListSharedPtr > m_fieldsALE
Definition: ALEHelper.h:88
SOLVER_UTILS_EXPORT void InitObject(int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Definition: ALEHelper.cpp:48
Array< OneD, Array< OneD, NekDouble > > m_gridVelocityTrace
Definition: ALEHelper.h:90
SOLVER_UTILS_EXPORT void ALEDoElmtInvMassBwdTrans(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: ALEHelper.cpp:148
SOLVER_UTILS_EXPORT void ExtraFldOutputGridVelocity(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
Definition: ALEHelper.cpp:389
SOLVER_UTILS_EXPORT void MoveMesh(const NekDouble &time, Array< OneD, Array< OneD, NekDouble > > &traceNormals)
Definition: ALEHelper.cpp:168
Array< OneD, Array< OneD, NekDouble > > m_gridVelocity
Definition: ALEHelper.h:89
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).
int m_expdim
Expansion dimension.
NekDouble m_timestep
Time step size.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
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.
SOLVER_UTILS_EXPORT int GetTotPoints()
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.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
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.
void v_ALEInitObject(int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields) override
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables) override
Array< OneD, NekDouble > m_traceVn
Array< OneD, NekDouble > & GetNormalVel(const Array< OneD, const Array< OneD, NekDouble > > &velfield)
Get the normal velocity based on input velfield.
static std::string className
Name of class.
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
void GetFluxVectorDeAlias(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point using dealiasing.
UnsteadyAdvection(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
bool v_PreIntegrate(int step) override
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< GJPStabilisation > GJPStabilisationSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
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()
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
static Array< OneD, NekDouble > NullNekDouble1DArray
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 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
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220