Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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
117  std::string diffName;
118  m_session->LoadSolverInfo("DiffusionType", diffName, "LDG");
120  CreateInstance(diffName, diffName);
121  m_diffusion->SetFluxVector(&UnsteadyAdvectionDiffusion::
122  GetFluxVectorDiff, this);
123  m_diffusion->InitObject(m_session, m_fields);
124 
125  ASSERTL0(m_subSteppingScheme == false,"SubSteppingScheme is not set up for DG projection");
126  break;
127  }
128  // Continuous field
131  {
132  // Advection term
133  std::string advName;
134  m_session->LoadSolverInfo("AdvectionType", advName,
135  "NonConservative");
137  CreateInstance(advName, advName);
138  m_advObject->SetFluxVector(&UnsteadyAdvectionDiffusion::
139  GetFluxVectorAdv, this);
140 
141  if(advName.compare("WeakDG") == 0)
142  {
143  string riemName;
144  m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
146  CreateInstance(riemName);
147  m_riemannSolver->SetScalar("Vn",
149  GetNormalVelocity, this);
150  m_advObject->SetRiemannSolver(m_riemannSolver);
151  m_advObject->InitObject (m_session, m_fields);
152  }
153 
154  // In case of Galerkin explicit diffusion gives an error
156  {
157  ASSERTL0(false, "Explicit Galerkin diffusion not set up.");
158  }
159  // In case of Galerkin implicit diffusion: do nothing
160  break;
161  }
162  default:
163  {
164  ASSERTL0(false, "Unsupported projection type.");
165  break;
166  }
167  }
168 
171 
172  if(m_subSteppingScheme) // Substepping
173  {
175  "Projection must be set to Mixed_CG_Discontinuous for "
176  "substepping");
178  m_intScheme->GetIntegrationMethod(), m_intScheme);
179 
180  }
181  else // Standard velocity correction scheme
182  {
184  }
185 
187  m_explicitDiffusion == 1)
188  {
190  }
191  }
192 
193  /**
194  * @brief Unsteady linear advection diffusion equation destructor.
195  */
197  {
198  }
199 
200  /**
201  * @brief Get the normal velocity for the unsteady linear advection
202  * diffusion equation.
203  */
205  {
206  return GetNormalVel(m_velocity);
207  }
208 
209 
211  const Array<OneD, const Array<OneD, NekDouble> > &velfield)
212  {
213  // Number of trace (interface) points
214  int i;
215  int nTracePts = GetTraceNpoints();
216 
217  // Auxiliary variable to compute the normal velocity
218  Array<OneD, NekDouble> tmp(nTracePts);
219  m_traceVn = Array<OneD, NekDouble>(nTracePts, 0.0);
220 
221  // Reset the normal velocity
222  Vmath::Zero(nTracePts, m_traceVn, 1);
223 
224  for (i = 0; i < velfield.num_elements(); ++i)
225  {
226  m_fields[0]->ExtractTracePhys(velfield[i], tmp);
227 
228  Vmath::Vvtvp(nTracePts,
229  m_traceNormals[i], 1,
230  tmp, 1,
231  m_traceVn, 1,
232  m_traceVn, 1);
233  }
234 
235  return m_traceVn;
236  }
237 
238  /**
239  * @brief Compute the right-hand side for the unsteady linear advection
240  * diffusion problem.
241  *
242  * @param inarray Given fields.
243  * @param outarray Calculated solution.
244  * @param time Time.
245  */
247  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
248  Array<OneD, Array<OneD, NekDouble> >&outarray,
249  const NekDouble time)
250  {
251  // Number of fields (variables of the problem)
252  int nVariables = inarray.num_elements();
253 
254  // Number of solution points
255  int nSolutionPts = GetNpoints();
256 
257  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nVariables);
258 
259  for (int i = 0; i < nVariables; ++i)
260  {
261  outarrayDiff[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
262  }
263 
264  // RHS computation using the new advection base class
265  m_advObject->Advect(nVariables, m_fields, m_velocity,
266  inarray, outarray, time);
267 
268  // Negate the RHS
269  for (int i = 0; i < nVariables; ++i)
270  {
271  Vmath::Neg(nSolutionPts, outarray[i], 1);
272  }
273 
274  // No explicit diffusion for CG
276  {
277  m_diffusion->Diffuse(nVariables, m_fields, inarray, outarrayDiff);
278 
279  for (int i = 0; i < nVariables; ++i)
280  {
281  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
282  &outarrayDiff[i][0], 1, &outarray[i][0], 1);
283  }
284  }
285 
286  }
287 
288  /**
289  * @brief Compute the projection for the unsteady advection
290  * diffusion problem.
291  *
292  * @param inarray Given fields.
293  * @param outarray Calculated solution.
294  * @param time Time.
295  */
297  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
298  Array<OneD, Array<OneD, NekDouble> > &outarray,
299  const NekDouble time)
300  {
301  int i;
302  int nvariables = inarray.num_elements();
303  SetBoundaryConditions(time);
304  switch(m_projectionType)
305  {
307  {
308  // Just copy over array
309  int npoints = GetNpoints();
310 
311  for(i = 0; i < nvariables; ++i)
312  {
313  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
314  }
315  break;
316  }
319  {
321 
322  for(i = 0; i < nvariables; ++i)
323  {
324  m_fields[i]->FwdTrans(inarray[i], coeffs);
325  m_fields[i]->BwdTrans_IterPerExp(coeffs, outarray[i]);
326  }
327  break;
328  }
329  default:
330  {
331  ASSERTL0(false, "Unknown projection scheme");
332  break;
333  }
334  }
335  }
336 
337 
338  /* @brief Compute the diffusion term implicitly.
339  *
340  * @param inarray Given fields.
341  * @param outarray Calculated solution.
342  * @param time Time.
343  * @param lambda Diffusion coefficient.
344  */
346  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
347  Array<OneD, Array<OneD, NekDouble> >&outarray,
348  const NekDouble time,
349  const NekDouble lambda)
350  {
351  int nvariables = inarray.num_elements();
352  int nq = m_fields[0]->GetNpoints();
353 
355  factors[StdRegions::eFactorLambda] = 1.0/lambda/m_epsilon;
356 
357  if(m_useSpecVanVisc)
358  {
361  }
362 
363  Array<OneD, Array< OneD, NekDouble> > F(nvariables);
364  F[0] = Array<OneD, NekDouble> (nq*nvariables);
365 
366  for (int n = 1; n < nvariables; ++n)
367  {
368  F[n] = F[n-1] + nq;
369  }
370 
371  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
372  // inarray = input: \hat{rhs} -> output: \hat{Y}
373  // outarray = output: nabla^2 \hat{Y}
374  // where \hat = modal coeffs
375  for (int i = 0; i < nvariables; ++i)
376  {
377  // Multiply 1.0/timestep/lambda
378  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda],
379  inarray[i], 1, F[i], 1);
380  }
381 
382  //Setting boundary conditions
383  SetBoundaryConditions(time);
384 
385  for (int i = 0; i < nvariables; ++i)
386  {
387  // Solve a system of equations with Helmholtz solver
388  m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(),
389  NullFlagList, factors);
390 
391  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
392  }
393  }
394 
395  /**
396  * @brief Return the flux vector for the advection part.
397  *
398  * @param physfield Fields.
399  * @param flux Resulting flux.
400  */
402  const Array<OneD, Array<OneD, NekDouble> > &physfield,
404  {
405  ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
406  "Dimension of flux array and velocity array do not match");
407 
408  const int nq = m_fields[0]->GetNpoints();
409 
410  for (int i = 0; i < flux.num_elements(); ++i)
411  {
412  for (int j = 0; j < flux[0].num_elements(); ++j)
413  {
414  Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
415  flux[i][j], 1);
416  }
417  }
418  }
419 
420  /**
421  * @brief Return the flux vector for the diffusion part.
422  *
423  * @param i Equation number.
424  * @param j Spatial direction.
425  * @param physfield Fields.
426  * @param derivatives First order derivatives.
427  * @param flux Resulting flux.
428  */
430  const int i,
431  const int j,
432  const Array<OneD, Array<OneD, NekDouble> > &physfield,
433  Array<OneD, Array<OneD, NekDouble> > &derivatives,
435  {
436  for (int k = 0; k < flux.num_elements(); ++k)
437  {
438  Vmath::Zero(GetNpoints(), flux[k], 1);
439  }
440  Vmath::Vcopy(GetNpoints(), physfield[i], 1, flux[j], 1);
441  }
442 
445  {
447  }
448 
449 
450  /**
451  * Perform the extrapolation.
452  */
454  {
456  {
458  }
459 
460  return false;
461  }
462 
463 
464  /**
465  *
466  */
468  const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln,
469  int nstep,
470  NekDouble time)
471  {
472  int n;
473  int nsubsteps;
474 
475  NekDouble dt;
476 
477  Array<OneD, Array<OneD, NekDouble> > fields, velfields;
478 
479  static int ncalls = 1;
480  int nint = min(ncalls++, m_intSteps);
481 
484 
485  LibUtilities::CommSharedPtr comm = m_session->GetComm();
486 
487  // Get the proper time step with CFL control
488  dt = GetSubstepTimeStep();
489 
490  nsubsteps = (m_timestep > dt)? ((int)(m_timestep/dt)+1):1;
491  nsubsteps = max(m_minsubsteps, nsubsteps);
492 
493  dt = m_timestep/nsubsteps;
494 
495  if (m_infosteps && !((nstep+1)%m_infosteps) && comm->GetRank() == 0)
496  {
497  cout << "Sub-integrating using "<< nsubsteps
498  << " steps over Dt = " << m_timestep
499  << " (SubStep CFL=" << m_cflSafetyFactor << ")"<< endl;
500  }
501 
502  for (int m = 0; m < nint; ++m)
503  {
504  // We need to update the fields held by the m_integrationSoln
505  fields = integrationSoln->UpdateSolutionVector()[m];
506 
507  // Initialise NS solver which is set up to use a GLM method
508  // with calls to EvaluateAdvection_SetPressureBCs and
509  // SolveUnsteadyStokesSystem
511  SubIntegrationSoln = m_subStepIntegrationScheme->
512  InitializeScheme(dt, fields, time, m_subStepIntegrationOps);
513 
514  for(n = 0; n < nsubsteps; ++n)
515  {
516  fields = m_subStepIntegrationScheme->TimeIntegrate(n, dt, SubIntegrationSoln,
518  }
519 
520  // Reset time integrated solution in m_integrationSoln
521  integrationSoln->SetSolVector(m,fields);
522  }
523  }
524 
525 
526  /**
527  *
528  */
530  {
531  int n_element = m_fields[0]->GetExpSize();
532 
533  const Array<OneD, int> ExpOrder=m_fields[0]->EvalBasisNumModesMaxPerExp();
534  Array<OneD, int> ExpOrderList (n_element, ExpOrder);
535 
536  const NekDouble cLambda = 0.2; // Spencer book pag. 317
537 
538  Array<OneD, NekDouble> tstep (n_element, 0.0);
539  Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
540 
541  stdVelocity = GetMaxStdVelocity(m_velocity);
542 
543  for(int el = 0; el < n_element; ++el)
544  {
545  tstep[el] = m_cflSafetyFactor /
546  (stdVelocity[el] * cLambda *
547  (ExpOrder[el]-1) * (ExpOrder[el]-1));
548  }
549 
550  NekDouble TimeStep = Vmath::Vmin(n_element, tstep, 1);
551  m_session->GetComm()->AllReduce(TimeStep,LibUtilities::ReduceMin);
552 
553  return TimeStep;
554  }
555 
557  int intMethod,
558  const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
559  {
560  // Set to 1 for first step and it will then be increased in
561  // time advance routines
562  switch(intMethod)
563  {
566  {
568 
569  }
570  break;
572  {
574  }
575  break;
576  default:
577  ASSERTL0(0,"Integration method not suitable: Options include BackwardEuler or BDFImplicitOrder1");
578  break;
579  }
580  m_intSteps = IntegrationScheme->GetIntegrationSteps();
581 
582  // set explicit time-integration class operators
585  }
586 
587  /**
588  * Explicit Advection terms used by SubStepAdvance time integration
589  */
591  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
592  Array<OneD, Array<OneD, NekDouble> > &outarray,
593  const NekDouble time)
594  {
595  int i;
596  int nVariables = inarray.num_elements();
597 
598  /// Get the number of coefficients
599  int ncoeffs = m_fields[0]->GetNcoeffs();
600 
601  /// Define an auxiliary variable to compute the RHS
602  Array<OneD, Array<OneD, NekDouble> > WeakAdv(nVariables);
603  WeakAdv[0] = Array<OneD, NekDouble> (ncoeffs*nVariables);
604  for(i = 1; i < nVariables; ++i)
605  {
606  WeakAdv[i] = WeakAdv[i-1] + ncoeffs;
607  }
608 
609  // Currently assume velocity field is time independent and does not therefore
610  // need extrapolating.
611  // RHS computation using the advection base class
612  m_advObject->Advect(nVariables, m_fields, m_velocity,
613  inarray, outarray, time);
614 
615  for(i = 0; i < nVariables; ++i)
616  {
617  m_fields[i]->IProductWRTBase(outarray[i],WeakAdv[i]);
618  // negation requried due to sign of DoAdvection term to be consistent
619  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
620  }
621 
622  AddAdvectionPenaltyFlux(m_velocity, inarray, WeakAdv);
623 
624 
625  /// Operations to compute the RHS
626  for(i = 0; i < nVariables; ++i)
627  {
628  // Negate the RHS
629  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
630 
631  /// Multiply the flux by the inverse of the mass matrix
632  m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
633 
634  /// Store in outarray the physical values of the RHS
635  m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
636  }
637  }
638 
639  /**
640  * Projection used by SubStepAdvance time integration
641  */
643  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
644  Array<OneD, Array<OneD, NekDouble> > &outarray,
645  const NekDouble time)
646  {
647  ASSERTL1(inarray.num_elements() == outarray.num_elements(),"Inarray and outarray of different sizes ");
648 
649  for(int i = 0; i < inarray.num_elements(); ++i)
650  {
651  Vmath::Vcopy(inarray[i].num_elements(),inarray[i],1,outarray[i],1);
652  }
653  }
654 
656  const Array<OneD, const Array<OneD, NekDouble> > &velfield,
657  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
658  Array<OneD, Array<OneD, NekDouble> > &Outarray)
659  {
660  ASSERTL1(physfield.num_elements() == Outarray.num_elements(),
661  "Physfield and outarray are of different dimensions");
662 
663  int i;
664 
665  /// Number of trace points
666  int nTracePts = m_fields[0]->GetTrace()->GetNpoints();
667 
668  /// Forward state array
669  Array<OneD, NekDouble> Fwd(3*nTracePts);
670 
671  /// Backward state array
672  Array<OneD, NekDouble> Bwd = Fwd + nTracePts;
673 
674  /// upwind numerical flux state array
675  Array<OneD, NekDouble> numflux = Bwd + nTracePts;
676 
677  /// Normal velocity array
678  Array<OneD, NekDouble> Vn = GetNormalVel(velfield);
679 
680  for(i = 0; i < physfield.num_elements(); ++i)
681  {
682  /// Extract forwards/backwards trace spaces
683  /// Note: Needs to have correct i value to get boundary conditions
684  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
685 
686  /// Upwind between elements
687  m_fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux);
688 
689  /// Construct difference between numflux and Fwd,Bwd
690  Vmath::Vsub(nTracePts, numflux, 1, Fwd, 1, Fwd, 1);
691  Vmath::Vsub(nTracePts, numflux, 1, Bwd, 1, Bwd, 1);
692 
693  /// Calculate the numerical fluxes multipling Fwd, Bwd and
694  /// numflux by the normal advection velocity
695  Vmath::Vmul(nTracePts, Fwd, 1, Vn, 1, Fwd, 1);
696  Vmath::Vmul(nTracePts, Bwd, 1, Vn, 1, Bwd, 1);
697 
698  m_fields[0]->AddFwdBwdTraceIntegral(Fwd,Bwd,Outarray[i]);
699  }
700  }
701 
702 
704  const Array<OneD, Array<OneD,NekDouble> > inarray)
705  {
706 
707  int n_points_0 = m_fields[0]->GetExp(0)->GetTotPoints();
708  int n_element = m_fields[0]->GetExpSize();
709  int nvel = inarray.num_elements();
710  int cnt;
711 
712  ASSERTL0(nvel >= 2, "Method not implemented for 1D");
713 
714  NekDouble pntVelocity;
715 
716  // Getting the standard velocity vector on the 2D normal space
717  Array<OneD, Array<OneD, NekDouble> > stdVelocity(nvel);
718  Array<OneD, NekDouble> maxV(n_element, 0.0);
720 
721  for (int i = 0; i < nvel; ++i)
722  {
723  stdVelocity[i] = Array<OneD, NekDouble>(n_points_0);
724  }
725 
726  if (nvel == 2)
727  {
728  cnt = 0.0;
729  for (int el = 0; el < n_element; ++el)
730  {
731  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
732  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
733 
734  // reset local space if necessary
735  if(n_points != n_points_0)
736  {
737  for (int j = 0; j < nvel; ++j)
738  {
739  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
740  }
741  n_points_0 = n_points;
742  }
743 
745  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
746 
747  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
749  {
750  for (int i = 0; i < n_points; i++)
751  {
752  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
753  + gmat[2][i]*inarray[1][i+cnt];
754 
755  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
756  + gmat[3][i]*inarray[1][i+cnt];
757  }
758  }
759  else
760  {
761  for (int i = 0; i < n_points; i++)
762  {
763  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
764  + gmat[2][0]*inarray[1][i+cnt];
765 
766  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
767  + gmat[3][0]*inarray[1][i+cnt];
768  }
769  }
770 
771  cnt += n_points;
772 
773 
774  for (int i = 0; i < n_points; i++)
775  {
776  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
777  + stdVelocity[1][i]*stdVelocity[1][i];
778 
779  if (pntVelocity>maxV[el])
780  {
781  maxV[el] = pntVelocity;
782  }
783  }
784  maxV[el] = sqrt(maxV[el]);
785  }
786  }
787  else
788  {
789  cnt = 0;
790  for (int el = 0; el < n_element; ++el)
791  {
792 
793  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
794  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
795 
796  // reset local space if necessary
797  if(n_points != n_points_0)
798  {
799  for (int j = 0; j < nvel; ++j)
800  {
801  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
802  }
803  n_points_0 = n_points;
804  }
805 
807  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
808 
809  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
811  {
812  for (int i = 0; i < n_points; i++)
813  {
814  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
815  + gmat[3][i]*inarray[1][i+cnt]
816  + gmat[6][i]*inarray[2][i+cnt];
817 
818  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
819  + gmat[4][i]*inarray[1][i+cnt]
820  + gmat[7][i]*inarray[2][i+cnt];
821 
822  stdVelocity[2][i] = gmat[2][i]*inarray[0][i+cnt]
823  + gmat[5][i]*inarray[1][i+cnt]
824  + gmat[8][i]*inarray[2][i+cnt];
825  }
826  }
827  else
828  {
829  for (int i = 0; i < n_points; i++)
830  {
831  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
832  + gmat[3][0]*inarray[1][i+cnt]
833  + gmat[6][0]*inarray[2][i+cnt];
834 
835  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
836  + gmat[4][0]*inarray[1][i+cnt]
837  + gmat[7][0]*inarray[2][i+cnt];
838 
839  stdVelocity[2][i] = gmat[2][0]*inarray[0][i+cnt]
840  + gmat[5][0]*inarray[1][i+cnt]
841  + gmat[8][0]*inarray[2][i+cnt];
842  }
843  }
844 
845  cnt += n_points;
846 
847  for (int i = 0; i < n_points; i++)
848  {
849  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
850  + stdVelocity[1][i]*stdVelocity[1][i]
851  + stdVelocity[2][i]*stdVelocity[2][i];
852 
853  if (pntVelocity > maxV[el])
854  {
855  maxV[el] = pntVelocity;
856  }
857  }
858 
859  maxV[el] = sqrt(maxV[el]);
860  //cout << maxV[el]*maxV[el] << endl;
861  }
862  }
863 
864  return maxV;
865  }
866 }
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:161
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:220
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:857
Array< OneD, Array< OneD, NekDouble > > m_velocity
boost::shared_ptr< TimeIntegrationWrapper > TimeIntegrationWrapperSharedPtr
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.
STL namespace.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:251
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:53
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:199
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:382
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:329
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:359
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:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Geometry is curved or has non-constant factors.
virtual void v_InitObject()
Initialise the object.
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
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
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.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215