Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
UnsteadyAdvectionDiffusion.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File UnsteadyAdvectionDiffusion.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 
40 using namespace std;
41 
42 namespace Nektar
43 {
44  string UnsteadyAdvectionDiffusion::className
46  "UnsteadyAdvectionDiffusion",
47  UnsteadyAdvectionDiffusion::create);
48 
49  UnsteadyAdvectionDiffusion::UnsteadyAdvectionDiffusion(
51  : UnsteadySystem(pSession),
52  AdvectionSystem(pSession)
53  {
54  m_planeNumber = 0;
55  }
56 
57  /**
58  * @brief Initialisation object for the unsteady linear advection
59  * diffusion equation.
60  */
62  {
64 
65  m_session->LoadParameter("wavefreq", m_waveFreq, 0.0);
66  m_session->LoadParameter("epsilon", m_epsilon, 0.0);
67 
68  // turn on substepping
69  m_session->MatchSolverInfo("Extrapolation", "SubStepping",
70  m_subSteppingScheme, false);
71 
72  // Define Velocity fields
74  std::vector<std::string> vel;
75  vel.push_back("Vx");
76  vel.push_back("Vy");
77  vel.push_back("Vz");
78  vel.resize(m_spacedim);
79 
80  EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
81 
82  m_session->MatchSolverInfo(
83  "SpectralVanishingViscosity", "True", m_useSpecVanVisc, false);
84 
86  {
87  m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
88  m_session->LoadParameter("SVVDiffCoeff",m_sVVDiffCoeff,0.1);
89  }
90 
91  // Type of advection and diffusion classes to be used
92  switch(m_projectionType)
93  {
94  // Discontinuous field
96  {
97  // Do not forwards transform initial condition
98  m_homoInitialFwd = false;
99 
100  // Advection term
101  string advName;
102  string riemName;
103  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
105  CreateInstance(advName, advName);
106  m_advObject->SetFluxVector(&UnsteadyAdvectionDiffusion::
107  GetFluxVectorAdv, this);
108  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
110  CreateInstance(riemName);
111  m_riemannSolver->SetScalar("Vn", &UnsteadyAdvectionDiffusion::
112  GetNormalVelocity, this);
113  m_advObject->SetRiemannSolver(m_riemannSolver);
114  m_advObject->InitObject (m_session, m_fields);
115 
116  // Diffusion term
118  {
119  std::string diffName;
120  m_session->LoadSolverInfo("DiffusionType", diffName, "LDG");
122  CreateInstance(diffName, diffName);
123  m_diffusion->SetFluxVector(&UnsteadyAdvectionDiffusion::
124  GetFluxVectorDiff, this);
125  m_diffusion->InitObject(m_session, m_fields);
126  }
127 
128  ASSERTL0(m_subSteppingScheme == false,"SubSteppingScheme is not set up for DG projection");
129  break;
130  }
131  // Continuous field
134  {
135  // Advection term
136  std::string advName;
137  m_session->LoadSolverInfo("AdvectionType", advName,
138  "NonConservative");
140  CreateInstance(advName, advName);
141  m_advObject->SetFluxVector(&UnsteadyAdvectionDiffusion::
142  GetFluxVectorAdv, this);
143 
144  if(advName.compare("WeakDG") == 0)
145  {
146  string riemName;
147  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
149  CreateInstance(riemName);
150  m_riemannSolver->SetScalar("Vn",
152  GetNormalVelocity, this);
153  m_advObject->SetRiemannSolver(m_riemannSolver);
154  m_advObject->InitObject (m_session, m_fields);
155  }
156 
157  // In case of Galerkin explicit diffusion gives an error
159  {
160  ASSERTL0(false, "Explicit Galerkin diffusion not set up.");
161  }
162  // In case of Galerkin implicit diffusion: do nothing
163  break;
164  }
165  default:
166  {
167  ASSERTL0(false, "Unsupported projection type.");
168  break;
169  }
170  }
171 
177 
178  if(m_subSteppingScheme) // Substepping
179  {
181  "Projection must be set to Mixed_CG_Discontinuous for "
182  "substepping");
184  m_intScheme->GetIntegrationMethod(), m_intScheme);
185  }
186  }
187 
188  /**
189  * @brief Unsteady linear advection diffusion equation destructor.
190  */
192  {
193  }
194 
195  /**
196  * @brief Get the normal velocity for the unsteady linear advection
197  * diffusion equation.
198  */
200  {
201  return GetNormalVel(m_velocity);
202  }
203 
204 
206  const Array<OneD, const Array<OneD, NekDouble> > &velfield)
207  {
208  // Number of trace (interface) points
209  int i;
210  int nTracePts = GetTraceNpoints();
211 
212  // Auxiliary variable to compute the normal velocity
213  Array<OneD, NekDouble> tmp(nTracePts);
214  m_traceVn = Array<OneD, NekDouble>(nTracePts, 0.0);
215 
216  // Reset the normal velocity
217  Vmath::Zero(nTracePts, m_traceVn, 1);
218 
219  for (i = 0; i < velfield.num_elements(); ++i)
220  {
221  m_fields[0]->ExtractTracePhys(velfield[i], tmp);
222 
223  Vmath::Vvtvp(nTracePts,
224  m_traceNormals[i], 1,
225  tmp, 1,
226  m_traceVn, 1,
227  m_traceVn, 1);
228  }
229 
230  return m_traceVn;
231  }
232 
233  /**
234  * @brief Compute the right-hand side for the unsteady linear advection
235  * diffusion problem.
236  *
237  * @param inarray Given fields.
238  * @param outarray Calculated solution.
239  * @param time Time.
240  */
242  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
243  Array<OneD, Array<OneD, NekDouble> >&outarray,
244  const NekDouble time)
245  {
246  // Number of fields (variables of the problem)
247  int nVariables = inarray.num_elements();
248 
249  // Number of solution points
250  int nSolutionPts = GetNpoints();
251 
252  // RHS computation using the new advection base class
253  m_advObject->Advect(nVariables, m_fields, m_velocity,
254  inarray, outarray, time);
255 
256  // Negate the RHS
257  for (int i = 0; i < nVariables; ++i)
258  {
259  Vmath::Neg(nSolutionPts, outarray[i], 1);
260  }
261 
263  {
264  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nVariables);
265  for (int i = 0; i < nVariables; ++i)
266  {
267  outarrayDiff[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
268  }
269 
270  m_diffusion->Diffuse(nVariables, m_fields, inarray, outarrayDiff);
271 
272  for (int i = 0; i < nVariables; ++i)
273  {
274  Vmath::Svtvp(nSolutionPts, m_epsilon, &outarrayDiff[i][0], 1,
275  &outarray[i][0], 1,
276  &outarray[i][0], 1);
277  }
278  }
279 
280  }
281 
282  /**
283  * @brief Compute the projection for the unsteady advection
284  * diffusion problem.
285  *
286  * @param inarray Given fields.
287  * @param outarray Calculated solution.
288  * @param time Time.
289  */
291  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
292  Array<OneD, Array<OneD, NekDouble> > &outarray,
293  const NekDouble time)
294  {
295  int i;
296  int nvariables = inarray.num_elements();
297  SetBoundaryConditions(time);
298  switch(m_projectionType)
299  {
301  {
302  // Just copy over array
303  int npoints = GetNpoints();
304 
305  for(i = 0; i < nvariables; ++i)
306  {
307  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
308  }
309  break;
310  }
313  {
315 
316  for(i = 0; i < nvariables; ++i)
317  {
318  m_fields[i]->FwdTrans(inarray[i], coeffs);
319  m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
320  }
321  break;
322  }
323  default:
324  {
325  ASSERTL0(false, "Unknown projection scheme");
326  break;
327  }
328  }
329  }
330 
331 
332  /* @brief Compute the diffusion term implicitly.
333  *
334  * @param inarray Given fields.
335  * @param outarray Calculated solution.
336  * @param time Time.
337  * @param lambda Diffusion coefficient.
338  */
340  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
341  Array<OneD, Array<OneD, NekDouble> >&outarray,
342  const NekDouble time,
343  const NekDouble lambda)
344  {
345  int nvariables = inarray.num_elements();
346  int nq = m_fields[0]->GetNpoints();
347 
349  factors[StdRegions::eFactorLambda] = 1.0/lambda/m_epsilon;
350 
351  if(m_useSpecVanVisc)
352  {
355  }
357  {
358  factors[StdRegions::eFactorTau] = 1.0;
359  }
360 
361  Array<OneD, Array< OneD, NekDouble> > F(nvariables);
362  for (int n = 0; n < nvariables; ++n)
363  {
364  F[n] = Array<OneD, NekDouble> (nq);
365  }
366 
367  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
368  // inarray = input: \hat{rhs} -> output: \hat{Y}
369  // outarray = output: nabla^2 \hat{Y}
370  // where \hat = modal coeffs
371  for (int i = 0; i < nvariables; ++i)
372  {
373  // Multiply 1.0/timestep/lambda
374  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda],
375  inarray[i], 1, F[i], 1);
376  }
377 
378  //Setting boundary conditions
379  SetBoundaryConditions(time);
380 
381  for (int i = 0; i < nvariables; ++i)
382  {
383  // Solve a system of equations with Helmholtz solver
384  m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(),
385  NullFlagList, factors);
386 
387  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
388  }
389  }
390 
391  /**
392  * @brief Return the flux vector for the advection part.
393  *
394  * @param physfield Fields.
395  * @param flux Resulting flux.
396  */
398  const Array<OneD, Array<OneD, NekDouble> > &physfield,
400  {
401  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
402  "Dimension of flux array and velocity array do not match");
403 
404  const int nq = m_fields[0]->GetNpoints();
405 
406  for (int i = 0; i < flux.num_elements(); ++i)
407  {
408  for (int j = 0; j < flux[0].num_elements(); ++j)
409  {
410  Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
411  flux[i][j], 1);
412  }
413  }
414  }
415 
416  /**
417  * @brief Return the flux vector for the diffusion part.
418  *
419  * @param i Equation number.
420  * @param j Spatial direction.
421  * @param physfield Fields.
422  * @param derivatives First order derivatives.
423  * @param flux Resulting flux.
424  */
426  const int i,
427  const int j,
428  const Array<OneD, Array<OneD, NekDouble> > &physfield,
429  Array<OneD, Array<OneD, NekDouble> > &derivatives,
431  {
432  for (int k = 0; k < flux.num_elements(); ++k)
433  {
434  Vmath::Zero(GetNpoints(), flux[k], 1);
435  }
436  Vmath::Vcopy(GetNpoints(), physfield[i], 1, flux[j], 1);
437  }
438 
441  {
443  }
444 
445 
446  /**
447  * Perform the extrapolation.
448  */
450  {
452  {
454  }
455 
456  return false;
457  }
458 
459 
460  /**
461  *
462  */
464  const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln,
465  int nstep,
466  NekDouble time)
467  {
468  int n;
469  int nsubsteps;
470 
471  NekDouble dt;
472 
473  Array<OneD, Array<OneD, NekDouble> > fields, velfields;
474 
475  static int ncalls = 1;
476  int nint = min(ncalls++, m_intSteps);
477 
480 
481  LibUtilities::CommSharedPtr comm = m_session->GetComm();
482 
483  // Get the proper time step with CFL control
484  dt = GetSubstepTimeStep();
485 
486  nsubsteps = (m_timestep > dt)? ((int)(m_timestep/dt)+1):1;
487  nsubsteps = max(m_minsubsteps, nsubsteps);
488 
489  dt = m_timestep/nsubsteps;
490 
491  if (m_infosteps && !((nstep+1)%m_infosteps) && comm->GetRank() == 0)
492  {
493  cout << "Sub-integrating using "<< nsubsteps
494  << " steps over Dt = " << m_timestep
495  << " (SubStep CFL=" << m_cflSafetyFactor << ")"<< endl;
496  }
497 
498  for (int m = 0; m < nint; ++m)
499  {
500  // We need to update the fields held by the m_integrationSoln
501  fields = integrationSoln->UpdateSolutionVector()[m];
502 
503  // Initialise NS solver which is set up to use a GLM method
504  // with calls to EvaluateAdvection_SetPressureBCs and
505  // SolveUnsteadyStokesSystem
507  SubIntegrationSoln = m_subStepIntegrationScheme->
508  InitializeScheme(dt, fields, time, m_subStepIntegrationOps);
509 
510  for(n = 0; n < nsubsteps; ++n)
511  {
512  fields = m_subStepIntegrationScheme->TimeIntegrate(n, dt, SubIntegrationSoln,
514  }
515 
516  // Reset time integrated solution in m_integrationSoln
517  integrationSoln->SetSolVector(m,fields);
518  }
519  }
520 
521 
522  /**
523  *
524  */
526  {
527  int n_element = m_fields[0]->GetExpSize();
528 
529  const Array<OneD, int> ExpOrder=m_fields[0]->EvalBasisNumModesMaxPerExp();
530  Array<OneD, int> ExpOrderList (n_element, ExpOrder);
531 
532  const NekDouble cLambda = 0.2; // Spencer book pag. 317
533 
534  Array<OneD, NekDouble> tstep (n_element, 0.0);
535  Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
536 
537  stdVelocity = GetMaxStdVelocity(m_velocity);
538 
539  for(int el = 0; el < n_element; ++el)
540  {
541  tstep[el] = m_cflSafetyFactor /
542  (stdVelocity[el] * cLambda *
543  (ExpOrder[el]-1) * (ExpOrder[el]-1));
544  }
545 
546  NekDouble TimeStep = Vmath::Vmin(n_element, tstep, 1);
547  m_session->GetComm()->AllReduce(TimeStep,LibUtilities::ReduceMin);
548 
549  return TimeStep;
550  }
551 
553  int intMethod,
554  const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
555  {
556  // Set to 1 for first step and it will then be increased in
557  // time advance routines
558  switch(intMethod)
559  {
562  {
564 
565  }
566  break;
568  {
570  }
571  break;
572  default:
573  ASSERTL0(0,"Integration method not suitable: Options include BackwardEuler or BDFImplicitOrder1");
574  break;
575  }
576  m_intSteps = IntegrationScheme->GetIntegrationSteps();
577 
578  // set explicit time-integration class operators
581  }
582 
583  /**
584  * Explicit Advection terms used by SubStepAdvance time integration
585  */
587  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
588  Array<OneD, Array<OneD, NekDouble> > &outarray,
589  const NekDouble time)
590  {
591  int i;
592  int nVariables = inarray.num_elements();
593 
594  /// Get the number of coefficients
595  int ncoeffs = m_fields[0]->GetNcoeffs();
596 
597  /// Define an auxiliary variable to compute the RHS
598  Array<OneD, Array<OneD, NekDouble> > WeakAdv(nVariables);
599  WeakAdv[0] = Array<OneD, NekDouble> (ncoeffs*nVariables);
600  for(i = 1; i < nVariables; ++i)
601  {
602  WeakAdv[i] = WeakAdv[i-1] + ncoeffs;
603  }
604 
605  // Currently assume velocity field is time independent and does not therefore
606  // need extrapolating.
607  // RHS computation using the advection base class
608  m_advObject->Advect(nVariables, m_fields, m_velocity,
609  inarray, outarray, time);
610 
611  for(i = 0; i < nVariables; ++i)
612  {
613  m_fields[i]->IProductWRTBase(outarray[i],WeakAdv[i]);
614  // negation requried due to sign of DoAdvection term to be consistent
615  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
616  }
617 
618  AddAdvectionPenaltyFlux(m_velocity, inarray, WeakAdv);
619 
620 
621  /// Operations to compute the RHS
622  for(i = 0; i < nVariables; ++i)
623  {
624  // Negate the RHS
625  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
626 
627  /// Multiply the flux by the inverse of the mass matrix
628  m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
629 
630  /// Store in outarray the physical values of the RHS
631  m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
632  }
633  }
634 
635  /**
636  * Projection used by SubStepAdvance time integration
637  */
639  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
640  Array<OneD, Array<OneD, NekDouble> > &outarray,
641  const NekDouble time)
642  {
643  ASSERTL1(inarray.num_elements() == outarray.num_elements(),"Inarray and outarray of different sizes ");
644 
645  for(int i = 0; i < inarray.num_elements(); ++i)
646  {
647  Vmath::Vcopy(inarray[i].num_elements(),inarray[i],1,outarray[i],1);
648  }
649  }
650 
652  const Array<OneD, const Array<OneD, NekDouble> > &velfield,
653  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
654  Array<OneD, Array<OneD, NekDouble> > &Outarray)
655  {
656  ASSERTL1(physfield.num_elements() == Outarray.num_elements(),
657  "Physfield and outarray are of different dimensions");
658 
659  int i;
660 
661  /// Number of trace points
662  int nTracePts = m_fields[0]->GetTrace()->GetNpoints();
663 
664  /// Forward state array
665  Array<OneD, NekDouble> Fwd(3*nTracePts);
666 
667  /// Backward state array
668  Array<OneD, NekDouble> Bwd = Fwd + nTracePts;
669 
670  /// upwind numerical flux state array
671  Array<OneD, NekDouble> numflux = Bwd + nTracePts;
672 
673  /// Normal velocity array
674  Array<OneD, NekDouble> Vn = GetNormalVel(velfield);
675 
676  for(i = 0; i < physfield.num_elements(); ++i)
677  {
678  /// Extract forwards/backwards trace spaces
679  /// Note: Needs to have correct i value to get boundary conditions
680  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
681 
682  /// Upwind between elements
683  m_fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux);
684 
685  /// Construct difference between numflux and Fwd,Bwd
686  Vmath::Vsub(nTracePts, numflux, 1, Fwd, 1, Fwd, 1);
687  Vmath::Vsub(nTracePts, numflux, 1, Bwd, 1, Bwd, 1);
688 
689  /// Calculate the numerical fluxes multipling Fwd, Bwd and
690  /// numflux by the normal advection velocity
691  Vmath::Vmul(nTracePts, Fwd, 1, Vn, 1, Fwd, 1);
692  Vmath::Vmul(nTracePts, Bwd, 1, Vn, 1, Bwd, 1);
693 
694  m_fields[0]->AddFwdBwdTraceIntegral(Fwd,Bwd,Outarray[i]);
695  }
696  }
697 
698 
700  const Array<OneD, Array<OneD,NekDouble> > inarray)
701  {
702 
703  int n_points_0 = m_fields[0]->GetExp(0)->GetTotPoints();
704  int n_element = m_fields[0]->GetExpSize();
705  int nvel = inarray.num_elements();
706  int cnt;
707 
708  ASSERTL0(nvel >= 2, "Method not implemented for 1D");
709 
710  NekDouble pntVelocity;
711 
712  // Getting the standard velocity vector on the 2D normal space
713  Array<OneD, Array<OneD, NekDouble> > stdVelocity(nvel);
714  Array<OneD, NekDouble> maxV(n_element, 0.0);
716 
717  for (int i = 0; i < nvel; ++i)
718  {
719  stdVelocity[i] = Array<OneD, NekDouble>(n_points_0);
720  }
721 
722  if (nvel == 2)
723  {
724  cnt = 0.0;
725  for (int el = 0; el < n_element; ++el)
726  {
727  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
728  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
729 
730  // reset local space if necessary
731  if(n_points != n_points_0)
732  {
733  for (int j = 0; j < nvel; ++j)
734  {
735  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
736  }
737  n_points_0 = n_points;
738  }
739 
741  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
742 
743  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
745  {
746  for (int i = 0; i < n_points; i++)
747  {
748  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
749  + gmat[2][i]*inarray[1][i+cnt];
750 
751  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
752  + gmat[3][i]*inarray[1][i+cnt];
753  }
754  }
755  else
756  {
757  for (int i = 0; i < n_points; i++)
758  {
759  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
760  + gmat[2][0]*inarray[1][i+cnt];
761 
762  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
763  + gmat[3][0]*inarray[1][i+cnt];
764  }
765  }
766 
767  cnt += n_points;
768 
769 
770  for (int i = 0; i < n_points; i++)
771  {
772  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
773  + stdVelocity[1][i]*stdVelocity[1][i];
774 
775  if (pntVelocity>maxV[el])
776  {
777  maxV[el] = pntVelocity;
778  }
779  }
780  maxV[el] = sqrt(maxV[el]);
781  }
782  }
783  else
784  {
785  cnt = 0;
786  for (int el = 0; el < n_element; ++el)
787  {
788 
789  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
790  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
791 
792  // reset local space if necessary
793  if(n_points != n_points_0)
794  {
795  for (int j = 0; j < nvel; ++j)
796  {
797  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
798  }
799  n_points_0 = n_points;
800  }
801 
803  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
804 
805  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
807  {
808  for (int i = 0; i < n_points; i++)
809  {
810  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
811  + gmat[3][i]*inarray[1][i+cnt]
812  + gmat[6][i]*inarray[2][i+cnt];
813 
814  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
815  + gmat[4][i]*inarray[1][i+cnt]
816  + gmat[7][i]*inarray[2][i+cnt];
817 
818  stdVelocity[2][i] = gmat[2][i]*inarray[0][i+cnt]
819  + gmat[5][i]*inarray[1][i+cnt]
820  + gmat[8][i]*inarray[2][i+cnt];
821  }
822  }
823  else
824  {
825  for (int i = 0; i < n_points; i++)
826  {
827  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
828  + gmat[3][0]*inarray[1][i+cnt]
829  + gmat[6][0]*inarray[2][i+cnt];
830 
831  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
832  + gmat[4][0]*inarray[1][i+cnt]
833  + gmat[7][0]*inarray[2][i+cnt];
834 
835  stdVelocity[2][i] = gmat[2][0]*inarray[0][i+cnt]
836  + gmat[5][0]*inarray[1][i+cnt]
837  + gmat[8][0]*inarray[2][i+cnt];
838  }
839  }
840 
841  cnt += n_points;
842 
843  for (int i = 0; i < n_points; i++)
844  {
845  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
846  + stdVelocity[1][i]*stdVelocity[1][i]
847  + stdVelocity[2][i]*stdVelocity[2][i];
848 
849  if (pntVelocity > maxV[el])
850  {
851  maxV[el] = pntVelocity;
852  }
853  }
854 
855  maxV[el] = sqrt(maxV[el]);
856  //cout << maxV[el]*maxV[el] << endl;
857  }
858  }
859 
860  return maxV;
861  }
862 }
void SubStepProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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.
Array< OneD, NekDouble > GetMaxStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
BDF multi-step scheme of order 1 (implicit)
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
NekDouble m_time
Current time of simulation.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
LibUtilities::TimeIntegrationWrapperSharedPtr m_subStepIntegrationScheme
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_timestep
Time step size.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:871
Array< OneD, Array< OneD, NekDouble > > m_velocity
boost::shared_ptr< TimeIntegrationWrapper > TimeIntegrationWrapperSharedPtr
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
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.
STL namespace.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
void SubStepAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Array< OneD, NekDouble > & GetNormalVel(const Array< OneD, const Array< OneD, NekDouble > > &velfield)
Get the normal velocity based on input velfield.
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity based on m_velocity.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Perform the projection.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
virtual void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
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
virtual bool v_PreIntegrate(int step)
PreIntegration step for substepping.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
RiemannSolverFactory & GetRiemannSolverFactory()
int m_spacedim
Spatial dimension (>= expansion dim).
void SetUpSubSteppingTimeIntegration(int intMethod, const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
SolverUtils::DiffusionSharedPtr m_diffusion
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
TimeIntegrationWrapperFactory & GetTimeIntegrationWrapperFactory()
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
EquationSystemFactory & GetEquationSystemFactory()
void SubStepAdvance(const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln, int nstep, NekDouble time)
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
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.cpp:343
BDF multi-step scheme of order 2 (implicit)
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.
SOLVER_UTILS_EXPORT int GetNpoints()
LibUtilities::TimeIntegrationSchemeOperators m_subStepIntegrationOps
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetExpSize()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
void AddAdvectionPenaltyFlux(const Array< OneD, const Array< OneD, NekDouble > > &velfield, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &Outarray)
SOLVER_UTILS_EXPORT int GetNcoeffs()
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
Wrapper to the time integration scheme.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
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.
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.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Geometry is curved or has non-constant factors.
virtual void v_InitObject()
Initialise the object.
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
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.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215