Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 namespace Nektar
42 {
45  "UnsteadyViscousBurgers",
47 
50  : UnsteadySystem(pSession),
51  AdvectionSystem(pSession),
52  m_varCoeffLap(StdRegions::NullVarCoeffMap)
53  {
54  m_planeNumber = 0;
55  }
56 
57  /**
58  * @brief Initialisation object for the unsteady linear advection
59  * diffusion equation.
60  */
62  {
63  AdvectionSystem::v_InitObject();
64 
65  m_session->LoadParameter("wavefreq", m_waveFreq, 0.0);
66  m_session->LoadParameter("epsilon", m_epsilon, 0.0);
67 
68  m_session->MatchSolverInfo(
69  "SpectralVanishingViscosity", "True", m_useSpecVanVisc, false);
70  m_session->MatchSolverInfo(
71  "SpectralVanishingViscosity", "VarDiff", m_useSpecVanViscVarDiff, false);
73  {
74  m_useSpecVanVisc = true;
75  }
76 
78  {
79  m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
80  m_session->LoadParameter("SVVDiffCoeff",m_sVVDiffCoeff,0.1);
81  }
82 
83  // Type of advection and diffusion classes to be used
84  switch(m_projectionType)
85  {
86  // Discontinuous field
88  {
89  ASSERTL0(false,"Need to implement for DG");
90  // Do not forwards transform initial condition
91  m_homoInitialFwd = false;
92 
93  // Advection term
94  string advName;
95  string riemName;
96  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
98  CreateInstance(advName, advName);
99  m_advObject->SetFluxVector(&UnsteadyViscousBurgers::
100  GetFluxVectorAdv, this);
101  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
103  CreateInstance(riemName);
104  m_advObject->SetRiemannSolver(m_riemannSolver);
105  m_advObject->InitObject (m_session, m_fields);
106 
107  // Diffusion term
108  std::string diffName;
109  m_session->LoadSolverInfo("DiffusionType", diffName, "LDG");
111  CreateInstance(diffName, diffName);
112  m_diffusion->SetFluxVector(&UnsteadyViscousBurgers::
113  GetFluxVectorDiff, this);
114  m_diffusion->InitObject(m_session, m_fields);
115  break;
116  }
117  // Continuous field
120  {
121  // Advection term
122  std::string advName;
123  m_session->LoadSolverInfo("AdvectionType", advName,
124  "NonConservative");
126  CreateInstance(advName, advName);
127  m_advObject->SetFluxVector(&UnsteadyViscousBurgers::
128  GetFluxVectorAdv, this);
129 
131  {
132  Array<OneD, Array<OneD, NekDouble> > vel(m_fields.num_elements());
133  for(int i = 0; i < m_fields.num_elements(); ++i)
134  {
135  vel[i] = m_fields[i]->UpdatePhys();
136  }
138  }
139 
140  // In case of Galerkin explicit diffusion gives an error
142  {
143  ASSERTL0(false, "Explicit Galerkin diffusion not set up.");
144  }
145  // In case of Galerkin implicit diffusion: do nothing
146  break;
147  }
148  default:
149  {
150  ASSERTL0(false, "Unsupported projection type.");
151  break;
152  }
153  }
154 
155  // Forcing terms
157  m_fields.num_elements());
158 
161 
163  m_explicitDiffusion == 1)
164  {
166  }
167  }
168 
169  /**
170  * @brief Unsteady linear advection diffusion equation destructor.
171  */
173  {
174  }
175 
176  /**
177  * @brief Get the normal velocity for the unsteady linear advection
178  * diffusion equation.
179  */
181  Array<OneD, Array<OneD, NekDouble> >&inarray)
182  {
183  // Number of trace (interface) points
184  int i;
185  int nTracePts = GetTraceNpoints();
186 
187  // Auxiliary variable to compute the normal velocity
188  Array<OneD, NekDouble> tmp(nTracePts);
189  m_traceVn = Array<OneD, NekDouble>(nTracePts, 0.0);
190 
191  // Reset the normal velocity
192  Vmath::Zero(nTracePts, m_traceVn, 1);
193 
194  for (i = 0; i < inarray.num_elements(); ++i)
195  {
196  m_fields[0]->ExtractTracePhys(inarray[i], tmp);
197 
198  Vmath::Vvtvp(nTracePts,
199  m_traceNormals[i], 1,
200  tmp, 1,
201  m_traceVn, 1,
202  m_traceVn, 1);
203  }
204 
205  return m_traceVn;
206  }
207 
208  /**
209  * @brief Compute the right-hand side for the unsteady linear advection
210  * diffusion problem.
211  *
212  * @param inarray Given fields.
213  * @param outarray Calculated solution.
214  * @param time Time.
215  */
217  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
218  Array<OneD, Array<OneD, NekDouble> >&outarray,
219  const NekDouble time)
220  {
221  // Number of fields (variables of the problem)
222  int nVariables = inarray.num_elements();
223 
224  // Number of solution points
225  int nSolutionPts = GetNpoints();
226 
227  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nVariables);
228 
229  for (int i = 0; i < nVariables; ++i)
230  {
231  outarrayDiff[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
232  }
233 
234  // RHS computation using the new advection base class
235  m_advObject->Advect(nVariables, m_fields, inarray,
236  inarray, outarray, time);
237 
238  // Negate the RHS
239  for (int i = 0; i < nVariables; ++i)
240  {
241  Vmath::Neg(nSolutionPts, outarray[i], 1);
242  }
243 
244  // No explicit diffusion for CG
246  {
247  m_diffusion->Diffuse(nVariables, m_fields, inarray, outarrayDiff);
248 
249  for (int i = 0; i < nVariables; ++i)
250  {
251  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
252  &outarrayDiff[i][0], 1, &outarray[i][0], 1);
253  }
254  }
255 
256  // Add forcing terms
257  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
258  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
259  {
260  // set up non-linear terms
261  (*x)->Apply(m_fields, inarray, outarray, time);
262  }
263  }
264 
265  /**
266  * @brief Compute the projection for the unsteady advection
267  * diffusion problem.
268  *
269  * @param inarray Given fields.
270  * @param outarray Calculated solution.
271  * @param time Time.
272  */
274  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
275  Array<OneD, Array<OneD, NekDouble> > &outarray,
276  const NekDouble time)
277  {
278  int i;
279  int nvariables = inarray.num_elements();
280  SetBoundaryConditions(time);
281  switch(m_projectionType)
282  {
284  {
285  // Just copy over array
286  int npoints = GetNpoints();
287 
288  for(i = 0; i < nvariables; ++i)
289  {
290  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
291  }
292  break;
293  }
296  {
297  // Do nothing for the moment.
298  }
299  default:
300  {
301  ASSERTL0(false, "Unknown projection scheme");
302  break;
303  }
304  }
305  }
306 
307 
308  /* @brief Compute the diffusion term implicitly.
309  *
310  * @param inarray Given fields.
311  * @param outarray Calculated solution.
312  * @param time Time.
313  * @param lambda Diffusion coefficient.
314  */
316  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
317  Array<OneD, Array<OneD, NekDouble> >&outarray,
318  const NekDouble time,
319  const NekDouble lambda)
320  {
321  int nvariables = inarray.num_elements();
322  int nq = m_fields[0]->GetNpoints();
323 
325  factors[StdRegions::eFactorLambda] = 1.0/lambda/m_epsilon;
326 
327  if(m_useSpecVanVisc)
328  {
331  }
332 
333  Array<OneD, Array< OneD, NekDouble> > F(nvariables);
334  F[0] = Array<OneD, NekDouble> (nq*nvariables);
335 
336  for (int n = 1; n < nvariables; ++n)
337  {
338  F[n] = F[n-1] + nq;
339  }
340 
341  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
342  // inarray = input: \hat{rhs} -> output: \hat{Y}
343  // outarray = output: nabla^2 \hat{Y}
344  // where \hat = modal coeffs
345  for (int i = 0; i < nvariables; ++i)
346  {
347  // Multiply 1.0/timestep/lambda
348  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda],
349  inarray[i], 1, F[i], 1);
350  }
351 
352  //Setting boundary conditions
353  SetBoundaryConditions(time);
354 
356  {
357  static int cnt = 0;
358 
359  if(cnt %10 == 0)
360  {
361  Array<OneD, Array<OneD, NekDouble> > vel(m_fields.num_elements());
362  for(int i = 0; i < m_fields.num_elements(); ++i)
363  {
364  m_fields[i]->ClearGlobalLinSysManager();
365  vel[i] = m_fields[i]->UpdatePhys();
366  }
368  }
369  ++cnt;
370  }
371  for (int i = 0; i < nvariables; ++i)
372  {
373  // Solve a system of equations with Helmholtz solver
374  m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(),
375  NullFlagList, factors, m_varCoeffLap);
376 
377  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
378  }
379  }
380 
381  /**
382  * @brief Return the flux vector for the advection part.
383  *
384  * @param physfield Fields.
385  * @param flux Resulting flux.
386  */
388  const Array<OneD, Array<OneD, NekDouble> > &physfield,
390  {
391 
392  const int nq = m_fields[0]->GetNpoints();
393 
394  for (int i = 0; i < flux.num_elements(); ++i)
395  {
396  for (int j = 0; j < flux[0].num_elements(); ++j)
397  {
398  Vmath::Vmul(nq, physfield[i], 1, physfield[j], 1,
399  flux[i][j], 1);
400  }
401  }
402  }
403 
404  /**
405  * @brief Return the flux vector for the diffusion part.
406  *
407  * @param i Equation number.
408  * @param j Spatial direction.
409  * @param physfield Fields.
410  * @param derivatives First order derivatives.
411  * @param flux Resulting flux.
412  */
414  const int i,
415  const int j,
416  const Array<OneD, Array<OneD, NekDouble> > &physfield,
417  Array<OneD, Array<OneD, NekDouble> > &derivatives,
419  {
420  for (int k = 0; k < flux.num_elements(); ++k)
421  {
422  Vmath::Zero(GetNpoints(), flux[k], 1);
423  }
424  Vmath::Vcopy(GetNpoints(), physfield[i], 1, flux[j], 1);
425  }
426 
429  {
430  AdvectionSystem::v_GenerateSummary(s);
431  if(m_useSpecVanVisc)
432  {
433  stringstream ss;
434  ss << "SVV (cut off = " << m_sVVCutoffRatio
435  << ", coeff = " << m_sVVDiffCoeff << ")";
436  AddSummaryItem(s, "Smoothing", ss.str());
437  }
438  }
439 }
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:161
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:428
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Array< OneD, NekDouble > m_traceVn
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:251
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:84
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
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)
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:199
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
UnsteadyViscousBurgers(const LibUtilities::SessionReaderSharedPtr &pSession)
Session reader.
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:382
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.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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:1047
static std::string className
Name of class.
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:285
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:169
static FlagList NullFlagList
An empty flag list.
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:227
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215