Nektar++
SubSteppingExtrapolate.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: SubSteppingExtrapolate.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: Abstract base class for SubSteppingExtrapolate.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
39 namespace Nektar
40 {
41  /**
42  * Registers the class with the Factory.
43  */
45  "SubStepping",
47  "SubStepping");
48 
53  const Array<OneD, int> pVel,
54  const SolverUtils::AdvectionSharedPtr advObject)
55  : Extrapolate(pSession,pFields,pPressure,pVel,advObject)
56  {
57  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
58  m_session->LoadParameter("CFL", m_cflSafetyFactor, 0.5);
59  m_session->LoadParameter("MinSubSteps", minsubsteps,0);
60  }
61 
63  {
64  }
65 
67  int intMethod,
68  const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
69  {
70  int i;
71 
72  // Set to 1 for first step and it will then be increased in
73  // time advance routines
74  switch(intMethod)
75  {
78  {
80 
81  // Fields for linear interpolation
83  int ntotpts = m_fields[0]->GetTotPoints();
84  m_previousVelFields[0] = Array<OneD, NekDouble>(2*m_fields.num_elements()*ntotpts);
85  for(i = 1; i < 2*m_fields.num_elements(); ++i)
86  {
87  m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
88  }
89 
90  }
91  break;
93  {
95 
96  int nvel = m_velocity.num_elements();
97 
98  // Fields for quadratic interpolation
100 
101  int ntotpts = m_fields[0]->GetTotPoints();
102  m_previousVelFields[0] = Array<OneD, NekDouble>(3*nvel*ntotpts);
103  for(i = 1; i < 3*nvel; ++i)
104  {
105  m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
106  }
107 
108  }
109  break;
110  default:
111  ASSERTL0(0,"Integration method not suitable: Options include BackwardEuler or BDFImplicitOrder1");
112  break;
113  }
114 
115  m_intSteps = IntegrationScheme->GetIntegrationSteps();
116 
117  // set explicit time-integration class operators
120  }
121 
122  /**
123  * Explicit Advection terms used by SubStepAdvance time integration
124  */
126  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
127  Array<OneD, Array<OneD, NekDouble> > &outarray,
128  const NekDouble time)
129  {
130  int i;
131  int nVariables = inarray.num_elements();
132  int nQuadraturePts = inarray[0].num_elements();
133 
134  /// Get the number of coefficients
135  int ncoeffs = m_fields[0]->GetNcoeffs();
136 
137  /// Define an auxiliary variable to compute the RHS
138  Array<OneD, Array<OneD, NekDouble> > WeakAdv(nVariables);
139  WeakAdv[0] = Array<OneD, NekDouble> (ncoeffs*nVariables);
140  for(i = 1; i < nVariables; ++i)
141  {
142  WeakAdv[i] = WeakAdv[i-1] + ncoeffs;
143  }
144 
145  Array<OneD, Array<OneD, NekDouble> > Velfields(m_velocity.num_elements());
146  Array<OneD, int> VelIds(m_velocity.num_elements());
147 
148  Velfields[0] = Array<OneD, NekDouble> (nQuadraturePts*m_velocity.num_elements());
149 
150  for(i = 1; i < m_velocity.num_elements(); ++i)
151  {
152  Velfields[i] = Velfields[i-1] + nQuadraturePts;
153  }
154 
155  SubStepExtrapolateField(fmod(time,m_timestep), Velfields);
156 
157  m_advObject->Advect(m_velocity.num_elements(), m_fields, Velfields, inarray, outarray, time);
158 
159  for(i = 0; i < nVariables; ++i)
160  {
161  m_fields[i]->IProductWRTBase(outarray[i],WeakAdv[i]);
162  // negation requried due to sign of DoAdvection term to be consistent
163  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
164  }
165 
166  AddAdvectionPenaltyFlux(Velfields, inarray, WeakAdv);
167 
168  /// Operations to compute the RHS
169  for(i = 0; i < nVariables; ++i)
170  {
171  // Negate the RHS
172  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
173 
174  /// Multiply the flux by the inverse of the mass matrix
175  m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
176 
177  /// Store in outarray the physical values of the RHS
178  m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
179  }
180 
181  //add the force
182  /*if(m_session->DefinesFunction("BodyForce"))
183  {
184  if(m_SingleMode || m_HalfMode)
185  {
186  for(int i = 0; i < m_nConvectiveFields; ++i)
187  {
188  m_forces[i]->SetWaveSpace(true);
189  m_forces[i]->BwdTrans(m_forces[i]->GetCoeffs(),
190  m_forces[i]->UpdatePhys());
191  }
192  }
193 
194  int nqtot = m_fields[0]->GetTotPoints();
195  for(int i = 0; i < m_nConvectiveFields; ++i)
196  {
197  Vmath::Vadd(nqtot,outarray[i],1,(m_forces[i]->GetPhys()),1,outarray[i],1);
198  }
199  }*/
200  }
201 
202  /**
203  * Projection used by SubStepAdvance time integration
204  */
206  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
207  Array<OneD, Array<OneD, NekDouble> > &outarray,
208  const NekDouble time)
209  {
210  ASSERTL1(inarray.num_elements() == outarray.num_elements(),"Inarray and outarray of different sizes ");
211 
212  for(int i = 0; i < inarray.num_elements(); ++i)
213  {
214  Vmath::Vcopy(inarray[i].num_elements(),inarray[i],1,outarray[i],1);
215  }
216  }
217 
218 
219  /**
220  *
221  */
223  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
224  const NekDouble Aii_Dt,
225  NekDouble kinvis)
226  {
227  int nConvectiveFields =m_fields.num_elements()-1;
228  Array<OneD, Array<OneD, NekDouble> > velfields(nConvectiveFields);
229 
230  for(int i = 0; i < nConvectiveFields; ++i)
231  {
232  velfields[i] = m_fields[m_velocity[i]]->GetPhys();
233  }
234 
235  EvaluatePressureBCs(velfields,inarray,kinvis);
236 
237  AddDuDt(inarray,Aii_Dt);
238  }
239 
240  /**
241  *
242  */
244  const int nstep)
245  {
246  int i,n;
247  int nvel = m_velocity.num_elements();
248  int npts = m_fields[0]->GetTotPoints();
249 
250  // rotate fields
251  int nblocks = m_previousVelFields.num_elements()/nvel;
253 
254  for(n = 0; n < nvel; ++n)
255  {
256  save = m_previousVelFields[(nblocks-1)*nvel+n];
257 
258  for(i = nblocks-1; i > 0; --i)
259  {
260  m_previousVelFields[i*nvel+n] = m_previousVelFields[(i-1)*nvel+n];
261  }
262 
263  m_previousVelFields[n] = save;
264  }
265 
266  // Put previous field
267  for(i = 0; i < nvel; ++i)
268  {
269  m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
270  m_fields[m_velocity[i]]->UpdatePhys());
271  Vmath::Vcopy(npts,m_fields[m_velocity[i]]->GetPhys(),1,
272  m_previousVelFields[i],1);
273  }
274 
275  if(nstep == 0)// initialise all levels with first field
276  {
277  for(n = 0; n < nvel; ++n)
278  {
279  for(i = 1; i < nblocks; ++i)
280  {
281  Vmath::Vcopy(npts,m_fields[m_velocity[n]]->GetPhys(),1,
282  m_previousVelFields[i*nvel+n],1);
283 
284  }
285  }
286  }
287  }
288 
289  /**
290  *
291  */
293  const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln,
294  int nstep,
295  NekDouble time)
296  {
297  int n;
298  int nsubsteps;
299 
300  NekDouble dt;
301 
302  Array<OneD, Array<OneD, NekDouble> > fields, velfields;
303 
304  static int ncalls = 1;
305  int nint = min(ncalls++, m_intSteps);
306 
307  Array<OneD, NekDouble> CFL(m_fields[0]->GetExpSize(),
309  //this needs to change
310  m_comm = m_fields[0]->GetComm()->GetRowComm();
311 
312  // Get the proper time step with CFL control
313  dt = GetSubstepTimeStep();
314 
315  nsubsteps = (m_timestep > dt)? ((int)(m_timestep/dt)+1):1;
316  nsubsteps = max(minsubsteps, nsubsteps);
317 
318  dt = m_timestep/nsubsteps;
319 
320  if (m_infosteps && !((nstep+1)%m_infosteps) && m_comm->GetRank() == 0)
321  {
322  cout << "Sub-integrating using "<< nsubsteps
323  << " steps over Dt = " << m_timestep
324  << " (SubStep CFL=" << m_cflSafetyFactor << ")"<< endl;
325  }
326 
327  for (int m = 0; m < nint; ++m)
328  {
329  // We need to update the fields held by the m_integrationSoln
330  fields = integrationSoln->UpdateSolutionVector()[m];
331 
332  // Initialise NS solver which is set up to use a GLM method
333  // with calls to EvaluateAdvection_SetPressureBCs and
334  // SolveUnsteadyStokesSystem
336  SubIntegrationSoln = m_subStepIntegrationScheme->
337  InitializeScheme(dt, fields, time, m_subStepIntegrationOps);
338 
339  for(n = 0; n < nsubsteps; ++n)
340  {
341  fields = m_subStepIntegrationScheme->TimeIntegrate(
342  n, dt, SubIntegrationSoln,
344  }
345 
346  // Reset time integrated solution in m_integrationSoln
347  integrationSoln->SetSolVector(m,fields);
348  }
349  }
350 
351 
352  /**
353  *
354  */
356  {
357  int n_element = m_fields[0]->GetExpSize();
358 
359  const Array<OneD, int> ExpOrder=m_fields[0]->EvalBasisNumModesMaxPerExp();
360  Array<OneD, int> ExpOrderList (n_element, ExpOrder);
361 
362  const NekDouble cLambda = 0.2; // Spencer book pag. 317
363 
364  Array<OneD, NekDouble> tstep (n_element, 0.0);
365  Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
366  Array<OneD, Array<OneD, NekDouble> > velfields(m_velocity.num_elements());
367 
368  for(int i = 0; i < m_velocity.num_elements(); ++i)
369  {
370  velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
371  }
372  stdVelocity = GetMaxStdVelocity(velfields);
373 
374  for(int el = 0; el < n_element; ++el)
375  {
376  tstep[el] = m_cflSafetyFactor /
377  (stdVelocity[el] * cLambda *
378  (ExpOrder[el]-1) * (ExpOrder[el]-1));
379  }
380 
381  NekDouble TimeStep = Vmath::Vmin(n_element, tstep, 1);
382  m_comm->AllReduce(TimeStep,LibUtilities::ReduceMin);
383 
384  return TimeStep;
385  }
386 
387 
388 
389 
391  const Array<OneD, const Array<OneD, NekDouble> > &velfield,
392  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
393  Array<OneD, Array<OneD, NekDouble> > &Outarray)
394  {
395  ASSERTL1(
396  physfield.num_elements() == Outarray.num_elements(),
397  "Physfield and outarray are of different dimensions");
398 
399  int i;
400 
401  /// Number of trace points
402  int nTracePts = m_fields[0]->GetTrace()->GetNpoints();
403 
404  /// Number of spatial dimensions
405  int nDimensions = m_curl_dim;
406 
408 
409  for(i = 0; i < nDimensions; ++i)
410  {
411  traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
412  }
413 
414  //trace normals
415  m_fields[0]->GetTrace()->GetNormals(traceNormals);
416 
417  /// Forward state array
418  Array<OneD, NekDouble> Fwd(3*nTracePts);
419 
420  /// Backward state array
421  Array<OneD, NekDouble> Bwd = Fwd + nTracePts;
422 
423  /// upwind numerical flux state array
424  Array<OneD, NekDouble> numflux = Bwd + nTracePts;
425 
426  /// Normal velocity array
427  Array<OneD, NekDouble> Vn (nTracePts, 0.0);
428 
429  // Extract velocity field along the trace space and multiply by trace normals
430  for(i = 0; i < nDimensions; ++i)
431  {
432  m_fields[0]->ExtractTracePhys(m_fields[m_velocity[i]]->GetPhys(), Fwd);
433  Vmath::Vvtvp(nTracePts, traceNormals[i], 1, Fwd, 1, Vn, 1, Vn, 1);
434  }
435 
436  for(i = 0; i < physfield.num_elements(); ++i)
437  {
438  /// Extract forwards/backwards trace spaces
439  /// Note: Needs to have correct i value to get boundary conditions
440  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
441 
442  /// Upwind between elements
443  m_fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux);
444 
445  /// Construct difference between numflux and Fwd,Bwd
446  Vmath::Vsub(nTracePts, numflux, 1, Fwd, 1, Fwd, 1);
447  Vmath::Vsub(nTracePts, numflux, 1, Bwd, 1, Bwd, 1);
448 
449  /// Calculate the numerical fluxes multipling Fwd, Bwd and
450  /// numflux by the normal advection velocity
451  Vmath::Vmul(nTracePts, Fwd, 1, Vn, 1, Fwd, 1);
452  Vmath::Vmul(nTracePts, Bwd, 1, Vn, 1, Bwd, 1);
453 
454  m_fields[0]->AddFwdBwdTraceIntegral(Fwd,Bwd,Outarray[i]);
455  }
456  }
457 
458  /**
459  * Extrapolate field using equally time spaced field un,un-1,un-2, (at
460  * dt intervals) to time n+t at order Ord
461  */
463  NekDouble toff,
464  Array< OneD, Array<OneD, NekDouble> > &ExtVel)
465  {
466  int npts = m_fields[0]->GetTotPoints();
467  int nvel = m_velocity.num_elements();
468  int i,j;
470 
471  int ord = m_intSteps;
472 
473  // calculate Lagrange interpolants
474  Vmath::Fill(4,1.0,l,1);
475 
476  for(i = 0; i <= ord; ++i)
477  {
478  for(j = 0; j <= ord; ++j)
479  {
480  if(i != j)
481  {
482  l[i] *= (j*m_timestep+toff);
483  l[i] /= (j*m_timestep-i*m_timestep);
484  }
485  }
486  }
487 
488  for(i = 0; i < nvel; ++i)
489  {
490  Vmath::Smul(npts,l[0],m_previousVelFields[i],1,ExtVel[i],1);
491 
492  for(j = 1; j <= ord; ++j)
493  {
494  Blas::Daxpy(npts,l[j],m_previousVelFields[j*nvel+i],1,
495  ExtVel[i],1);
496  }
497  }
498  }
499 
500  /**
501  *
502  */
504  int HBCdata,
505  NekDouble kinvis,
507  Array<OneD, const NekDouble> &Advection)
508  {
509  Vmath::Smul(HBCdata,-kinvis,Q,1,Q,1);
510  }
511 
512  /**
513  *
514  */
516  const Array<OneD, const Array<OneD, NekDouble> > &N,
517  NekDouble Aii_Dt)
518  {
519  switch(m_velocity.num_elements())
520  {
521  case 1:
522  ASSERTL0(
523  false,
524  "Velocity correction scheme not designed to have just one velocity component");
525  break;
526  case 2:
527  AddDuDt2D(N,Aii_Dt);
528  break;
529  case 3:
530  AddDuDt3D(N,Aii_Dt);
531  break;
532  }
533  }
534 
535  /**
536  *
537  */
539  const Array<OneD, const Array<OneD, NekDouble> > &N,
540  NekDouble Aii_Dt)
541  {
542  int i,n;
543  ASSERTL0(m_velocity.num_elements() == 2," Routine currently only set up for 2D");
544 
545  int pindex=2;
546 
548  Array<OneD, MultiRegions::ExpListSharedPtr> PBndExp,UBndExp,VBndExp;
549 
550  PBndConds = m_fields[pindex]->GetBndConditions();
551  PBndExp = m_fields[pindex]->GetBndCondExpansions();
552 
553  UBndExp = m_fields[m_velocity[0]]->GetBndCondExpansions();
554  VBndExp = m_fields[m_velocity[1]]->GetBndCondExpansions();
555 
558 
559 
560  int cnt,elmtid,nq,offset, boundary,ncoeffs;
561 
565  Array<OneD, NekDouble> Nu,Nv,Ptmp;
566 
567  for(cnt = n = 0; n < PBndConds.num_elements(); ++n)
568  {
569  std::string type = PBndConds[n]->GetUserDefined();
570 
571  if(boost::iequals(type,"H"))
572  {
573  for(i = 0; i < PBndExp[n]->GetExpSize(); ++i,cnt++)
574  {
575  // find element and edge of this expansion.
576  elmtid = m_pressureBCtoElmtID[cnt];
577  elmt = m_fields[0]->GetExp(elmtid);
578  offset = m_fields[0]->GetPhys_Offset(elmtid);
579 
580  Nu = N[0] + offset;
581  Nv = N[1] + offset;
582 
583  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion1D> (PBndExp[n]->GetExp(i));
584  nq = Pbc->GetTotPoints();
585  ncoeffs = Pbc->GetNcoeffs();
586  boundary = m_pressureBCtoTraceID[cnt];
587 
588  // Get velocity bc
589  UBndExp[n]->GetExp(i)->BwdTrans(UBndExp[n]->GetCoeffs() + UBndExp[n]->GetCoeff_Offset(i),ubc);
590  VBndExp[n]->GetExp(i)->BwdTrans(VBndExp[n]->GetCoeffs() + VBndExp[n]->GetCoeff_Offset(i),vbc);
591 
592 
593  // Get edge values and put into Nu,Nv
594  elmt->GetEdgePhysVals(boundary,Pbc,Nu,N1);
595  elmt->GetEdgePhysVals(boundary,Pbc,Nv,N2);
596 
597 
598  // Take different as Forward Euler but N1,N2
599  // actually contain the integration of the
600  // previous steps from the time integration
601  // scheme.
602  Vmath::Vsub(nq,ubc,1,N1,1,ubc,1);
603  Vmath::Vsub(nq,vbc,1,N2,1,vbc,1);
604 
605 
606  // Divide by aii_Dt to get correct Du/Dt. This is
607  // because all coefficients in the integration
608  // scheme are normalised so u^{n+1} has unit
609  // coefficient and N is already multiplied by
610  // local coefficient when taken from integration
611  // scheme
612  Blas::Dscal(nq,1.0/Aii_Dt,&ubc[0],1);
613  Blas::Dscal(nq,1.0/Aii_Dt,&vbc[0],1);
614 
615  // subtrace off du/dt derivative
616  Pbc->NormVectorIProductWRTBase(ubc,vbc,Pvals);
617 
618  Vmath::Vsub(ncoeffs,Ptmp = PBndExp[n]->UpdateCoeffs()
619  +PBndExp[n]->GetCoeff_Offset(i),1, Pvals,1,
620  Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),1);
621  }
622  }
623  else // No High order
624  {
625  cnt += PBndExp[n]->GetExpSize();
626  }
627  }
628  }
629 
630  /**
631  *
632  */
634  const Array<OneD, const Array<OneD, NekDouble> > &N,
635  NekDouble Aii_Dt)
636  {
637  int i,n;
638  ASSERTL0(m_velocity.num_elements() == 3," Routine currently only set up for 3D");
639  int pindex=3;
640 
642  Array<OneD, MultiRegions::ExpListSharedPtr> PBndExp,UBndExp,VBndExp,WBndExp;
643 
644  PBndConds = m_fields[pindex]->GetBndConditions();
645  PBndExp = m_fields[pindex]->GetBndCondExpansions();
646 
647  UBndExp = m_fields[m_velocity[0]]->GetBndCondExpansions();
648  VBndExp = m_fields[m_velocity[1]]->GetBndCondExpansions();
649  WBndExp = m_fields[m_velocity[2]]->GetBndCondExpansions();
650 
653 
654  int cnt,elmtid,nq,offset, boundary,ncoeffs;
655 
659  wbc(m_pressureBCsMaxPts);
661  Array<OneD, NekDouble> Nu,Nv,Nw,Ptmp;
662 
663  for(cnt = n = 0; n < PBndConds.num_elements(); ++n)
664  {
665  std::string type = PBndConds[n]->GetUserDefined();
666 
667  if(boost::iequals(type,"H"))
668  {
669  for(i = 0; i < PBndExp[n]->GetExpSize(); ++i,cnt++)
670  {
671  // find element and face of this expansion.
672  elmtid = m_pressureBCtoElmtID[cnt];
673  elmt = m_fields[0]->GetExp(elmtid);
674  offset = m_fields[0]->GetPhys_Offset(elmtid);
675 
676  Nu = N[0] + offset;
677  Nv = N[1] + offset;
678  Nw = N[2] + offset;
679 
680  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (PBndExp[n]->GetExp(i));
681  nq = Pbc->GetTotPoints();
682  ncoeffs = Pbc->GetNcoeffs();
683  boundary = m_pressureBCtoTraceID[cnt];
684 
685  // Get velocity bc
686  UBndExp[n]->GetExp(i)->BwdTrans(UBndExp[n]->GetCoeffs() +
687  UBndExp[n]->GetCoeff_Offset(i),ubc);
688  VBndExp[n]->GetExp(i)->BwdTrans(VBndExp[n]->GetCoeffs() +
689  VBndExp[n]->GetCoeff_Offset(i),vbc);
690  WBndExp[n]->GetExp(i)->BwdTrans(WBndExp[n]->GetCoeffs() +
691  WBndExp[n]->GetCoeff_Offset(i),wbc);
692 
693  // Get edge values and put into N1,N2,N3
694  elmt->GetFacePhysVals(boundary,Pbc,Nu,N1);
695  elmt->GetFacePhysVals(boundary,Pbc,Nv,N2);
696  elmt->GetFacePhysVals(boundary,Pbc,Nw,N3);
697 
698 
699  // Take different as Forward Euler but N1,N2,N3
700  // actually contain the integration of the
701  // previous steps from the time integration
702  // scheme.
703  Vmath::Vsub(nq,ubc,1,N1,1,ubc,1);
704  Vmath::Vsub(nq,vbc,1,N2,1,vbc,1);
705  Vmath::Vsub(nq,wbc,1,N3,1,wbc,1);
706 
707  // Divide by aii_Dt to get correct Du/Dt. This is
708  // because all coefficients in the integration
709  // scheme are normalised so u^{n+1} has unit
710  // coefficient and N is already multiplied by
711  // local coefficient when taken from integration
712  // scheme
713  Blas::Dscal(nq,1.0/Aii_Dt,&ubc[0],1);
714  Blas::Dscal(nq,1.0/Aii_Dt,&vbc[0],1);
715  Blas::Dscal(nq,1.0/Aii_Dt,&wbc[0],1);
716 
717  // subtrace off du/dt derivative
718  Pbc->NormVectorIProductWRTBase(ubc,vbc,wbc,Pvals);
719 
720  Vmath::Vsub(
721  ncoeffs,
722  Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),
723  1,Pvals,1,
724  Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),
725  1);
726  }
727  }
728  else // No High order
729  {
730  cnt += PBndExp[n]->GetExpSize();
731  }
732  }
733  }
734 }
735 
void SubStepProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void SubStepExtrapolateField(NekDouble toff, Array< OneD, Array< OneD, NekDouble > > &ExtVel)
LibUtilities::TimeIntegrationWrapperSharedPtr m_subStepIntegrationScheme
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void AddDuDt3D(const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble Aii_Dt)
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:170
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition: Extrapolate.h:181
BDF multi-step scheme of order 1 (implicit)
int m_pressureBCsMaxPts
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:203
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
virtual void v_SubStepSetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, NekDouble Aii_Dt, NekDouble kinvis)
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:48
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:848
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
boost::shared_ptr< TimeIntegrationWrapper > TimeIntegrationWrapperSharedPtr
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:158
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
void AddAdvectionPenaltyFlux(const Array< OneD, const Array< OneD, NekDouble > > &velfield, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &outarray)
void AddDuDt(const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble Aii_Dt)
LibUtilities::CommSharedPtr m_comm
Definition: Extrapolate.h:172
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Extrapolate.h:174
virtual void v_SubStepSaveFields(int nstep)
Array< OneD, Array< OneD, NekDouble > > m_previousVelFields
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
virtual void v_SubSteppingTimeIntegration(int intMethod, const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
Array< OneD, int > m_pressureBCtoElmtID
Id of element to which pressure boundary condition belongs.
Definition: Extrapolate.h:227
static std::string className
Name of class.
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
Array< OneD, NekDouble > GetMaxStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:209
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
static std::string npts
Definition: InputFld.cpp:43
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
virtual void v_SubStepAdvance(const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln, int nstep, NekDouble time)
double NekDouble
TimeIntegrationWrapperFactory & GetTimeIntegrationWrapperFactory()
Array< OneD, int > m_pressureBCtoTraceID
Id of edge (2D) or face (3D) to which pressure boundary condition belongs.
Definition: Extrapolate.h:230
SubSteppingExtrapolate(const LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields, MultiRegions::ExpListSharedPtr pPressure, const Array< OneD, int > pVel, const SolverUtils::AdvectionSharedPtr advObject)
void SubStepAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble 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
virtual void v_MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
BDF multi-step scheme of order 2 (implicit)
static ExtrapolateSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, MultiRegions::ExpListSharedPtr &pPressure, const Array< OneD, int > &pVel, const SolverUtils::AdvectionSharedPtr &advObject)
Creates an instance of this class.
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
void AddDuDt2D(const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble Aii_Dt)
SolverUtils::AdvectionSharedPtr m_advObject
Definition: Extrapolate.h:183
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:188
void EvaluatePressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
Definition: Extrapolate.cpp:90
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#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:1038
NekDouble m_timestep
Definition: Extrapolate.h:211
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
LibUtilities::TimeIntegrationSchemeOperators m_subStepIntegrationOps
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215