Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Extrapolate.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Extrapolate.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 Extrapolate.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
39 using namespace std;
40 
41 namespace Nektar
42 {
43  NekDouble Extrapolate::StifflyStable_Betaq_Coeffs[3][3] = {
44  { 1.0, 0.0, 0.0},{ 2.0, -1.0, 0.0},{ 3.0, -3.0, 1.0}};
45  NekDouble Extrapolate::StifflyStable_Alpha_Coeffs[3][3] = {
46  { 1.0, 0.0, 0.0},{ 2.0, -0.5, 0.0},{ 3.0, -1.5, 1.0/3.0}};
47  NekDouble Extrapolate::StifflyStable_Gamma0_Coeffs[3] = {
48  1.0, 1.5, 11.0/6.0};
49 
51  {
52  typedef Loki::SingletonHolder<ExtrapolateFactory,
53  Loki::CreateUsingNew,
54  Loki::NoDestroy,
55  Loki::SingleThreaded > Type;
56  return Type::Instance();
57  }
58 
59  Extrapolate::Extrapolate(
63  const Array<OneD, int> pVel,
64  const SolverUtils::AdvectionSharedPtr advObject)
65  : m_session(pSession),
66  m_fields(pFields),
67  m_pressure(pPressure),
68  m_velocity(pVel),
69  m_advObject(advObject)
70  {
71  m_session->LoadParameter("TimeStep", m_timestep, 0.01);
72  m_comm = m_session->GetComm();
73  }
74 
76  {
77  }
78 
79  std::string Extrapolate::def =
81  "StandardExtrapolate", "StandardExtrapolate");
82 
83  /**
84  *
85  */
87  {
88  if(m_numHBCDof)
89  {
90  // Update velocity BF at n+1 (actually only needs doing if
91  // velocity is time dependent on HBCs)
93 
94  //Calculate acceleration term at level n based on previous steps
96 
97  // Subtract acceleration term off m_pressureHBCs[nlevels-1]
99  m_iprodnormvel[m_intSteps], 1,
100  m_pressureHBCs[m_intSteps-1], 1,
101  m_pressureHBCs[m_intSteps-1], 1);
102  }
103  }
104 
105  /**
106  *
107  */
109  {
110  if(m_numHBCDof)
111  {
112  int order = std::min(m_pressureCalls,m_intSteps);
113 
114  // Update velocity BF at n+1 (actually only needs doing if
115  // velocity is time dependent on HBCs)
117 
118  // Subtract acceleration term off m_pressureHBCs[nlevels-1]
121  m_iprodnormvel[0], 1,
123  m_pressureHBCs[m_intSteps-1], 1);
124  }
125  }
126 
127  /**
128  * Unified routine for calculation high-oder terms
129  */
131  const Array<OneD, const Array<OneD, NekDouble> > &fields,
132  const Array<OneD, const Array<OneD, NekDouble> > &N,
133  NekDouble kinvis)
134  {
135  int n, cnt;
136 
138 
141 
144 
146  for(n = cnt = 0; n < m_PBndConds.num_elements(); ++n)
147  {
148  // High order boundary condition;
149  if((m_hbcType[n] == eHBCNeumann)||(m_hbcType[n] == eConvectiveOBC))
150  {
151  m_fields[0]->GetBndElmtExpansion(n, BndElmtExp, false);
152  int nqb = m_PBndExp[n]->GetTotPoints();
153  int nq = BndElmtExp->GetTotPoints();
154 
155  for(int i = 0; i < m_bnd_dim; i++)
156  {
157  BndValues[i] = Array<OneD, NekDouble> (nqb,0.0);
158  }
159 
160  for(int i = 0; i < m_curl_dim; i++)
161  {
162  Q[i] = Array<OneD, NekDouble> (nq,0.0);
163  }
164 
165  // Obtaining fields on BndElmtExp
166  for(int i = 0; i < m_curl_dim; i++)
167  {
168  m_fields[0]->ExtractPhysToBndElmt(n, fields[i],Velocity[i]);
169  }
170 
171  if(N.num_elements()) // not required for some extrapolation
172  {
173  for(int i = 0; i < m_bnd_dim; i++)
174  {
175  m_fields[0]->ExtractPhysToBndElmt(n, N[i], Advection[i]);
176  }
177  }
178 
179  // CurlCurl
180  BndElmtExp->CurlCurl(Velocity, Q);
181 
182  // Mounting advection component into the high-order condition
183  for(int i = 0; i < m_bnd_dim; i++)
184  {
185  MountHOPBCs(nq, kinvis,Q[i],Advection[i]);
186  }
187 
188  Pvals = (m_pressureHBCs[m_intSteps-1]) + cnt;
189 
190  // Getting values on the boundary and filling the pressure bnd
191  // expansion. Multiplication by the normal is required
192  for(int i = 0; i < m_bnd_dim; i++)
193  {
194  m_fields[0]->ExtractElmtToBndPhys(n, Q[i],BndValues[i]);
195  }
196 
197  m_PBndExp[n]->NormVectorIProductWRTBase(BndValues, Pvals);
198 
199  // Get offset for next terms
200  cnt += m_PBndExp[n]->GetNcoeffs();
201  }
202  }
203  }
204 
205  // do nothing unless otherwise defined.
207  {
208  }
209 
210  // do nothing unless otherwise defined.
211  void Extrapolate::v_AddNormVelOnOBC(const int nbcoeffs, const int nreg,
213  {
214  }
215 
217  const Array<OneD, const Array<OneD, NekDouble> > &fields,
218  NekDouble kinvis)
219  {
220  if(!m_houtflow.get())
221  {
222  return;
223  }
224 
226 
228  int cnt = 0;
229 
230  // Evaluate robin primitive coefficient here so they can be
231  // updated whem m_int > 1 Currently not using this update
232  // since we only using u^n at outflow instead of BDF rule.
234 
235  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
236  {
237  if((m_hbcType[n] == eOBC)||(m_hbcType[n] == eConvectiveOBC))
238  {
239  // Get expansion with element on this boundary
240  m_fields[0]->GetBndElmtExpansion(n, BndElmtExp, false);
241  int nqb = m_PBndExp[n]->GetTotPoints();
242  int nq = BndElmtExp->GetTotPoints();
243 
244  // Get velocity and extrapolate
245  for(int i = 0; i < m_curl_dim; i++)
246  {
247  m_fields[0]->ExtractPhysToBndElmt(n, fields[i],
248  m_houtflow->m_outflowVel[cnt][i][m_intSteps-1]);
249  ExtrapolateArray(m_houtflow->m_outflowVel[cnt][i]);
250  Velocity[i] = m_houtflow->m_outflowVel[cnt][i][m_intSteps-1];
251 
252  }
253 
254  // Homogeneous case needs conversion to physical space
255  if ( m_fields[0]->GetWaveSpace())
256  {
257  for(int i = 0; i < m_curl_dim; i++)
258  {
259  BndElmtExp->HomogeneousBwdTrans(Velocity[i],
260  Velocity[i]);
261  }
262  BndElmtExp->SetWaveSpace(false);
263  }
264 
265  // Get normal vector
267  m_fields[0]->GetBoundaryNormals(n, normals);
268 
269  // Calculate n.gradU.n, div(U)
270  Array<OneD, NekDouble> nGradUn (nqb, 0.0);
271  Array<OneD, NekDouble> divU (nqb, 0.0);
272  Array<OneD, Array<OneD, NekDouble> > grad(m_curl_dim);
273  Array<OneD, NekDouble> bndVal (nqb, 0.0);
274  for( int i = 0; i < m_curl_dim; i++)
275  {
276  grad[i] = Array<OneD, NekDouble> (nq, 0.0);
277  }
278  for( int i = 0; i < m_curl_dim; i++)
279  {
280  if( m_curl_dim == 2)
281  {
282  BndElmtExp->PhysDeriv(Velocity[i], grad[0], grad[1]);
283  }
284  else
285  {
286  BndElmtExp->PhysDeriv(Velocity[i], grad[0], grad[1],
287  grad[2]);
288  }
289 
290  for( int j = 0; j < m_curl_dim; j++)
291  {
292  m_fields[0]->ExtractElmtToBndPhys(n, grad[j],bndVal);
293  // div(U) = gradU_ii
294  if ( i == j)
295  {
296  Vmath::Vadd(nqb , divU, 1, bndVal, 1, divU, 1);
297  }
298  // n.gradU.n = gradU_ij n_i n_j
299  Vmath::Vmul(nqb , normals[i], 1, bndVal, 1,
300  bndVal, 1);
301  Vmath::Vvtvp(nqb , normals[j], 1, bndVal, 1,
302  nGradUn, 1, nGradUn, 1);
303  }
304  }
305 
306  // Obtain u at the boundary
307  Array<OneD, Array<OneD, NekDouble> > u (m_curl_dim);
308  for( int i = 0; i < m_curl_dim; i++)
309  {
310  u[i] = Array<OneD, NekDouble> (nqb, 0.0);
311  m_fields[0]->ExtractElmtToBndPhys(n, Velocity[i],u[i]);
312  }
313 
314  // Calculate u.n and u^2
315  Array<OneD, NekDouble> un (nqb, 0.0);
316  Array<OneD, NekDouble> u2 (nqb, 0.0);
317  for( int i = 0; i < m_curl_dim; i++)
318  {
319  Vmath::Vvtvp(nqb, normals[i], 1, u[i], 1,
320  un, 1, un, 1);
321  Vmath::Vvtvp(nqb, u[i], 1, u[i], 1,
322  u2, 1, u2, 1);
323  }
324 
325  // Calculate S_0(u.n) = 0.5*(1-tanh(u.n/*U0*delta))
326  Array<OneD, NekDouble> S0 (nqb, 0.0);
327  for( int i = 0; i < nqb; i++)
328  {
329  S0[i] = 0.5*(1.0-tanh(un[i]/(m_houtflow->m_U0*m_houtflow->m_delta)));
330  }
331 
332  // Calculate E(n,u) = ((theta+alpha2)*0.5*(u^2)n +
333  // (1-theta+alpha1)*0.5*(n.u)u ) * S_0(u.n)
334  NekDouble k1 = 0.5*(m_houtflow->m_obcTheta +
335  m_houtflow->m_obcAlpha2);
336  NekDouble k2 = 0.5*(1-m_houtflow->m_obcTheta +
337  m_houtflow->m_obcAlpha1);
338 
339  Array<OneD, Array<OneD, NekDouble> > E (m_curl_dim);
340  for( int i = 0; i < m_curl_dim; i++)
341  {
342  E[i] = Array<OneD, NekDouble> (nqb, 0.0);
343 
344  Vmath::Smul(nqb, k1, u2, 1, E[i], 1);
345  Vmath::Vmul(nqb, E[i], 1, normals[i], 1, E[i], 1);
346  // Use bndVal as a temporary storage
347  Vmath::Smul(nqb, k2, un, 1, bndVal, 1);
348  Vmath::Vvtvp(nqb, u[i], 1, bndVal, 1, E[i], 1, E[i], 1);
349  Vmath::Vmul(nqb, E[i], 1, S0, 1, E[i], 1);
350  }
351 
352  // if non-zero forcing is provided we want to subtract
353  // value if we want to interpret values as being the
354  // desired pressure value. This is now precribed from
355  // the velocity forcing to be consistent with the
356  // paper except f_b = -f_b
357 
358  // Calculate (E(n,u) + f_b).n
359  Array<OneD, NekDouble> En (nqb, 0.0);
360  for( int i = 0; i < m_bnd_dim; i++)
361  {
362  // Use bndVal as temporary
363  Vmath::Vsub(nqb,E[i],1,m_houtflow->
364  m_UBndExp[i][n]->GetPhys(),
365  1, bndVal, 1);
366 
367  Vmath::Vvtvp(nqb, normals[i], 1, bndVal, 1,
368  En, 1, En, 1);
369 
370  }
371 
372  // Calculate pressure bc = kinvis*n.gradU.n - E.n + f_b.n
373  Array<OneD, NekDouble> pbc (nqb, 0.0);
374  Vmath::Svtvm( nqb, kinvis, nGradUn, 1, En, 1, pbc, 1);
375 
376  if(m_hbcType[n] == eOBC)
377  {
378 
379  if ( m_PBndExp[n]->GetWaveSpace())
380  {
381  m_PBndExp[n]->HomogeneousFwdTrans(pbc, bndVal);
382  m_PBndExp[n]->FwdTrans(bndVal,
383  m_PBndExp[n]->UpdateCoeffs());
384  }
385  else
386  {
387  m_PBndExp[n]->FwdTrans(pbc,
388  m_PBndExp[n]->UpdateCoeffs());
389  }
390  }
391  else if(m_hbcType[n] == eConvectiveOBC) // add outflow values to calculation from HBC
392  {
393  int nbcoeffs = m_PBndExp[n]->GetNcoeffs();
394  Array<OneD, NekDouble> bndCoeffs (nbcoeffs, 0.0);
395  if ( m_PBndExp[n]->GetWaveSpace())
396  {
397  m_PBndExp[n]->HomogeneousFwdTrans(pbc, bndVal);
398  m_PBndExp[n]->IProductWRTBase(bndVal,bndCoeffs);
399  }
400  else
401  {
402  m_PBndExp[n]->IProductWRTBase(pbc,bndCoeffs);
403  }
404  // Note we have the negative of what is in the Dong paper in bndVal
405  Vmath::Svtvp(nbcoeffs,m_houtflow->m_pressurePrimCoeff[n],
406  bndCoeffs, 1,m_PBndExp[n]->UpdateCoeffs(),1,
407  m_PBndExp[n]->UpdateCoeffs(),1);
408 
409  // evaluate u^n at outflow boundary for velocity BC
410  for( int i = 0; i < m_curl_dim; i++)
411  {
412  m_fields[0]->ExtractElmtToBndPhys(n,
413  m_houtflow->
414  m_outflowVel[cnt][i][0],
415  m_houtflow->
416  m_outflowVelBnd[cnt][i][m_intSteps-1]);
417 
418  EvaluateBDFArray(m_houtflow->m_outflowVelBnd[cnt][i]);
419 
420  // point u[i] to BDF evalauted value \hat{u}
421  u[i] = m_houtflow->m_outflowVelBnd[cnt][i]
422  [m_intSteps-1];
423  }
424 
425  // Add normal velocity if weak pressure
426  // formulation. In this case there is an
427  // additional \int \hat{u}.n ds on the outflow
428  // boundary since we use the inner product wrt
429  // deriv of basis in pressure solve.
430  AddNormVelOnOBC(cnt, n, u);
431  }
432 
433  // Calculate velocity boundary conditions
434  if(m_hbcType[n] == eOBC)
435  {
436  // = (pbc n - kinvis divU n)
437  Vmath::Smul(nqb, kinvis, divU, 1, divU, 1);
438  Vmath::Vsub(nqb, pbc, 1, divU, 1, bndVal, 1);
439  }
440  else if (m_hbcType[n] == eConvectiveOBC)
441  {
442  // = (-kinvis divU n)
443  Vmath::Smul(nqb, -1.0*kinvis, divU, 1, bndVal, 1);
444 
445  // pbc needs to be added after pressure solve
446  }
447 
448  for(int i = 0; i < m_curl_dim; ++i)
449  {
450  // Reuse divU -> En
451  Vmath::Vvtvp( nqb, normals[i], 1, bndVal, 1, E[i], 1,
452  divU, 1);
453  // - f_b
454  Vmath::Vsub( nqb, divU, 1,
455  m_houtflow->m_UBndExp[i][n]->GetPhys(),
456  1, divU, 1);
457  // * 1/kinvis
458  Vmath::Smul(nqb, 1.0/kinvis, divU, 1, divU, 1);
459 
460  if(m_hbcType[n] == eConvectiveOBC)
461  {
462  Vmath::Svtvp(nqb,m_houtflow->m_velocityPrimCoeff[i][n],
463  u[i], 1,divU,1,divU,1);
464  }
465 
466  if ( m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
467  {
468  m_houtflow->m_UBndExp[i][n]->HomogeneousFwdTrans(divU,
469  divU);
470  }
471 
472  m_houtflow->m_UBndExp[i][n]->IProductWRTBase(divU,
473  m_houtflow->m_UBndExp[i][n]->UpdateCoeffs());
474 
475  }
476 
477  // Get offset for next terms
478  cnt++;
479  }
480  }
481  }
482 
483 
485  {
486  if(!m_houtflow.get())
487  {
488  return;
489  }
490 
491 
492  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
493  {
494  if(m_hbcType[n] == eConvectiveOBC)
495  {
496  int nqb = m_PBndExp[n]->GetTotPoints();
497  int ncb = m_PBndExp[n]->GetNcoeffs();
498 
499  m_pressure->FillBndCondFromField(n);
500  Array<OneD, NekDouble> pbc(nqb);
501 
502  m_PBndExp[n]->BwdTrans(m_PBndExp[n]->GetCoeffs(), pbc);
503 
504  if (m_PBndExp[n]->GetWaveSpace())
505  {
506  m_PBndExp[n]->HomogeneousBwdTrans(pbc, pbc);
507  }
508 
509  Array<OneD, NekDouble> wk(nqb);
510  Array<OneD, NekDouble> wk1(ncb);
511 
512  // Get normal vector
514  m_fields[0]->GetBoundaryNormals(n, normals);
515 
516  // Add 1/kinvis * (pbc n )
517  for(int i = 0; i < m_curl_dim; ++i)
518  {
519  Vmath::Vmul(nqb, normals[i], 1, pbc, 1, wk, 1);
520 
521  Vmath::Smul(nqb, 1.0/kinvis, wk, 1, wk, 1);
522 
523  if (m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
524  {
525  m_houtflow->m_UBndExp[i][n]->
526  HomogeneousFwdTrans(wk, wk);
527  }
528  m_houtflow->m_UBndExp[i][n]->IProductWRTBase(wk,wk1);
529 
530  Vmath::Vadd(ncb, wk1,1,
531  m_houtflow->m_UBndExp[i][n]->GetCoeffs(), 1,
532  m_houtflow->m_UBndExp[i][n]->UpdateCoeffs(),1);
533 
534  }
535  }
536  }
537  }
538 
539 
541  const Array<OneD, const Array<OneD, NekDouble> > &Vel,
542  Array<OneD, NekDouble> &IProdVn)
543  {
544  int i,n,cnt;
545  Array<OneD, NekDouble> IProdVnTmp;
547 
548  for(n = cnt = 0; n < m_PBndConds.num_elements(); ++n)
549  {
550  // High order boundary condition;
551  if(m_hbcType[n] == eHBCNeumann)
552  {
553  for(i = 0; i < m_bnd_dim; ++i)
554  {
555  m_fields[0]->ExtractPhysToBnd(n, Vel[i], velbc[i]);
556  }
557  IProdVnTmp = IProdVn + cnt;
558  m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
559  cnt += m_PBndExp[n]->GetNcoeffs();
560  }
561  else if(m_hbcType[n] == eConvectiveOBC) // skip over conective OBC
562  {
563  cnt += m_PBndExp[n]->GetNcoeffs();
564  }
565  }
566  }
567 
569  {
570 
571  if(!m_HBCnumber)
572  {
573  return;
574  }
575  int i,n,cnt;
576  Array<OneD, NekDouble> IProdVnTmp;
579  for(i = 0; i < m_bnd_dim; ++i)
580  {
581  VelBndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
582  }
583 
584  for(n = cnt = 0; n < m_PBndConds.num_elements(); ++n)
585  {
586  // High order boundary condition;
587  if(m_hbcType[n] == eHBCNeumann)
588  {
589  for(i = 0; i < m_bnd_dim; ++i)
590  {
591  velbc[i] = Array<OneD, NekDouble>
592  (VelBndExp[i][n]->GetTotPoints(), 0.0);
593  VelBndExp[i][n]->SetWaveSpace(
594  m_fields[m_velocity[i]]->GetWaveSpace());
595  VelBndExp[i][n]->BwdTrans(VelBndExp[i][n]->GetCoeffs(),
596  velbc[i]);
597  }
598  IProdVnTmp = IProdVn + cnt;
599  m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
600  cnt += m_PBndExp[n]->GetNcoeffs();
601  }
602  else if(m_hbcType[n] == eConvectiveOBC)
603  {
604  // skip over convective OBC
605  cnt += m_PBndExp[n]->GetNcoeffs();
606  }
607  }
608  }
609 
610  /**
611  * Function to roll time-level storages to the next step layout.
612  * The stored data associated with the oldest time-level
613  * (not required anymore) are moved to the top, where they will
614  * be overwritten as the solution process progresses.
615  */
617  {
618  int nlevels = input.num_elements();
619 
621 
622  tmp = input[nlevels-1];
623 
624  for(int n = nlevels-1; n > 0; --n)
625  {
626  input[n] = input[n-1];
627  }
628 
629  input[0] = tmp;
630  }
631 
632 
633  /**
634  * Initialize HOBCs
635  */
637  {
638  m_PBndConds = m_pressure->GetBndConditions();
639  m_PBndExp = m_pressure->GetBndCondExpansions();
640 
641  int cnt, n;
642 
643  // Storage array for high order pressure BCs
646 
647  // Get useful values for HOBCs
648  m_HBCnumber = 0;
649  m_numHBCDof = 0;
650 
651  int outHBCnumber = 0;
652  int numOutHBCPts = 0;
653 
655  for( n = 0; n < m_PBndConds.num_elements(); ++n)
656  {
657  // High order boundary Neumann Condiiton
658  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"H"))
659  {
660  m_hbcType[n] = eHBCNeumann;
661  m_numHBCDof += m_PBndExp[n]->GetNcoeffs();
662  m_HBCnumber += m_PBndExp[n]->GetExpSize();
663  }
664 
665  // High order outflow convective condition
666  if(m_PBndConds[n]->GetBoundaryConditionType() ==
668  boost::iequals(m_PBndConds[n]->GetUserDefined(),
669  "HOutflow"))
670  {
672  m_numHBCDof += m_PBndExp[n]->GetNcoeffs();
673  m_HBCnumber += m_PBndExp[n]->GetExpSize();
674  numOutHBCPts += m_PBndExp[n]->GetTotPoints();
675  outHBCnumber++;
676  }
677  // High order outflow boundary condition;
678  else if(boost::iequals(m_PBndConds[n]->GetUserDefined(),
679  "HOutflow"))
680  {
681  m_hbcType[n] = eOBC;
682  numOutHBCPts += m_PBndExp[n]->GetTotPoints();
683  outHBCnumber++;
684  }
685  }
686 
687  m_iprodnormvel[0] = Array<OneD, NekDouble>(m_numHBCDof, 0.0);
688  for(n = 0; n < m_intSteps; ++n)
689  {
691  m_iprodnormvel[n+1] = Array<OneD, NekDouble>(m_numHBCDof, 0.0);
692  }
693 
694  m_pressureCalls = 0;
695 
696  switch(m_pressure->GetExpType())
697  {
698  case MultiRegions::e2D:
699  {
700  m_curl_dim = 2;
701  m_bnd_dim = 2;
702  }
703  break;
705  {
706  m_curl_dim = 3;
707  m_bnd_dim = 2;
708  }
709  break;
711  {
712  m_curl_dim = 3;
713  m_bnd_dim = 1;
714  }
715  break;
716  case MultiRegions::e3D:
717  {
718  m_curl_dim = 3;
719  m_bnd_dim = 3;
720  }
721  break;
722  default:
723  ASSERTL0(0,"Dimension not supported");
724  break;
725  }
726 
727  // Initialise storage for outflow HOBCs
728  if(numOutHBCPts > 0)
729  {
730  m_houtflow = MemoryManager<HighOrderOutflow>::AllocateSharedPtr(numOutHBCPts, outHBCnumber, m_curl_dim, pSession);
731 
733 
734  // set up boundary expansions link
735  for (int i = 0; i < m_curl_dim; ++i)
736  {
737  m_houtflow->m_UBndExp[i] =
738  m_fields[m_velocity[i]]->GetBndCondExpansions();
739  }
740 
741  for(n = 0, cnt = 0; n < m_PBndConds.num_elements(); ++n)
742  {
743  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"HOutflow"))
744  {
745  m_houtflow->m_outflowVel[cnt] =
746  Array<OneD, Array<OneD,
748 
749  m_houtflow->m_outflowVelBnd[cnt] =
750  Array<OneD, Array<OneD,
752 
753  m_fields[0]->GetBndElmtExpansion(n, BndElmtExp, false);
754  int nqb = m_PBndExp[n]->GetTotPoints();
755  int nq = BndElmtExp->GetTotPoints();
756  for(int j = 0; j < m_curl_dim; ++j)
757  {
758  m_houtflow->m_outflowVel[cnt][j] =
760 
761  m_houtflow->m_outflowVelBnd[cnt][j] =
763 
764  for(int k = 0; k < m_intSteps; ++k)
765  {
766  m_houtflow->m_outflowVel[cnt][j][k] =
767  Array<OneD, NekDouble>(nq,0.0);
768  m_houtflow->m_outflowVelBnd[cnt][j][k] =
769  Array<OneD, NekDouble>(nqb,0.0);
770  }
771  }
772  cnt++;
773  }
774 
775  // evaluate convective primitive coefficient if
776  // convective OBCs are used
777  if(m_hbcType[n] == eConvectiveOBC)
778  {
779  // initialise convective members of
780  // HighOrderOutflow struct
781  if(m_houtflow->m_pressurePrimCoeff.num_elements() == 0)
782  {
783  m_houtflow->m_pressurePrimCoeff =
785  (m_PBndConds.num_elements(),0.0);
786  m_houtflow->m_velocityPrimCoeff =
788 
789 
790  for(int i = 0; i < m_curl_dim; ++i)
791  {
792  m_houtflow->m_velocityPrimCoeff[i] =
794  num_elements(),0.0);
795  }
796  }
797 
798  LibUtilities::Equation coeff =
799  boost::static_pointer_cast<
801  >(m_PBndConds[n])->m_robinPrimitiveCoeff;
802 
803  // checkout equation evaluation options!!
804  m_houtflow->m_pressurePrimCoeff[n] = coeff.Evaluate();
805 
806  for (int i = 0; i < m_curl_dim; ++i)
807  {
809  UBndConds = m_fields[m_velocity[i]]->GetBndConditions();
810 
811  LibUtilities::Equation coeff1 =
812  boost::static_pointer_cast<
814  >(UBndConds[n])->m_robinPrimitiveCoeff;
815 
816 
817  m_houtflow->m_defVelPrimCoeff[i] = coeff1.GetExpression();
818 
819  ASSERTL1(UBndConds[n]->GetBoundaryConditionType()
820  == SpatialDomains::eRobin,"Require Velocity "
821  "conditions to be of Robin type when pressure"
822  "outflow is specticied as Robin Boundary type");
823 
824 
825  // checkout equation evaluation options!!
826  m_houtflow->m_velocityPrimCoeff[i][n] = coeff1.Evaluate();
827  }
828  }
829  }
830 
831  }
832  }
833 
835  {
836 
837  if((m_pressureCalls == 1) || (m_pressureCalls > m_intSteps))
838  {
839  return;
840  }
841 
842  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
843  {
844  // Get expansion with element on this boundary
845  if(m_hbcType[n] == eConvectiveOBC)
846  {
847  for (int i = 0; i < m_curl_dim; ++i)
848  {
850  UBndConds = m_fields[m_velocity[i]]->UpdateBndConditions();
851 
852  std::string primcoeff = m_houtflow->m_defVelPrimCoeff[i] + "*" +
853  boost::lexical_cast<std::string>(StifflyStable_Gamma0_Coeffs
854  [m_pressureCalls-1]);
855 
857  boost::dynamic_pointer_cast<
859 
862  m_session,rcond->m_robinFunction.GetExpression(),
863  primcoeff,
864  rcond->GetUserDefined(),
865  rcond->m_filename);
866 
867  UBndConds[n] = bcond;
868  }
869 
870  }
871  }
872  }
873 
874  /**
875  *
876  */
878  const Array<OneD, Array<OneD,NekDouble> > inarray)
879  {
880  // Checking if the problem is 2D
881  ASSERTL0(m_curl_dim >= 2, "Method not implemented for 1D");
882 
883  int n_points_0 = m_fields[0]->GetExp(0)->GetTotPoints();
884  int n_element = m_fields[0]->GetExpSize();
885  int nvel = inarray.num_elements();
886  int cnt;
887 
888  NekDouble pntVelocity;
889 
890  // Getting the standard velocity vector
891  Array<OneD, Array<OneD, NekDouble> > stdVelocity(nvel);
893  Array<OneD, NekDouble> maxV(n_element, 0.0);
895 
896  for (int i = 0; i < nvel; ++i)
897  {
898  stdVelocity[i] = Array<OneD, NekDouble>(n_points_0);
899  }
900 
901  cnt = 0.0;
902  for (int el = 0; el < n_element; ++el)
903  {
904  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
905  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
906 
907  // reset local space
908  if(n_points != n_points_0)
909  {
910  for (int j = 0; j < nvel; ++j)
911  {
912  stdVelocity[j] = Array<OneD, NekDouble>(n_points, 0.0);
913  }
914  n_points_0 = n_points;
915  }
916  else
917  {
918  for (int j = 0; j < nvel; ++j)
919  {
920  Vmath::Zero( n_points, stdVelocity[j], 1);
921  }
922  }
923 
925  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
926 
927  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
929  {
930  for(int j = 0; j < nvel; ++j)
931  {
932  for(int k = 0; k < nvel; ++k)
933  {
934  Vmath::Vvtvp( n_points, gmat[k*nvel + j], 1,
935  tmp = inarray[k] + cnt, 1,
936  stdVelocity[j], 1,
937  stdVelocity[j], 1);
938  }
939  }
940  }
941  else
942  {
943  for(int j = 0; j < nvel; ++j)
944  {
945  for(int k = 0; k < nvel; ++k)
946  {
947  Vmath::Svtvp( n_points, gmat[k*nvel + j][0],
948  tmp = inarray[k] + cnt, 1,
949  stdVelocity[j], 1,
950  stdVelocity[j], 1);
951  }
952  }
953  }
954  cnt += n_points;
955 
956  // Calculate total velocity in stdVelocity[0]
957  Vmath::Vmul( n_points, stdVelocity[0], 1, stdVelocity[0], 1,
958  stdVelocity[0], 1);
959  for(int k = 1; k < nvel; ++k)
960  {
961  Vmath::Vvtvp( n_points, stdVelocity[k], 1,
962  stdVelocity[k], 1,
963  stdVelocity[0], 1,
964  stdVelocity[0], 1);
965  }
966  pntVelocity = Vmath::Vmax( n_points, stdVelocity[0], 1);
967  maxV[el] = sqrt(pntVelocity);
968  }
969 
970  return maxV;
971  }
972 
973 
975  {
977  }
978 
979  /**
980  * At the start, the newest value is stored in array[nlevels-1]
981  * and the previous values in the first positions
982  * At the end, the extrapolated value is stored in array[nlevels-1]
983  * and the storage has been updated to included the new value
984  */
987  {
988  int nint = min(m_pressureCalls,m_intSteps);
989  int nlevels = array.num_elements();
990  int nPts = array[0].num_elements();
991 
992  // Update array
993  RollOver(array);
994 
995  // Extrapolate to outarray
996  Vmath::Smul(nPts, StifflyStable_Betaq_Coeffs[nint-1][nint-1],
997  array[nint-1], 1,
998  array[nlevels-1], 1);
999 
1000  for(int n = 0; n < nint-1; ++n)
1001  {
1002  Vmath::Svtvp(nPts, StifflyStable_Betaq_Coeffs[nint-1][n],
1003  array[n],1, array[nlevels-1],1,
1004  array[nlevels-1],1);
1005  }
1006  }
1007 
1008 
1009  /**
1010  * At the start, the newest value is stored in array[nlevels-1]
1011  * and the previous values in the first positions
1012  * At the end, the value of the bdf explicit part is stored in array[nlevels-1]
1013  * and the storage has been updated to included the new value
1014  */
1016  Array<OneD, Array<OneD, NekDouble> > &array)
1017  {
1018  int nint = min(m_pressureCalls,m_intSteps);
1019  int nlevels = array.num_elements();
1020  int nPts = array[0].num_elements();
1021 
1022  // Update array
1023  RollOver(array);
1024 
1025  // Extrapolate to outarray
1026  Vmath::Smul(nPts, StifflyStable_Alpha_Coeffs[nint-1][nint-1],
1027  array[nint-1], 1,
1028  array[nlevels-1], 1);
1029 
1030  for(int n = 0; n < nint-1; ++n)
1031  {
1032  Vmath::Svtvp(nPts, StifflyStable_Alpha_Coeffs[nint-1][n],
1033  array[n],1, array[nlevels-1],1,
1034  array[nlevels-1],1);
1035  }
1036  }
1037 
1038  /**
1039  * At the start, the newest value is stored in array[nlevels-1]
1040  * and the previous values in the first positions
1041  * At the end, the acceleration from BDF is stored in array[nlevels-1]
1042  * and the storage has been updated to included the new value
1043  */
1045  Array<OneD, Array<OneD, NekDouble> > &array)
1046  {
1047  int nlevels = array.num_elements();
1048  int nPts = array[0].num_elements();
1049 
1050 
1051  if(nPts)
1052  {
1053  // Update array
1054  RollOver(array);
1055 
1056  // Calculate acceleration using Backward Differentiation Formula
1057  Array<OneD, NekDouble> accelerationTerm (nPts, 0.0);
1058  if (m_pressureCalls > 2)
1059  {
1060  int acc_order = min(m_pressureCalls-2,m_intSteps);
1061  Vmath::Smul(nPts,
1062  StifflyStable_Gamma0_Coeffs[acc_order-1],
1063  array[0], 1,
1064  accelerationTerm, 1);
1065 
1066  for(int i = 0; i < acc_order; i++)
1067  {
1068  Vmath::Svtvp(nPts,
1069  -1*StifflyStable_Alpha_Coeffs[acc_order-1][i],
1070  array[i+1], 1,
1071  accelerationTerm, 1,
1072  accelerationTerm, 1);
1073  }
1074  }
1075  array[nlevels-1] = accelerationTerm;
1076  }
1077  }
1078 
1080  {
1081  int n, cnt;
1082  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
1083  {
1084  if((m_hbcType[n] == eHBCNeumann)||(m_hbcType[n] == eConvectiveOBC))
1085  {
1086  int nq = m_PBndExp[n]->GetNcoeffs();
1087  Vmath::Vcopy(nq, &(m_pressureHBCs[m_intSteps-1])[cnt], 1,
1088  &(m_PBndExp[n]->UpdateCoeffs()[0]), 1);
1089  cnt += nq;
1090  }
1091  }
1092  }
1093 }
void UpdateRobinPrimCoeff(void)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:222
void IProductNormVelocityOnHBC(const Array< OneD, const Array< OneD, NekDouble > > &Vel, Array< OneD, NekDouble > &IprodVn)
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:211
void ExtrapolateArray(Array< OneD, Array< OneD, NekDouble > > &array)
virtual void v_AddNormVelOnOBC(const int nbcoeffs, const int nreg, Array< OneD, Array< OneD, NekDouble > > &u)
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition: Extrapolate.h:226
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:242
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:50
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
std::string GetExpression(void) const
Definition: Equation.h:206
void AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > &array)
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
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:165
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
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
STL namespace.
HighOrderOutflowSharedPtr m_houtflow
Definition: Extrapolate.h:272
void CopyPressureHBCsToPbndExp(void)
LibUtilities::CommSharedPtr m_comm
Definition: Extrapolate.h:213
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
void AddPressureToOutflowBCs(NekDouble kinvis)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Velocity fields.
Definition: Extrapolate.h:219
static NekDouble StifflyStable_Gamma0_Coeffs[3]
Definition: Extrapolate.h:269
void EvaluateBDFArray(Array< OneD, Array< OneD, NekDouble > > &array)
static std::string def
Definition: Extrapolate.h:275
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:236
void MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
Definition: Extrapolate.h:392
void GenerateHOPBCMap(const LibUtilities::SessionReaderSharedPtr &pSsession)
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
Array< OneD, NekDouble > GetMaxStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
Array< OneD, HBCType > m_hbcType
Array of type of high order BCs for splitting shemes.
Definition: Extrapolate.h:216
virtual LibUtilities::TimeIntegrationMethod v_GetSubStepIntegrationMethod(void)
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:254
NekDouble Evaluate() const
Definition: Equation.h:102
LibUtilities::NekFactory< std::string, Extrapolate, const LibUtilities::SessionReaderSharedPtr &, Array< OneD, MultiRegions::ExpListSharedPtr > &, MultiRegions::ExpListSharedPtr &, const Array< OneD, int > &, const SolverUtils::AdvectionSharedPtr & > ExtrapolateFactory
Definition: Extrapolate.h:67
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
void Svtvm(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:518
double NekDouble
int m_pressureCalls
number of times the high-order pressure BCs have been called
Definition: Extrapolate.h:245
virtual void v_CalcNeumannPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
void IProductNormVelocityBCOnHBC(Array< OneD, NekDouble > &IprodVn)
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
virtual void v_CorrectPressureBCs(const Array< OneD, NekDouble > &pressure)
Array< OneD, Array< OneD, NekDouble > > m_pressureHBCs
Storage for current and previous levels of high order pressure boundary conditions.
Definition: Extrapolate.h:259
static NekDouble StifflyStable_Alpha_Coeffs[3][3]
Definition: Extrapolate.h:268
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:233
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:219
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
void CalcOutflowBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, NekDouble kinvis)
Array< OneD, Array< OneD, NekDouble > > m_iprodnormvel
Storage for current and previous levels of the inner product of normal velocity.
Definition: Extrapolate.h:262
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:239
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.
NekDouble m_timestep
Definition: Extrapolate.h:256
virtual ~Extrapolate()
Definition: Extrapolate.cpp:75
boost::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:222
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:299
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
Defines a callback function which evaluates the flux vector.
Definition: Advection.h:69
static NekDouble StifflyStable_Betaq_Coeffs[3][3]
Definition: Extrapolate.h:267
Provides a generic Factory class.
Definition: NekFactory.hpp:116
void AddNormVelOnOBC(const int nbcoeffs, const int nreg, Array< OneD, Array< OneD, NekDouble > > &u)
Definition: Extrapolate.h:423