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