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