Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // 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 advection-diffusion solve routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
39 #include <StdRegions/StdQuadExp.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45  string UnsteadyViscousBurgers::className
47  "UnsteadyViscousBurgers",
48  UnsteadyViscousBurgers::create);
49 
50  UnsteadyViscousBurgers::UnsteadyViscousBurgers(
52  : UnsteadySystem(pSession),
53  AdvectionSystem(pSession),
54  m_varCoeffLap(StdRegions::NullVarCoeffMap)
55  {
56  m_planeNumber = 0;
57  }
58 
59  /**
60  * @brief Initialisation object for the unsteady linear advection
61  * diffusion equation.
62  */
64  {
66 
67  m_session->LoadParameter("wavefreq", m_waveFreq, 0.0);
68  m_session->LoadParameter("epsilon", m_epsilon, 0.0);
69 
70  m_session->MatchSolverInfo(
71  "SpectralVanishingViscosity", "True", m_useSpecVanVisc, false);
72  m_session->MatchSolverInfo(
73  "SpectralVanishingViscosity", "VarDiff", m_useSpecVanViscVarDiff, false);
75  {
76  m_useSpecVanVisc = true;
77  }
78 
80  {
81  m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
82  m_session->LoadParameter("SVVDiffCoeff",m_sVVDiffCoeff,0.1);
83  }
84 
85  // Type of advection and diffusion classes to be used
86  switch(m_projectionType)
87  {
88  // Discontinuous field
90  {
91  ASSERTL0(false,"Need to implement for DG");
92  // Do not forwards transform initial condition
93  m_homoInitialFwd = false;
94 
95  // Advection term
96  string advName;
97  string riemName;
98  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
100  CreateInstance(advName, advName);
101  m_advObject->SetFluxVector(&UnsteadyViscousBurgers::
102  GetFluxVectorAdv, this);
103  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
105  CreateInstance(riemName);
106  m_advObject->SetRiemannSolver(m_riemannSolver);
107  m_advObject->InitObject (m_session, m_fields);
108 
109  // Diffusion term
110  std::string diffName;
111  m_session->LoadSolverInfo("DiffusionType", diffName, "LDG");
113  CreateInstance(diffName, diffName);
114  m_diffusion->SetFluxVector(&UnsteadyViscousBurgers::
115  GetFluxVectorDiff, this);
116  m_diffusion->InitObject(m_session, m_fields);
117  break;
118  }
119  // Continuous field
122  {
123  // Advection term
124  std::string advName;
125  m_session->LoadSolverInfo("AdvectionType", advName,
126  "NonConservative");
128  CreateInstance(advName, advName);
129  m_advObject->SetFluxVector(&UnsteadyViscousBurgers::
130  GetFluxVectorAdv, this);
131 
133  {
134  Array<OneD, Array<OneD, NekDouble> > vel(m_fields.num_elements());
135  for(int i = 0; i < m_fields.num_elements(); ++i)
136  {
137  vel[i] = m_fields[i]->UpdatePhys();
138  }
140  }
141 
142  // In case of Galerkin explicit diffusion gives an error
144  {
145  ASSERTL0(false, "Explicit Galerkin diffusion not set up.");
146  }
147  // In case of Galerkin implicit diffusion: do nothing
148  break;
149  }
150  default:
151  {
152  ASSERTL0(false, "Unsupported projection type.");
153  break;
154  }
155  }
156 
157  // Forcing terms
159  m_fields.num_elements());
160 
163 
165  m_explicitDiffusion == 1)
166  {
168  }
169  }
170 
171  /**
172  * @brief Unsteady linear advection diffusion equation destructor.
173  */
175  {
176  }
177 
178  /**
179  * @brief Get the normal velocity for the unsteady linear advection
180  * diffusion equation.
181  */
183  Array<OneD, Array<OneD, NekDouble> >&inarray)
184  {
185  // Number of trace (interface) points
186  int i;
187  int nTracePts = GetTraceNpoints();
188 
189  // Auxiliary variable to compute the normal velocity
190  Array<OneD, NekDouble> tmp(nTracePts);
191  m_traceVn = Array<OneD, NekDouble>(nTracePts, 0.0);
192 
193  // Reset the normal velocity
194  Vmath::Zero(nTracePts, m_traceVn, 1);
195 
196  for (i = 0; i < inarray.num_elements(); ++i)
197  {
198  m_fields[0]->ExtractTracePhys(inarray[i], tmp);
199 
200  Vmath::Vvtvp(nTracePts,
201  m_traceNormals[i], 1,
202  tmp, 1,
203  m_traceVn, 1,
204  m_traceVn, 1);
205  }
206 
207  return m_traceVn;
208  }
209 
210  /**
211  * @brief Compute the right-hand side for the unsteady linear advection
212  * diffusion problem.
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  // Number of fields (variables of the problem)
224  int nVariables = inarray.num_elements();
225 
226  // Number of solution points
227  int nSolutionPts = GetNpoints();
228 
229  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nVariables);
230 
231  for (int i = 0; i < nVariables; ++i)
232  {
233  outarrayDiff[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
234  }
235 
236  // RHS computation using the new advection base class
237  m_advObject->Advect(nVariables, m_fields, inarray,
238  inarray, outarray, time);
239 
240  // Negate the RHS
241  for (int i = 0; i < nVariables; ++i)
242  {
243  Vmath::Neg(nSolutionPts, outarray[i], 1);
244  }
245 
246  // No explicit diffusion for CG
248  {
249  m_diffusion->Diffuse(nVariables, m_fields, inarray, outarrayDiff);
250 
251  for (int i = 0; i < nVariables; ++i)
252  {
253  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
254  &outarrayDiff[i][0], 1, &outarray[i][0], 1);
255  }
256  }
257 
258  // Add forcing terms
259  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
260  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
261  {
262  // set up non-linear terms
263  (*x)->Apply(m_fields, inarray, outarray, time);
264  }
265  }
266 
267  /**
268  * @brief Compute the projection for the unsteady advection
269  * diffusion problem.
270  *
271  * @param inarray Given fields.
272  * @param outarray Calculated solution.
273  * @param time Time.
274  */
276  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
277  Array<OneD, Array<OneD, NekDouble> > &outarray,
278  const NekDouble time)
279  {
280  int i;
281  int nvariables = inarray.num_elements();
282  SetBoundaryConditions(time);
283  switch(m_projectionType)
284  {
286  {
287  // Just copy over array
288  int npoints = GetNpoints();
289 
290  for(i = 0; i < nvariables; ++i)
291  {
292  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
293  }
294  break;
295  }
298  {
299  // Do nothing for the moment.
300  }
301  default:
302  {
303  ASSERTL0(false, "Unknown projection scheme");
304  break;
305  }
306  }
307  }
308 
309 
310  /* @brief Compute the diffusion term implicitly.
311  *
312  * @param inarray Given fields.
313  * @param outarray Calculated solution.
314  * @param time Time.
315  * @param lambda Diffusion coefficient.
316  */
318  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
319  Array<OneD, Array<OneD, NekDouble> >&outarray,
320  const NekDouble time,
321  const NekDouble lambda)
322  {
323  int nvariables = inarray.num_elements();
324  int nq = m_fields[0]->GetNpoints();
325 
327  factors[StdRegions::eFactorLambda] = 1.0/lambda/m_epsilon;
328 
329  if(m_useSpecVanVisc)
330  {
333  }
334 
335  Array<OneD, Array< OneD, NekDouble> > F(nvariables);
336  F[0] = Array<OneD, NekDouble> (nq*nvariables);
337 
338  for (int n = 1; n < nvariables; ++n)
339  {
340  F[n] = F[n-1] + nq;
341  }
342 
343  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
344  // inarray = input: \hat{rhs} -> output: \hat{Y}
345  // outarray = output: nabla^2 \hat{Y}
346  // where \hat = modal coeffs
347  for (int i = 0; i < nvariables; ++i)
348  {
349  // Multiply 1.0/timestep/lambda
350  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda],
351  inarray[i], 1, F[i], 1);
352  }
353 
354  //Setting boundary conditions
355  SetBoundaryConditions(time);
356 
358  {
359  static int cnt = 0;
360 
361  if(cnt %10 == 0)
362  {
363  Array<OneD, Array<OneD, NekDouble> > vel(m_fields.num_elements());
364  for(int i = 0; i < m_fields.num_elements(); ++i)
365  {
366  m_fields[i]->ClearGlobalLinSysManager();
367  vel[i] = m_fields[i]->UpdatePhys();
368  }
370  }
371  ++cnt;
372  }
373  for (int i = 0; i < nvariables; ++i)
374  {
375  // Solve a system of equations with Helmholtz solver
376  m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(),
377  NullFlagList, factors, m_varCoeffLap);
378 
379  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
380  }
381  }
382 
383  /**
384  * @brief Return the flux vector for the advection part.
385  *
386  * @param physfield Fields.
387  * @param flux Resulting flux.
388  */
390  const Array<OneD, Array<OneD, NekDouble> > &physfield,
392  {
393 
394  const int nq = m_fields[0]->GetNpoints();
395 
396  for (int i = 0; i < flux.num_elements(); ++i)
397  {
398  for (int j = 0; j < flux[0].num_elements(); ++j)
399  {
400  Vmath::Vmul(nq, physfield[i], 1, physfield[j], 1,
401  flux[i][j], 1);
402  }
403  }
404  }
405 
406  /**
407  * @brief Return the flux vector for the diffusion part.
408  *
409  * @param i Equation number.
410  * @param j Spatial direction.
411  * @param physfield Fields.
412  * @param derivatives First order derivatives.
413  * @param flux Resulting flux.
414  */
416  const int i,
417  const int j,
418  const Array<OneD, Array<OneD, NekDouble> > &physfield,
419  Array<OneD, Array<OneD, NekDouble> > &derivatives,
421  {
422  for (int k = 0; k < flux.num_elements(); ++k)
423  {
424  Vmath::Zero(GetNpoints(), flux[k], 1);
425  }
426  Vmath::Vcopy(GetNpoints(), physfield[i], 1, flux[j], 1);
427  }
428 
431  {
433  if(m_useSpecVanVisc)
434  {
435  stringstream ss;
436  ss << "SVV (cut off = " << m_sVVCutoffRatio
437  << ", coeff = " << m_sVVDiffCoeff << ")";
438  AddSummaryItem(s, "Smoothing", ss.str());
439  }
440  }
441 }
void GetFluxVectorDiff(const int i, const int j, const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &derivatives, Array< OneD, Array< OneD, NekDouble > > &flux)
Evaluate the flux at each solution point for the diffusion part.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
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:442
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Array< OneD, NekDouble > m_traceVn
STL namespace.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
StdRegions::VarCoeffMap m_varCoeffLap
Variable Coefficient map for the Laplacian which can be activated as part of SVV or otherwise...
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:86
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
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...
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
virtual ~UnsteadyViscousBurgers()
Destructor.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
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:213
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
RiemannSolverFactory & GetRiemannSolverFactory()
void GetFluxVectorAdv(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point for the advection part.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Perform the projection.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
EquationSystemFactory & GetEquationSystemFactory()
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
virtual void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
Array< OneD, NekDouble > & GetNormalVelocity(Array< OneD, Array< OneD, NekDouble > > &inarray)
Get the normal velocity.
SolverUtils::DiffusionSharedPtr m_diffusion
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
virtual void v_InitObject()
Initialise the object.
virtual 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.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
A base class for PDEs which include an advection component.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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.cpp:299
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
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:183
static FlagList NullFlagList
An empty flag list.
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:228
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215