Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
36 #include <IncNavierStokesSolver/EquationSystems/Extrapolate.h>
37 #include <LibUtilities/Communication/Comm.h>
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 
50  ExtrapolateFactory& GetExtrapolateFactory()
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(
60  const LibUtilities::SessionReaderSharedPtr pSession,
61  Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
62  MultiRegions::ExpListSharedPtr pPressure,
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 
75  Extrapolate::~Extrapolate()
76  {
77  }
78 
79  std::string Extrapolate::def =
80  LibUtilities::SessionReader::RegisterDefaultSolverInfo(
81  "StandardExtrapolate", "StandardExtrapolate");
82 
83  void Extrapolate::CalcExplicitDuDt(const Array<OneD, const Array<OneD, NekDouble> > &fields)
84  {
85  int nHBCs = m_acceleration[0].num_elements();
86 
87  // Adding extrapolated acceleration term to HOPBCs
88  Array<OneD, NekDouble> accelerationTerm(nHBCs, 0.0);
89 
90  // Rotate acceleration term
91  RollOver(m_acceleration);
92 
93  // update current normal of field on bc to m_acceleration[0];
94  IProductNormVelocityOnHBC(fields,m_acceleration[0]);
95 
96  //Calculate acceleration term at level n based on previous steps
97  if (m_pressureCalls > 2)
98  {
99  int acc_order = min(m_pressureCalls-2,m_intSteps);
100  Vmath::Smul(nHBCs, StifflyStable_Gamma0_Coeffs[acc_order-1],
101  m_acceleration[0], 1,
102  accelerationTerm, 1);
103 
104  for(int i = 0; i < acc_order; i++)
105  {
106  Vmath::Svtvp(nHBCs,
107  -1*StifflyStable_Alpha_Coeffs[acc_order-1][i],
108  m_acceleration[i+1], 1,
109  accelerationTerm, 1,
110  accelerationTerm, 1);
111  }
112  }
113 
114  Vmath::Svtvp(nHBCs, -1.0/m_timestep,
115  accelerationTerm, 1,
116  m_pressureHBCs[0], 1,
117  m_pressureHBCs[0], 1);
118  }
119 
120 
121  // Extrapolate to n+1
122  void Extrapolate::ExtrapolatePressureHBCs(void)
123  {
124  int nint = min(m_pressureCalls,m_intSteps);
125  int nlevels = m_pressureHBCs.num_elements();
126  int nHBCs = m_pressureHBCs[0].num_elements();
127  Vmath::Smul(nHBCs, StifflyStable_Betaq_Coeffs[nint-1][nint-1],
128  m_pressureHBCs[nint-1], 1,
129  m_pressureHBCs[nlevels-1], 1);
130 
131  for(int n = 0; n < nint-1; ++n)
132  {
133  Vmath::Svtvp(nHBCs, StifflyStable_Betaq_Coeffs[nint-1][n],
134  m_pressureHBCs[n],1,m_pressureHBCs[nlevels-1],1,
135  m_pressureHBCs[nlevels-1],1);
136  }
137  }
138 
139  // Copy values of [dP/dn]^{n+1} in the pressure bcs storage.
140  // m_pressureHBCS[nlevels-1] will be cancelled at next time step
141  void Extrapolate::CopyPressureHBCsToPbndExp(void)
142  {
143  int n,cnt;
144  int nlevels = m_pressureHBCs.num_elements();
145 
146  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
147  {
148  // High order boundary condition;
149  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"H"))
150  {
151  int nq = m_PBndExp[n]->GetNcoeffs();
152  Vmath::Vcopy(nq, &(m_pressureHBCs[nlevels-1])[cnt], 1,
153  &(m_PBndExp[n]->UpdateCoeffs()[0]), 1);
154  cnt += nq;
155  }
156 
157  }
158  }
159 
160  /**
161  * Unified routine for calculation high-oder terms
162  */
163  void Extrapolate::v_CalcNeumannPressureBCs(
164  const Array<OneD, const Array<OneD, NekDouble> > &fields,
165  const Array<OneD, const Array<OneD, NekDouble> > &N,
166  NekDouble kinvis)
167  {
168  Array<OneD, NekDouble> Pvals;
169  StdRegions::StdExpansionSharedPtr Pbc;
170  StdRegions::StdExpansionSharedPtr elmt;
171 
172  int pindex=N.num_elements();
173 
174  Array<OneD, Array<OneD, const NekDouble> > Velocity(m_curl_dim);
175  Array<OneD, Array<OneD, const NekDouble> > Advection(m_bnd_dim);
176 
177  Array<OneD, Array<OneD, NekDouble> > BndValues(m_bnd_dim);
178  Array<OneD, Array<OneD, NekDouble> > Q(m_bnd_dim);
179 
180  for(int i = 0; i < m_bnd_dim; i++)
181  {
182  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,0.0);
183  Q[i] = Array<OneD, NekDouble> (m_pressureBCsElmtMaxPts,0.0);
184  }
185 
186  for(int j = 0 ; j < m_HBCdata.num_elements() ; j++)
187  {
188  /// Casting the boundary expansion to the specific case
189  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion>
190  (m_PBndExp[m_HBCdata[j].m_bndryID]->GetExp(m_HBCdata[j].m_bndElmtID));
191 
192  /// Picking up the element where the HOPBc is located
193  elmt = m_pressure->GetExp(m_HBCdata[j].m_globalElmtID);
194 
195  /// Assigning
196  for(int i = 0; i < m_bnd_dim; i++)
197  {
198  Velocity[i] = fields[i] + m_HBCdata[j].m_physOffset;
199  Advection[i] = N[i] + m_HBCdata[j].m_physOffset;
200  }
201 
202  // for the 3DH1D case we need to grab the conjugate mode
203  if(m_pressure->GetExpType() == MultiRegions::e3DH1D)
204  {
205  Velocity[2] = fields[2] + m_HBCdata[j].m_assPhysOffset;
206  }
207 
208  /// Calculating the curl-curl and storing it in Q
209  CurlCurl(Velocity,Q,j);
210 
211  // Mounting advection component into the high-order condition
212  for(int i = 0; i < m_bnd_dim; i++)
213  {
214  MountHOPBCs(m_HBCdata[j].m_ptsInElmt,kinvis,Q[i],Advection[i]);
215  }
216 
217  // put in m_pressureBCs[0]
218  Pvals = m_pressureHBCs[0] + m_HBCdata[j].m_coeffOffset;
219 
220  // Getting values on the edge and filling the pressure boundary expansion
221  // and the acceleration term. Multiplication by the normal is required
222  switch(m_fields[pindex]->GetExpType())
223  {
224  case MultiRegions::e2D:
225  case MultiRegions::e3DH1D:
226  {
227  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID,Pbc,Q[0],BndValues[0]);
228  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID,Pbc,Q[1],BndValues[1]);
229  Pbc->NormVectorIProductWRTBase(BndValues[0],BndValues[1],Pvals);
230  }
231  break;
232  case MultiRegions::e3DH2D:
233  {
234  if(m_HBCdata[j].m_elmtTraceID == 0)
235  {
236  Pvals[0] = -1.0*Q[0][0];
237  }
238  else if (m_HBCdata[j].m_elmtTraceID == 1)
239  {
240  Pvals[0] = Q[0][m_HBCdata[j].m_ptsInElmt-1];
241  }
242  else
243  {
244  ASSERTL0(false,
245  "In the 3D homogeneous 2D approach BCs edge "
246  "ID can be just 0 or 1 ");
247  }
248  }
249  break;
250  case MultiRegions::e3D:
251  {
252  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID,Pbc,Q[0],BndValues[0]);
253  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID,Pbc,Q[1],BndValues[1]);
254  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID,Pbc,Q[2],BndValues[2]);
255  Pbc->NormVectorIProductWRTBase(BndValues[0],BndValues[1],BndValues[2],Pvals);
256  }
257  break;
258  default:
259  ASSERTL0(0,"Dimension not supported");
260  break;
261  }
262  }
263  }
264 
265  // do nothing unless otherwise defined.
266  void Extrapolate::v_CorrectPressureBCs( const Array<OneD, NekDouble> &pressure)
267  {
268  }
269 
270  void Extrapolate::CalcOutflowBCs(
271  const Array<OneD, const Array<OneD, NekDouble> > &fields,
272  NekDouble kinvis)
273  {
274  static bool init = true;
275  static bool noHOBC = false;
276 
277  if(noHOBC == true)
278  {
279  return;
280  }
281 
282  if(init) // set up storage for boundary velocity at outflow
283  {
284  init = false;
285  int totbndpts = 0;
286  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
287  {
288  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"HOutflow"))
289  {
290  totbndpts += m_PBndExp[n]->GetTotPoints();
291  }
292  }
293 
294  if(totbndpts == 0)
295  {
296  noHOBC = true;
297  return;
298  }
299 
300  m_outflowVel = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (m_curl_dim);
301  for(int i = 0; i < m_curl_dim; ++i)
302  {
303  m_outflowVel[i] = Array<OneD, Array<OneD, NekDouble> >(m_curl_dim);
304  for(int j = 0; j < m_curl_dim; ++j)
305  {
306  // currently just set up for 2nd order extrapolation
307  m_outflowVel[i][j] = Array<OneD, NekDouble>(totbndpts,0.0);
308  }
309  }
310 
311  if (m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
312  {
313  m_PhyoutfVel = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (m_curl_dim);
314 
315  for(int i = 0; i < m_curl_dim; ++i)
316  {
317  m_PhyoutfVel[i] = Array<OneD, Array<OneD, NekDouble> > (m_curl_dim);
318  for(int j = 0; j < m_curl_dim; ++j)
319  {
320  // currently just set up for 2nd order extrapolation
321  m_PhyoutfVel[i][j] = Array<OneD, NekDouble> (totbndpts,0.0);
322  }
323  }
324 
325  m_nonlinearterm_phys = Array<OneD, NekDouble> (totbndpts,0.0);
326  m_nonlinearterm_coeffs = Array<OneD, NekDouble> (totbndpts,0.0);
327 
328  m_PBndCoeffs = Array<OneD, NekDouble> (totbndpts,0.0);
329  m_UBndCoeffs = Array<OneD, Array<OneD, NekDouble> > (m_curl_dim);
330  for(int i = 0; i < m_curl_dim; ++i)
331  {
332  m_UBndCoeffs[i] = Array<OneD, NekDouble> (totbndpts);
333  }
334  Array<OneD, unsigned int> planes;
335  planes = m_pressure->GetZIDs();
336  int num_planes = planes.num_elements();
337  m_expsize_per_plane = Array<OneD, unsigned int> (m_PBndConds.num_elements());
338  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
339  {
340  int exp_size = m_PBndExp[n]->GetExpSize();
341  m_expsize_per_plane[n] = exp_size/num_planes;
342  }
343  m_totexps_per_plane = 0;
344  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
345  {
346  m_totexps_per_plane += m_PBndExp[n]->GetExpSize()/num_planes;
347  }
348  }
349  }
350 
351  StdRegions::StdExpansionSharedPtr Bc,Pbc;
352  Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr> >
353  UBndConds(m_curl_dim);
354  Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
355  UBndExp(m_curl_dim);
356 
357  for (int i = 0; i < m_curl_dim; ++i)
358  {
359  UBndConds[i] = m_fields[m_velocity[i]]->GetBndConditions();
360  UBndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
361  }
362 
363  Array<OneD, Array<OneD, NekDouble> > BndValues(m_curl_dim);
364  Array<OneD, Array<OneD, NekDouble> > BndElmt (m_curl_dim);
365  Array<OneD, Array<OneD, NekDouble> > nGradu (m_curl_dim);
366  Array<OneD, NekDouble > gradtmp (m_pressureBCsElmtMaxPts),
367  fgradtmp(m_pressureBCsElmtMaxPts);
368 
369  nGradu[0] = Array<OneD, NekDouble>(m_curl_dim*m_pressureBCsMaxPts);
370  for(int i = 0; i < m_curl_dim; ++i)
371  {
372  BndElmt[i] = Array<OneD, NekDouble> (m_pressureBCsElmtMaxPts,
373  0.0);
374  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,0.0);
375  nGradu[i] = nGradu[0] + i*m_pressureBCsMaxPts;
376  RollOver(m_outflowVel[i]);
377 
378  if (m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
379  {
380  RollOver(m_PhyoutfVel[i]);
381  }
382  }
383 
384  int nbc,cnt,cnt_start;
385  int veloffset = 0;
386  int nint = min(m_pressureCalls,m_intSteps);
387 
388  StdRegions::StdExpansionSharedPtr elmt;
389  Array<OneD, NekDouble> PBCvals, UBCvals;
390  Array<OneD, Array<OneD, NekDouble> > ubc(m_curl_dim);
391  Array<OneD, Array<OneD, NekDouble> > normals;
392 
393  cnt = 0;
394  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
395  {
396  // Do outflow boundary conditions if they exist
397  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"HOutflow"))
398  {
399  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i,cnt++)
400  {
401  cnt = max(cnt,m_PBndExp[n]->GetTotPoints());
402  }
403  }
404  }
405 
406  for(int i =0; i < m_curl_dim; ++i)
407  {
408  ubc[i] = Array<OneD, NekDouble>(cnt);
409  }
410 
411  NekDouble U0,delta;
412  m_session->LoadParameter("U0_HighOrderBC",U0,1.0);
413  m_session->LoadParameter("Delta_HighOrderBC",delta,1/20.0);
414 
415  cnt = 0;
416  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
417  {
418  cnt_start = cnt;
419 
420  // Do outflow boundary conditions if they exist
421  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"HOutflow"))
422  {
423 
424  if (m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
425  {
426  int veloffset = 0;
427  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i, cnt++)
428  {
429  // find element and edge of this expansion.
430  Bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion>
431  (m_PBndExp[n]->GetExp(i));
432 
433  int elmtid = m_pressureBCtoElmtID[cnt];
434  elmt = m_fields[0]->GetExp(elmtid);
435  int offset = m_fields[0]->GetPhys_Offset(elmtid);
436 
437  int boundary = m_pressureBCtoTraceID[cnt];
438 
439  // Determine extrapolated U,V values
440  int nq = elmt->GetTotPoints();
441  int nbc = m_PBndExp[n]->GetExp(i)->GetTotPoints();
442  // currently just using first order approximation here.
443  // previously have obtained value from m_integrationSoln
444  Array<OneD, NekDouble> veltmp;
445 
446  for(int j = 0; j < m_curl_dim; ++j)
447  {
448  Vmath::Vcopy(nq, &fields[m_velocity[j]][offset], 1,
449  &BndElmt[j][0], 1);
450  elmt->GetTracePhysVals(boundary,Bc,BndElmt[j],
451  veltmp = m_outflowVel[j][0] + veloffset);
452  }
453  veloffset += nbc;
454  }
455 
456  // for velocity on the outflow boundary in e3DH1D,
457  // we need to make a backward fourier transformation
458  // to get the physical coeffs at the outflow BCs.
459  for(int j = 0; j < m_curl_dim; ++j)
460  {
461  m_PBndExp[n]->HomogeneousBwdTrans(
462  m_outflowVel[j][0],
463  m_PhyoutfVel[j][0]);
464  }
465 
466  cnt = cnt_start;
467  veloffset = 0;
468  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i, cnt++)
469  {
470 
471  int elmtid = m_pressureBCtoElmtID[cnt];
472  elmt = m_fields[0]->GetExp(elmtid);
473  int nbc = m_PBndExp[n]->GetExp(i)->GetTotPoints();
474 
475  Array<OneD, NekDouble> veltmp(nbc,0.0),
476  normDotu(nbc,0.0), utot(nbc,0.0);
477  int boundary = m_pressureBCtoTraceID[cnt];
478  normals=elmt->GetSurfaceNormal(boundary);
479 
480  // extrapolate velocity
481  if(nint <= 1)
482  {
483  for(int j = 0; j < m_curl_dim; ++j)
484  {
485  Vmath::Vcopy(nbc,
486  veltmp = m_PhyoutfVel[j][0] +veloffset, 1,
487  BndValues[j], 1);
488  }
489  }
490  else // only set up for 2nd order extrapolation
491  {
492  for(int j = 0; j < m_curl_dim; ++j)
493  {
494  Vmath::Smul(nbc, 2.0,
495  veltmp = m_PhyoutfVel[j][0] + veloffset, 1,
496  BndValues[j], 1);
497  Vmath::Svtvp(nbc, -1.0,
498  veltmp = m_PhyoutfVel[j][1] + veloffset, 1,
499  BndValues[j], 1,
500  BndValues[j], 1);
501  }
502  }
503 
504  // Set up |u|^2, n.u in physical space
505  for(int j = 0; j < m_curl_dim; ++j)
506  {
507  Vmath::Vvtvp(nbc, BndValues[j], 1, BndValues[j], 1,
508  utot, 1, utot, 1);
509  }
510  for(int j = 0; j < m_bnd_dim; ++j)
511  {
512  Vmath::Vvtvp(nbc, normals[j], 1, BndValues[j], 1,
513  normDotu, 1, normDotu, 1);
514  }
515 
516  int Offset = m_PBndExp[n]->GetPhys_Offset(i);
517 
518  for(int k = 0; k < nbc; ++k)
519  {
520  // calculate the nonlinear term (kinetic energy
521  // multiplies step function) in physical space
522  NekDouble fac = 0.5*(1.0-tanh(normDotu[k]/(U0*delta)));
523  m_nonlinearterm_phys[k + Offset] = 0.5 * utot[k] * fac;
524  }
525 
526  veloffset += nbc;
527  }
528 
529  // for e3DH1D, we need to make a forward fourier transformation
530  // for the kinetic energy term (nonlinear)
531  UBndExp[0][n]->HomogeneousFwdTrans(
532  m_nonlinearterm_phys,
533  m_nonlinearterm_coeffs);
534 
535  // for e3DH1D, we need to make a forward fourier transformation
536  // for Dirichlet pressure boundary condition that is from input file
537  m_PBndExp[n]->HomogeneousFwdTrans(
538  m_PBndExp[n]->UpdatePhys(),
539  m_PBndCoeffs);
540  // for e3DH1D, we need to make a forward fourier transformation
541  // for Neumann velocity boundary condition that is from input file
542  for (int j = 0; j < m_curl_dim; ++j)
543  {
544  UBndExp[j][n]->HomogeneousFwdTrans(
545  UBndExp[j][n]->UpdatePhys(),
546  m_UBndCoeffs[j]);
547  }
548  }
549 
550  cnt = cnt_start;
551  veloffset = 0;
552  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i,cnt++)
553  {
554  // find element and edge of this expansion.
555  Bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion>
556  (m_PBndExp[n]->GetExp(i));
557 
558  int elmtid = m_pressureBCtoElmtID[cnt];
559  elmt = m_fields[0]->GetExp(elmtid);
560  int offset = m_fields[0]->GetPhys_Offset(elmtid);
561 
562  // Determine extrapolated U,V values
563  int nq = elmt->GetTotPoints();
564 
565  // currently just using first order approximation here.
566  // previously have obtained value from m_integrationSoln
567  for(int j = 0; j < m_bnd_dim; ++j)
568  {
569  Vmath::Vcopy(nq, &fields[m_velocity[j]][offset], 1,
570  &BndElmt[j][0], 1);
571  }
572 
573  int nbc = m_PBndExp[n]->GetExp(i)->GetTotPoints();
574  int boundary = m_pressureBCtoTraceID[cnt];
575 
576  Array<OneD, NekDouble> ptmp(nbc,0.0),
577  divU(nbc,0.0);
578 
579  normals=elmt->GetSurfaceNormal(boundary);
580  Vmath::Zero(m_bnd_dim*m_pressureBCsMaxPts,nGradu[0],1);
581 
582  for (int j = 0; j < m_bnd_dim; j++)
583  {
584  // Calculate Grad u = du/dx, du/dy, du/dz, etc.
585  for (int k = 0; k< m_bnd_dim; k++)
586  {
587  elmt->PhysDeriv(MultiRegions::DirCartesianMap[k],
588  BndElmt[j], gradtmp);
589  elmt->GetTracePhysVals(boundary, Bc, gradtmp,
590  fgradtmp);
591  Vmath::Vvtvp(nbc,normals[k], 1, fgradtmp, 1,
592  nGradu[j], 1, nGradu[j],1);
593  if(j == k)
594  {
595  Vmath::Vadd(nbc,fgradtmp, 1, divU, 1, divU, 1);
596  }
597  }
598  }
599 
600  if (m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
601  {
602  // Set up |u|^2, n.u, div(u), and (n.grad(u) . n) for
603  // pressure condition
604  for(int j = 0; j < m_bnd_dim; ++j)
605  {
606  Vmath::Vvtvp(nbc, normals[j], 1, nGradu[j], 1,
607  ptmp, 1, ptmp, 1);
608  }
609  int p_offset = m_PBndExp[n]->GetPhys_Offset(i);
610 
611  for(int k = 0; k < nbc; ++k)
612  {
613  // Set up Dirichlet pressure condition and
614  // store in ptmp (m_UBndCoeffs contains Fourier Coeffs of the
615  // function from the input file )
616 
617  ptmp[k] = kinvis * ptmp[k]
618  - m_nonlinearterm_coeffs[k + p_offset]
619  - m_PBndCoeffs[k + p_offset];
620  }
621 
622  int u_offset = UBndExp[0][n]->GetPhys_Offset(i);
623 
624  for(int j = 0; j < m_bnd_dim; ++j)
625  {
626  for(int k = 0; k < nbc; ++k)
627  {
628  ubc[j][k + u_offset] = (1.0 / kinvis)
629  * (m_UBndCoeffs[j][k + u_offset]
630  + m_nonlinearterm_coeffs[k + u_offset]
631  * normals[j][k]);
632  }
633  }
634 
635  // boundary condition for velocity in homogenous direction
636  for(int k = 0; k < nbc; ++k)
637  {
638  ubc[m_bnd_dim][k + u_offset] = (1.0 / kinvis)
639  * m_UBndCoeffs[m_bnd_dim][k + u_offset];
640  }
641 
642  u_offset = UBndExp[m_bnd_dim][n]->GetPhys_Offset(i);
643  UBCvals = UBndExp[m_bnd_dim][n]->UpdateCoeffs()
644  + UBndExp[m_bnd_dim][n]->GetCoeff_Offset(i);
645  Bc->IProductWRTBase(ubc[m_bnd_dim] + u_offset, UBCvals);
646  }
647  else
648  {
649 
650  Array<OneD, NekDouble> veltmp, utot(nbc,0.0),
651  normDotu(nbc,0.0);
652  // extract velocity and store
653  for(int j = 0; j < m_bnd_dim; ++j)
654  {
655  elmt->GetTracePhysVals(boundary,Bc,BndElmt[j],
656  veltmp = m_outflowVel[j][0] + veloffset);
657  }
658 
659  // extrapolate velocity
660  if(nint <= 1)
661  {
662  for(int j = 0; j < m_bnd_dim; ++j)
663  {
664  Vmath::Vcopy(nbc,
665  veltmp = m_outflowVel[j][0]
666  +veloffset, 1,
667  BndValues[j],1);
668  }
669  }
670  else // only set up for 2nd order extrapolation
671  {
672  for(int j = 0; j < m_bnd_dim; ++j)
673  {
674  Vmath::Smul(nbc, 2.0,
675  veltmp = m_outflowVel[j][0]
676  + veloffset, 1,
677  BndValues[j], 1);
678  Vmath::Svtvp(nbc, -1.0,
679  veltmp = m_outflowVel[j][1]
680  + veloffset, 1,
681  BndValues[j], 1,
682  BndValues[j], 1);
683  }
684  }
685 
686  // Set up |u|^2, n.u, div(u), and (n.grad(u) . n) for
687  // pressure condition
688  for(int j = 0; j < m_bnd_dim; ++j)
689  {
690  Vmath::Vvtvp(nbc, BndValues[j], 1, BndValues[j], 1,
691  utot, 1, utot, 1);
692  Vmath::Vvtvp(nbc, normals[j], 1, BndValues[j], 1,
693  normDotu, 1, normDotu, 1);
694  Vmath::Vvtvp(nbc, normals[j], 1, nGradu[j], 1,
695  ptmp, 1, ptmp, 1);
696  }
697 
698  PBCvals = m_PBndExp[n]->GetPhys() +
699  m_PBndExp[n]->GetPhys_Offset(i);
700 
701  for(int k = 0; k < nbc; ++k)
702  {
703  NekDouble fac = 0.5*(1.0-tanh(normDotu[k]/(U0*delta)));
704 
705  // Set up Dirichlet pressure condition and
706  // store in ptmp (PBCvals contains a
707  // function from the input file )
708  ptmp[k] = kinvis * ptmp[k] - 0.5 * utot[k] * fac
709  + PBCvals[k];
710  }
711 
712  int u_offset = UBndExp[0][n]->GetPhys_Offset(i);
713 
714  for(int j = 0; j < m_bnd_dim; ++j)
715  {
716  UBCvals = UBndExp[j][n]->GetPhys()
717  + UBndExp[j][n]->GetPhys_Offset(i);
718 
719  for(int k = 0; k < nbc; ++k)
720  {
721  NekDouble fac = 0.5 * (1.0 - tanh(normDotu[k]
722  / (U0 * delta)));
723  ubc[j][k + u_offset] = (1.0 / kinvis)
724  * (UBCvals[k] + 0.5 * utot[k] * fac
725  * normals[j][k]);
726  }
727  }
728  }
729 
730  // set up pressure boundary condition
731  PBCvals = m_PBndExp[n]->UpdateCoeffs()
732  + m_PBndExp[n]->GetCoeff_Offset(i);
733  m_PBndExp[n]->GetExp(i)->FwdTrans(ptmp,PBCvals);
734 
735  veloffset += nbc;
736  }
737 
738  // Now set up Velocity conditions.
739  for(int j = 0; j < m_bnd_dim; j++)
740  {
741  if(boost::iequals(UBndConds[j][n]->GetUserDefined(),"HOutflow"))
742  {
743  cnt = cnt_start;
744  for(int i = 0; i < UBndExp[0][n]->GetExpSize();
745  ++i, cnt++)
746  {
747  Pbc = StdRegions::StdExpansionSharedPtr
748  (m_PBndExp[n]->GetExp(i));
749  Bc = StdRegions::StdExpansionSharedPtr
750  (UBndExp[0][n]->GetExp(i));
751 
752  nbc = UBndExp[0][n]->GetExp(i)->GetTotPoints();
753  int boundary = m_pressureBCtoTraceID[cnt];
754 
755  Array<OneD, NekDouble> pb(nbc), ub(nbc);
756  int elmtid = m_pressureBCtoElmtID[cnt];
757 
758  elmt = m_fields[0]->GetExp(elmtid);
759 
760  normals = elmt->GetSurfaceNormal(boundary);
761 
762  // Get p from projected boundary condition
763  PBCvals = m_PBndExp[n]->UpdateCoeffs()
764  + m_PBndExp[n]->GetCoeff_Offset(i);
765  Pbc->BwdTrans(PBCvals,pb);
766 
767  int u_offset = UBndExp[j][n]->GetPhys_Offset(i);
768 
769  for(int k = 0; k < nbc; ++k)
770  {
771  ub[k] = ubc[j][k + u_offset]
772  + pb[k] * normals[j][k] / kinvis;
773  }
774 
775  UBCvals = UBndExp[j][n]->UpdateCoeffs()
776  + UBndExp[j][n]->GetCoeff_Offset(i);
777  Bc->IProductWRTBase(ub,UBCvals);
778  }
779  }
780  }
781  }
782  else
783  {
784  cnt += m_PBndExp[n]->GetExpSize();
785  }
786  }
787  }
788 
789  /**
790  * Curl Curl routine - dimension dependent
791  */
792  void Extrapolate::CurlCurl(
793  Array<OneD, Array<OneD, const NekDouble> > &Vel,
794  Array<OneD, Array<OneD, NekDouble> > &Q,
795  const int j)
796  {
797  StdRegions::StdExpansionSharedPtr elmt
798  = m_fields[0]->GetExp(m_HBCdata[j].m_globalElmtID);
799 
800  Array<OneD,NekDouble> Vx(m_pressureBCsElmtMaxPts);
801  Array<OneD,NekDouble> Uy(m_pressureBCsElmtMaxPts);
802 
803  switch(m_fields[0]->GetExpType())
804  {
805  case MultiRegions::e2D:
806  {
807  Array<OneD,NekDouble> Dummy(m_pressureBCsElmtMaxPts);
808 
809  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[1], Vx);
810  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Vel[0], Uy);
811 
812  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Uy, 1, Dummy, 1);
813 
814  elmt->PhysDeriv(Dummy,Q[1],Q[0]);
815 
816  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, -1.0, Q[1], 1, Q[1], 1);
817  }
818  break;
819 
820  case MultiRegions::e3DH1D:
821  {
822  Array<OneD,NekDouble> Wz(m_pressureBCsElmtMaxPts);
823 
824  Array<OneD,NekDouble> Dummy1(m_pressureBCsElmtMaxPts);
825  Array<OneD,NekDouble> Dummy2(m_pressureBCsElmtMaxPts);
826 
827  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[1], Vx);
828  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Vel[0], Uy);
829  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_wavenumber[j],
830  Vel[2], 1, Wz, 1);
831 
832  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Vx, Dummy1);
833  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Uy, Dummy2);
834  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Dummy1, 1, Dummy2, 1,
835  Q[0], 1);
836  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
837  Vel[0], 1, Dummy1, 1);
838  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Wz, Dummy2);
839  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Q[0], 1, Dummy1, 1,
840  Q[0], 1);
841  Vmath::Vadd(m_HBCdata[j].m_ptsInElmt, Q[0], 1, Dummy2, 1,
842  Q[0], 1);
843 
844  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Wz, Dummy1);
845  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
846  Vel[1], 1, Dummy2, 1);
847  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Dummy1, 1, Dummy2, 1,
848  Q[1], 1);
849  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vx, Dummy1);
850  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Uy, Dummy2);
851  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Q[1], 1, Dummy1, 1,
852  Q[1], 1);
853  Vmath::Vadd(m_HBCdata[j].m_ptsInElmt, Q[1], 1, Dummy2, 1,
854  Q[1], 1);
855  }
856  break;
857  case MultiRegions::e3DH2D:
858  {
859  Array<OneD,NekDouble> Wx(m_pressureBCsElmtMaxPts);
860  Array<OneD,NekDouble> Wz(m_pressureBCsElmtMaxPts);
861  Array<OneD,NekDouble> Uz(m_pressureBCsElmtMaxPts);
862  Array<OneD,NekDouble> qz(m_pressureBCsElmtMaxPts);
863  Array<OneD,NekDouble> qy(m_pressureBCsElmtMaxPts);
864 
865  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[2], Wx);
866  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[1], Vx);
867 
868  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
869  Vel[0], 1, Uy, 1);
870 
871  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_wavenumber[j],
872  Vel[0], 1, Uz, 1);
873 
874  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Wz, 1, Wx, 1,
875  qy, 1);
876  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Uy, 1,
877  qz, 1);
878 
879  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
880  qz, 1, Uy, 1);
881 
882  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_wavenumber[j],
883  qy, 1, Uz, 1);
884 
885  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Uy, 1, Uz, 1,
886  Q[0], 1);
887  }
888  break;
889  case MultiRegions::e3D:
890  {
891  Array<OneD,NekDouble> Dummy(m_pressureBCsElmtMaxPts);
892  Array<OneD,NekDouble> Vz(m_pressureBCsElmtMaxPts);
893  Array<OneD,NekDouble> Uz(m_pressureBCsElmtMaxPts);
894  Array<OneD,NekDouble> Wx(m_pressureBCsElmtMaxPts);
895  Array<OneD,NekDouble> Wy(m_pressureBCsElmtMaxPts);
896 
897  elmt->PhysDeriv(Vel[0], Dummy, Uy, Uz);
898  elmt->PhysDeriv(Vel[1], Vx, Dummy, Vz);
899  elmt->PhysDeriv(Vel[2], Wx, Wy, Dummy);
900 
901  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Wy, 1, Vz, 1, Q[0], 1);
902  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Uz, 1, Wx, 1, Q[1], 1);
903  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Uy, 1, Q[2], 1);
904 
905  elmt->PhysDeriv(Q[0], Dummy, Wy, Vx);
906  elmt->PhysDeriv(Q[1], Wx, Dummy, Uz);
907  elmt->PhysDeriv(Q[2], Vz, Uy, Dummy);
908 
909  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Uy, 1, Uz, 1, Q[0], 1);
910  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Vz, 1, Q[1], 1);
911  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Wx, 1, Wy, 1, Q[2], 1);
912  }
913  break;
914  default:
915  ASSERTL0(0,"Dimension not supported");
916  break;
917  }
918  }
919 
920 
921  void Extrapolate::IProductNormVelocityOnHBC(
922  const Array<OneD, const Array<OneD, NekDouble> > &Vel,
923  Array<OneD, NekDouble> &IProdVn)
924  {
925 
926  int i,j,n;
927 
928  int pindex=m_fields.num_elements()-1;
929  StdRegions::StdExpansionSharedPtr Pbc,elmt;
930 
931  Array<OneD, Array<OneD, NekDouble> > velbc(m_bnd_dim);
932  velbc[0] = Array<OneD, NekDouble>(m_bnd_dim*m_pressureBCsMaxPts);
933  for(i = 1; i < m_bnd_dim; ++i)
934  {
935  velbc[i] = velbc[i-1] + m_pressureBCsMaxPts;
936  }
937 
938  Array<OneD, NekDouble> Velmt,IProdVnTmp;
939 
940  for(n = 0; n < m_HBCdata.num_elements(); ++n)
941  {
942  Pbc = m_PBndExp[m_HBCdata[n].m_bndryID]->GetExp(m_HBCdata[n].m_bndElmtID);
943 
944  /// Picking up the element where the HOPBc is located
945  elmt = m_fields[pindex]->GetExp(m_HBCdata[n].m_globalElmtID);
946 
947  // Get velocity bc
948  for(j = 0; j < m_bnd_dim; ++j)
949  {
950  // Get edge values and put into velbc
951  Velmt = Vel[j] + m_HBCdata[n].m_physOffset;
952  elmt->GetTracePhysVals(m_HBCdata[n].m_elmtTraceID,Pbc,Velmt,velbc[j]);
953  }
954 
955  IProdVnTmp = IProdVn + m_HBCdata[n].m_coeffOffset;
956 
957  // Evaluate Iproduct wrt norm term
958  Pbc->NormVectorIProductWRTBase(velbc,IProdVnTmp);
959  }
960  }
961 
962  void Extrapolate::IProductNormVelocityBCOnHBC(Array<OneD, NekDouble> &IProdVn)
963  {
964 
965  int i,j,n;
966 
967  StdRegions::StdExpansionSharedPtr Pbc;
968 
969  Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr> > VelBndExp(m_bnd_dim);
970 
971  for(i = 0; i < m_bnd_dim; ++i)
972  {
973  VelBndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
974  }
975 
976  Array<OneD, Array<OneD, NekDouble> > velbc(m_bnd_dim);
977  velbc[0] = Array<OneD, NekDouble>(m_bnd_dim*m_pressureBCsMaxPts);
978  for(i = 1; i < m_bnd_dim; ++i)
979  {
980  velbc[i] = velbc[i-1] + m_pressureBCsMaxPts;
981  }
982 
983  Array<OneD, NekDouble> Velmt,IProdVnTmp;
984  int bndid,elmtid;
985 
986  for(n = 0; n < m_HBCdata.num_elements(); ++n)
987  {
988  bndid = m_HBCdata[n].m_bndryID;
989  elmtid = m_HBCdata[n].m_bndElmtID;
990 
991  // Get velocity bc
992  for(j = 0; j < m_bnd_dim; ++j)
993  {
994  VelBndExp[j][bndid]->GetExp(elmtid)->BwdTrans(VelBndExp[j][bndid]->GetCoeffs() +
995  VelBndExp[j][bndid]->GetCoeff_Offset(elmtid),
996  velbc[j]);
997  }
998 
999  IProdVnTmp = IProdVn + m_HBCdata[n].m_coeffOffset;
1000  // Evaluate Iproduct wrt norm term
1001  Pbc = m_PBndExp[m_HBCdata[n].m_bndryID]->GetExp(m_HBCdata[n].m_bndElmtID);
1002  Pbc->NormVectorIProductWRTBase(velbc,IProdVnTmp);
1003  }
1004  }
1005 
1006  /**
1007  * Function to roll time-level storages to the next step layout.
1008  * The stored data associated with the oldest time-level
1009  * (not required anymore) are moved to the top, where they will
1010  * be overwritten as the solution process progresses.
1011  */
1012  void Extrapolate::RollOver(Array<OneD, Array<OneD, NekDouble> > &input)
1013  {
1014  int nlevels = input.num_elements();
1015 
1016  Array<OneD, NekDouble> tmp;
1017 
1018  tmp = input[nlevels-1];
1019 
1020  for(int n = nlevels-1; n > 0; --n)
1021  {
1022  input[n] = input[n-1];
1023  }
1024 
1025  input[0] = tmp;
1026  }
1027 
1028 
1029  /**
1030  * Map to directly locate HOPBCs position and offsets in all scenarios
1031  */
1032  void Extrapolate::GenerateHOPBCMap()
1033  {
1034 
1035  int pindex=m_fields.num_elements()-1;
1036 
1037  m_PBndConds = m_pressure->GetBndConditions();
1038  m_PBndExp = m_pressure->GetBndCondExpansions();
1039 
1040  // Set up mapping from pressure boundary condition to pressure element
1041  // details.
1042  m_pressure->GetBoundaryToElmtMap(m_pressureBCtoElmtID,
1043  m_pressureBCtoTraceID);
1044 
1045  // find the maximum values of points for pressure BC evaluation
1046  m_pressureBCsMaxPts = 0;
1047  m_pressureBCsElmtMaxPts = 0;
1048  int cnt, n;
1049  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
1050  {
1051  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i)
1052  {
1053  m_pressureBCsMaxPts = max(m_pressureBCsMaxPts,
1054  m_PBndExp[n]->GetExp(i)->GetTotPoints());
1055  m_pressureBCsElmtMaxPts = max(m_pressureBCsElmtMaxPts,
1056  m_pressure->GetExp(m_pressureBCtoElmtID[cnt++])
1057  ->GetTotPoints());
1058  }
1059  }
1060 
1061  // Storage array for high order pressure BCs
1062  m_pressureHBCs = Array<OneD, Array<OneD, NekDouble> > (m_intSteps);
1063  m_acceleration = Array<OneD, Array<OneD, NekDouble> > (m_intSteps + 1);
1064 
1065  int HBCnumber = 0;
1066  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
1067  {
1068  // High order boundary condition;
1069  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"H"))
1070  {
1071  cnt += m_PBndExp[n]->GetNcoeffs();
1072  HBCnumber += m_PBndExp[n]->GetExpSize();
1073  }
1074  }
1075 
1076  int checkHBC = HBCnumber;
1077  m_comm->AllReduce(checkHBC,LibUtilities::ReduceSum);
1078  //ASSERTL0(checkHBC > 0 ,"At least one high-order pressure boundary "
1079  //"condition is required for scheme "
1080  //"consistency");
1081 
1082  m_acceleration[0] = Array<OneD, NekDouble>(cnt, 0.0);
1083  for(n = 0; n < m_intSteps; ++n)
1084  {
1085  m_pressureHBCs[n] = Array<OneD, NekDouble>(cnt, 0.0);
1086  m_acceleration[n+1] = Array<OneD, NekDouble>(cnt, 0.0);
1087  }
1088 
1089  m_pressureCalls = 0;
1090 
1091  switch(m_pressure->GetExpType())
1092  {
1093  case MultiRegions::e2D:
1094  {
1095  m_curl_dim = 2;
1096  m_bnd_dim = 2;
1097  }
1098  break;
1099  case MultiRegions::e3DH1D:
1100  {
1101  m_curl_dim = 3;
1102  m_bnd_dim = 2;
1103  }
1104  break;
1105  case MultiRegions::e3DH2D:
1106  {
1107  m_curl_dim = 3;
1108  m_bnd_dim = 1;
1109  }
1110  break;
1111  case MultiRegions::e3D:
1112  {
1113  m_curl_dim = 3;
1114  m_bnd_dim = 3;
1115  }
1116  break;
1117  default:
1118  ASSERTL0(0,"Dimension not supported");
1119  break;
1120  }
1121 
1122 
1123  m_HBCdata = Array<OneD, HBCInfo>(HBCnumber);
1124  StdRegions::StdExpansionSharedPtr elmt;
1125 
1126  switch(m_pressure->GetExpType())
1127  {
1128  case MultiRegions::e2D:
1129  case MultiRegions::e3D:
1130  {
1131  int coeff_count = 0;
1132  int exp_size;
1133  int j=0;
1134  int cnt = 0;
1135  for(int n = 0 ; n < m_PBndConds.num_elements(); ++n)
1136  {
1137  exp_size = m_PBndExp[n]->GetExpSize();
1138 
1139  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"H"))
1140  {
1141  for(int i = 0; i < exp_size; ++i,cnt++)
1142  {
1143  m_HBCdata[j].m_globalElmtID = m_pressureBCtoElmtID[cnt];
1144  elmt = m_fields[pindex]->GetExp(m_HBCdata[j].m_globalElmtID);
1145  m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
1146  m_HBCdata[j].m_physOffset = m_fields[pindex]->GetPhys_Offset(m_HBCdata[j].m_globalElmtID);
1147  m_HBCdata[j].m_bndElmtID = i;
1148  m_HBCdata[j].m_elmtTraceID = m_pressureBCtoTraceID[cnt];
1149  m_HBCdata[j].m_bndryID = n;
1150  m_HBCdata[j].m_coeffOffset = coeff_count;
1151  coeff_count += elmt->GetTraceNcoeffs(m_HBCdata[j].m_elmtTraceID);
1152  j = j+1;
1153  }
1154  }
1155  else // setting if just standard BC no High order
1156  {
1157  cnt += exp_size;
1158  }
1159  }
1160  }
1161  break;
1162 
1163  case MultiRegions::e3DH1D:
1164  {
1165  Array<OneD, unsigned int> planes;
1166  planes = m_pressure->GetZIDs();
1167  int num_planes = planes.num_elements();
1168  int num_elm_per_plane = (m_pressure->GetExpSize())/num_planes;
1169 
1170  m_wavenumber = Array<OneD, NekDouble>(HBCnumber);
1171  m_negWavenumberSq = Array<OneD, NekDouble>(HBCnumber);
1172 
1173  int exp_size, exp_size_per_plane;
1174  int i, j, k, n;
1175  int K;
1176  NekDouble sign = -1.0;
1177  int cnt = 0;
1178 
1179  m_session->MatchSolverInfo("ModeType", "SingleMode",
1180  m_SingleMode, false);
1181  m_session->MatchSolverInfo("ModeType", "HalfMode",
1182  m_HalfMode, false);
1183  m_session->MatchSolverInfo("ModeType", "MultipleModes",
1184  m_MultipleModes, false);
1185  m_session->LoadParameter("LZ", m_LhomZ);
1186 
1187  // Stability Analysis flags
1188  if(m_session->DefinesSolverInfo("ModeType"))
1189  {
1190  if(m_SingleMode)
1191  {
1192  m_npointsZ = 2;
1193  }
1194  else if(m_HalfMode)
1195  {
1196  m_npointsZ = 1;
1197  }
1198  else if(m_MultipleModes)
1199  {
1200  m_npointsZ = m_session->GetParameter("HomModesZ");
1201  }
1202  else
1203  {
1204  ASSERTL0(false, "SolverInfo ModeType not valid");
1205  }
1206  }
1207  else
1208  {
1209  m_npointsZ = m_session->GetParameter("HomModesZ");
1210  }
1211 
1212  int coeff_count = 0;
1213 
1214  for(n = 0, j= 0, cnt = 0; n < m_PBndConds.num_elements(); ++n)
1215  {
1216  exp_size = m_PBndExp[n]->GetExpSize();
1217  exp_size_per_plane = exp_size/num_planes;
1218 
1219  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"H"))
1220  {
1221  for(k = 0; k < num_planes; k++)
1222  {
1223  K = planes[k]/2;
1224  for(i = 0; i < exp_size_per_plane; ++i, ++j, ++cnt)
1225  {
1226  m_HBCdata[j].m_globalElmtID = m_pressureBCtoElmtID[cnt];
1227  elmt = m_fields[pindex]->GetExp(m_HBCdata[j].m_globalElmtID);
1228  m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
1229  m_HBCdata[j].m_physOffset = m_fields[pindex]->GetPhys_Offset(m_HBCdata[j].m_globalElmtID);
1230  m_HBCdata[j].m_bndElmtID = i+k*exp_size_per_plane;
1231  m_HBCdata[j].m_elmtTraceID = m_pressureBCtoTraceID[cnt];
1232  m_HBCdata[j].m_bndryID = n;
1233  m_HBCdata[j].m_coeffOffset = coeff_count;
1234  coeff_count += elmt->GetEdgeNcoeffs(m_HBCdata[j].m_elmtTraceID);
1235 
1236  if(m_SingleMode)
1237  {
1238  m_wavenumber[j] = -2*M_PI/m_LhomZ;
1239  m_negWavenumberSq[j] = -1.0*m_wavenumber[j]*m_wavenumber[j];
1240  }
1241  else if(m_HalfMode || m_MultipleModes)
1242  {
1243  m_wavenumber[j] = 2*M_PI/m_LhomZ;
1244  m_negWavenumberSq[j] = -1.0*m_wavenumber[j]*m_wavenumber[j];
1245  }
1246  else
1247  {
1248  m_wavenumber[j] = 2*M_PI*sign*(NekDouble(K))/m_LhomZ;
1249  m_negWavenumberSq[j] = -1.0*m_wavenumber[j]*m_wavenumber[j];
1250  }
1251 
1252  int assElmtID;
1253 
1254  if(k%2==0)
1255  {
1256  if(m_HalfMode)
1257  {
1258  assElmtID = m_HBCdata[j].m_globalElmtID;
1259 
1260  }
1261  else
1262  {
1263  assElmtID = m_HBCdata[j].m_globalElmtID + num_elm_per_plane;
1264  }
1265  }
1266  else
1267  {
1268  assElmtID = m_HBCdata[j].m_globalElmtID - num_elm_per_plane;
1269  }
1270 
1271  m_HBCdata[j].m_assPhysOffset = m_pressure->GetPhys_Offset(assElmtID);
1272  }
1273  sign = -1.0*sign;
1274  }
1275  }
1276  else
1277  {
1278  cnt += exp_size;
1279  }
1280  }
1281  }
1282  break;
1283 
1284  case MultiRegions::e3DH2D:
1285  {
1286  int cnt = 0;
1287  int exp_size, exp_size_per_line;
1288  int j=0;
1289 
1290  for(int k1 = 0; k1 < m_npointsZ; k1++)
1291  {
1292  for(int k2 = 0; k2 < m_npointsY; k2++)
1293  {
1294  for(int n = 0 ; n < m_PBndConds.num_elements(); ++n)
1295  {
1296  exp_size = m_PBndExp[n]->GetExpSize();
1297 
1298  exp_size_per_line = exp_size/(m_npointsZ*m_npointsY);
1299 
1300  if(boost::iequals(m_PBndConds[n]->GetUserDefined(),"H"))
1301  {
1302  for(int i = 0; i < exp_size_per_line; ++i,cnt++)
1303  {
1304  // find element and edge of this expansion.
1305  m_HBCdata[j].m_globalElmtID = m_pressureBCtoElmtID[cnt];
1306  elmt = m_fields[pindex]->GetExp(m_HBCdata[j].m_globalElmtID);
1307  m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
1308  m_HBCdata[j].m_physOffset = m_fields[pindex]->GetPhys_Offset(m_HBCdata[j].m_globalElmtID);
1309  m_HBCdata[j].m_bndElmtID = i+(k1*m_npointsY+k2)*exp_size_per_line;
1310  m_HBCdata[j].m_elmtTraceID = m_pressureBCtoTraceID[cnt];
1311  m_HBCdata[j].m_bndryID = n;
1312  //m_wavenumber[j] = 2*M_PI*sign*(NekDouble(k1))/m_LhomZ;
1313  //m_negWavenumberSq[j] = 2*M_PI*sign*(NekDouble(k2))/m_LhomY;
1314  }
1315  }
1316  else
1317  {
1318  cnt += exp_size_per_line;
1319  }
1320  }
1321  }
1322  }
1323  }
1324  break;
1325  default:
1326  ASSERTL0(0,"Dimension not supported");
1327  break;
1328  }
1329  }
1330 
1331  /**
1332  *
1333  */
1334  Array<OneD, NekDouble> Extrapolate::GetMaxStdVelocity(
1335  const Array<OneD, Array<OneD,NekDouble> > inarray)
1336  {
1337  // Checking if the problem is 2D
1338  ASSERTL0(m_curl_dim >= 2, "Method not implemented for 1D");
1339 
1340  int n_points_0 = m_fields[0]->GetExp(0)->GetTotPoints();
1341  int n_element = m_fields[0]->GetExpSize();
1342  int nvel = inarray.num_elements();
1343  int cnt;
1344 
1345  NekDouble pntVelocity;
1346 
1347  // Getting the standard velocity vector on the 2D normal space
1348  Array<OneD, Array<OneD, NekDouble> > stdVelocity(nvel);
1349  Array<OneD, NekDouble> maxV(n_element, 0.0);
1350  LibUtilities::PointsKeyVector ptsKeys;
1351 
1352  for (int i = 0; i < nvel; ++i)
1353  {
1354  stdVelocity[i] = Array<OneD, NekDouble>(n_points_0);
1355  }
1356 
1357  if (nvel == 2)
1358  {
1359  cnt = 0.0;
1360  for (int el = 0; el < n_element; ++el)
1361  {
1362  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
1363  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
1364 
1365  // reset local space if necessary
1366  if(n_points != n_points_0)
1367  {
1368  for (int j = 0; j < nvel; ++j)
1369  {
1370  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
1371  }
1372  n_points_0 = n_points;
1373  }
1374 
1375  Array<TwoD, const NekDouble> gmat =
1376  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1377 
1378  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
1379  == SpatialDomains::eDeformed)
1380  {
1381  for (int i = 0; i < n_points; i++)
1382  {
1383  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
1384  + gmat[2][i]*inarray[1][i+cnt];
1385 
1386  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
1387  + gmat[3][i]*inarray[1][i+cnt];
1388  }
1389  }
1390  else
1391  {
1392  for (int i = 0; i < n_points; i++)
1393  {
1394  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
1395  + gmat[2][0]*inarray[1][i+cnt];
1396 
1397  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
1398  + gmat[3][0]*inarray[1][i+cnt];
1399  }
1400  }
1401 
1402  cnt += n_points;
1403 
1404 
1405  for (int i = 0; i < n_points; i++)
1406  {
1407  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
1408  + stdVelocity[1][i]*stdVelocity[1][i];
1409 
1410  if (pntVelocity>maxV[el])
1411  {
1412  maxV[el] = pntVelocity;
1413  }
1414  }
1415  maxV[el] = sqrt(maxV[el]);
1416  }
1417  }
1418  else
1419  {
1420  cnt = 0;
1421  for (int el = 0; el < n_element; ++el)
1422  {
1423 
1424  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
1425  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
1426 
1427  // reset local space if necessary
1428  if(n_points != n_points_0)
1429  {
1430  for (int j = 0; j < nvel; ++j)
1431  {
1432  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
1433  }
1434  n_points_0 = n_points;
1435  }
1436 
1437  Array<TwoD, const NekDouble> gmat =
1438  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1439 
1440  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
1441  == SpatialDomains::eDeformed)
1442  {
1443  for (int i = 0; i < n_points; i++)
1444  {
1445  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
1446  + gmat[3][i]*inarray[1][i+cnt]
1447  + gmat[6][i]*inarray[2][i+cnt];
1448 
1449  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
1450  + gmat[4][i]*inarray[1][i+cnt]
1451  + gmat[7][i]*inarray[2][i+cnt];
1452 
1453  stdVelocity[2][i] = gmat[2][i]*inarray[0][i+cnt]
1454  + gmat[5][i]*inarray[1][i+cnt]
1455  + gmat[8][i]*inarray[2][i+cnt];
1456  }
1457  }
1458  else
1459  {
1460  for (int i = 0; i < n_points; i++)
1461  {
1462  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
1463  + gmat[3][0]*inarray[1][i+cnt]
1464  + gmat[6][0]*inarray[2][i+cnt];
1465 
1466  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
1467  + gmat[4][0]*inarray[1][i+cnt]
1468  + gmat[7][0]*inarray[2][i+cnt];
1469 
1470  stdVelocity[2][i] = gmat[2][0]*inarray[0][i+cnt]
1471  + gmat[5][0]*inarray[1][i+cnt]
1472  + gmat[8][0]*inarray[2][i+cnt];
1473  }
1474  }
1475 
1476  cnt += n_points;
1477 
1478  for (int i = 0; i < n_points; i++)
1479  {
1480  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
1481  + stdVelocity[1][i]*stdVelocity[1][i]
1482  + stdVelocity[2][i]*stdVelocity[2][i];
1483 
1484  if (pntVelocity > maxV[el])
1485  {
1486  maxV[el] = pntVelocity;
1487  }
1488  }
1489 
1490  maxV[el] = sqrt(maxV[el]);
1491  }
1492  }
1493 
1494  return maxV;
1495  }
1496 
1497 
1498  LibUtilities::TimeIntegrationMethod Extrapolate::v_GetSubStepIntegrationMethod(void)
1499  {
1500  return LibUtilities::eNoTimeIntegrationMethod;
1501  }
1502 
1503  /**
1504  * Update oldarrays to include newarray and
1505  * extrapolate result to outarray
1506  */
1507  void Extrapolate::ExtrapolateArray(
1508  Array<OneD, Array<OneD, NekDouble> > &oldarrays,
1509  Array<OneD, NekDouble> &newarray,
1510  Array<OneD, NekDouble> &outarray)
1511  {
1512  int nint = min(m_pressureCalls,m_intSteps);
1513  int nPts = newarray.num_elements();
1514 
1515  // Update oldarrays
1516  RollOver(oldarrays);
1517  Vmath::Vcopy(nPts, newarray, 1, oldarrays[0], 1);
1518 
1519  // Extrapolate to outarray
1520  Vmath::Smul(nPts, StifflyStable_Betaq_Coeffs[nint-1][nint-1],
1521  oldarrays[nint-1], 1,
1522  outarray, 1);
1523 
1524  for(int n = 0; n < nint-1; ++n)
1525  {
1526  Vmath::Svtvp(nPts, StifflyStable_Betaq_Coeffs[nint-1][n],
1527  oldarrays[n],1,outarray,1,
1528  outarray,1);
1529  }
1530  }
1531 }
1532 
Extrapolate.h
Nektar::Extrapolate::m_wavenumber
Array< OneD, NekDouble > m_wavenumber
wave number 2 pi k /Lz
Definition: Extrapolate.h:280
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::Extrapolate::m_pressure
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:215
Nektar::Extrapolate::IProductNormVelocityOnHBC
void IProductNormVelocityOnHBC(const Array< OneD, const Array< OneD, NekDouble > > &Vel, Array< OneD, NekDouble > &IprodVn)
Definition: Extrapolate.cpp:921
Nektar::Extrapolate::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:208
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::Extrapolate::m_velocity
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition: Extrapolate.h:219
Nektar::Extrapolate::m_pressureBCsMaxPts
int m_pressureBCsMaxPts
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:241
Nektar::Extrapolate::m_PBndExp
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:235
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
Nektar::GetExtrapolateFactory
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:50
Nektar::Extrapolate::m_SingleMode
bool m_SingleMode
Flag to determine if single homogeneous mode is used.
Definition: Extrapolate.h:252
Nektar::Extrapolate::m_pressureBCsElmtMaxPts
int m_pressureBCsElmtMaxPts
Maximum points used in Element adjacent to pressure BC evaluation.
Definition: Extrapolate.h:244
Vmath::Svtvp
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:471
Nektar::SolverUtils::AdvectionSharedPtr
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:158
Vmath::Vvtvp
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
Nektar::Extrapolate::RollOver
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
Definition: Extrapolate.cpp:1012
Nektar::Extrapolate::m_PBndCoeffs
Array< OneD, NekDouble > m_PBndCoeffs
(if homogeneous)
Definition: Extrapolate.h:304
std
STL namespace.
Nektar::Extrapolate::m_UBndCoeffs
Array< OneD, Array< OneD, NekDouble > > m_UBndCoeffs
(if homogeneous)
Definition: Extrapolate.h:307
Nektar::Extrapolate::CopyPressureHBCsToPbndExp
void CopyPressureHBCsToPbndExp(void)
Definition: Extrapolate.cpp:141
Nektar::Extrapolate::m_nonlinearterm_phys
Array< OneD, NekDouble > m_nonlinearterm_phys
(if homogeneous)
Definition: Extrapolate.h:294
Nektar::Extrapolate::m_npointsY
int m_npointsY
number of points in Y direction (if homogeneous)
Definition: Extrapolate.h:261
Nektar::Extrapolate::ExtrapolatePressureHBCs
void ExtrapolatePressureHBCs(void)
Definition: Extrapolate.cpp:122
Nektar::Extrapolate::m_comm
LibUtilities::CommSharedPtr m_comm
Definition: Extrapolate.h:210
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::Extrapolate::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Extrapolate.h:212
Nektar::Extrapolate::StifflyStable_Gamma0_Coeffs
static NekDouble StifflyStable_Gamma0_Coeffs[3]
Definition: Extrapolate.h:314
Nektar::Extrapolate::m_negWavenumberSq
Array< OneD, NekDouble > m_negWavenumberSq
minus Square of wavenumber
Definition: Extrapolate.h:283
Nektar::Extrapolate::m_nonlinearterm_coeffs
Array< OneD, NekDouble > m_nonlinearterm_coeffs
(if homogeneous)
Definition: Extrapolate.h:297
Nektar::Extrapolate::m_totexps_per_plane
int m_totexps_per_plane
(if homogeneous)
Definition: Extrapolate.h:309
Nektar::Extrapolate::m_outflowVel
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_outflowVel
Storage for current and previous velocity fields at the otuflow for high order outflow BCs...
Definition: Extrapolate.h:286
Nektar::Extrapolate::def
static std::string def
Definition: Extrapolate.h:317
Nektar::Extrapolate::m_bnd_dim
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:229
Nektar::Extrapolate::m_pressureBCtoElmtID
Array< OneD, int > m_pressureBCtoElmtID
Id of element to which pressure boundary condition belongs.
Definition: Extrapolate.h:265
Nektar::Extrapolate::MountHOPBCs
void MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
Definition: Extrapolate.h:376
Nektar::StdRegions::StdExpansion
The base class for all shapes.
Definition: StdExpansion.h:69
Vmath::Smul
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
Nektar::Extrapolate::m_acceleration
Array< OneD, Array< OneD, NekDouble > > m_acceleration
Storage for current and previous levels of the acceleration term.
Definition: Extrapolate.h:274
Nektar::Extrapolate::GetMaxStdVelocity
Array< OneD, NekDouble > GetMaxStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
Definition: Extrapolate.cpp:1334
Nektar::Extrapolate::v_GetSubStepIntegrationMethod
virtual LibUtilities::TimeIntegrationMethod v_GetSubStepIntegrationMethod(void)
Definition: Extrapolate.cpp:1498
Nektar::Extrapolate::m_intSteps
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:247
Nektar::Extrapolate::m_LhomZ
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition: Extrapolate.h:258
Nektar::Extrapolate::m_MultipleModes
bool m_MultipleModes
Flag to determine if use multiple homogenenous modes are used.
Definition: Extrapolate.h:256
Nektar::ExtrapolateFactory
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:79
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::Extrapolate::m_expsize_per_plane
Array< OneD, unsigned int > m_expsize_per_plane
(if homogeneous)
Definition: Extrapolate.h:301
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::Extrapolate::CalcExplicitDuDt
void CalcExplicitDuDt(const Array< OneD, const Array< OneD, NekDouble > > &fields)
Definition: Extrapolate.cpp:83
Nektar::Extrapolate::m_pressureBCtoTraceID
Array< OneD, int > m_pressureBCtoTraceID
Id of edge (2D) or face (3D) to which pressure boundary condition belongs.
Definition: Extrapolate.h:268
Nektar::Extrapolate::m_pressureCalls
int m_pressureCalls
number of times the high-order pressure BCs have been called
Definition: Extrapolate.h:238
Nektar::Extrapolate::v_CalcNeumannPressureBCs
virtual void v_CalcNeumannPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
Definition: Extrapolate.cpp:163
Nektar::Extrapolate::IProductNormVelocityBCOnHBC
void IProductNormVelocityBCOnHBC(Array< OneD, NekDouble > &IprodVn)
Definition: Extrapolate.cpp:962
Vmath::Vsub
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
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Comm.h
Nektar::Extrapolate::ExtrapolateArray
void ExtrapolateArray(Array< OneD, Array< OneD, NekDouble > > &oldarrays, Array< OneD, NekDouble > &newarray, Array< OneD, NekDouble > &outarray)
Definition: Extrapolate.cpp:1507
Nektar::Extrapolate::m_PhyoutfVel
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_PhyoutfVel
Storage for current and previous velocity fields in physical space at the otuflow for high order outf...
Definition: Extrapolate.h:291
Nektar::Extrapolate::v_CorrectPressureBCs
virtual void v_CorrectPressureBCs(const Array< OneD, NekDouble > &pressure)
Definition: Extrapolate.cpp:266
Nektar::Extrapolate::m_pressureHBCs
Array< OneD, Array< OneD, NekDouble > > m_pressureHBCs
Storage for current and previous levels of high order pressure boundary conditions.
Definition: Extrapolate.h:271
Nektar::Extrapolate::m_HBCdata
Array< OneD, HBCInfo > m_HBCdata
data structure to old all the information regarding High order pressure BCs
Definition: Extrapolate.h:277
Nektar::Extrapolate::StifflyStable_Alpha_Coeffs
static NekDouble StifflyStable_Alpha_Coeffs[3][3]
Definition: Extrapolate.h:313
Nektar::Extrapolate::m_npointsZ
int m_npointsZ
number of points in Z direction (if homogeneous)
Definition: Extrapolate.h:262
Nektar::LibUtilities::SessionReader::RegisterDefaultSolverInfo
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
Definition: SessionReader.h:664
Nektar::Extrapolate::CurlCurl
void CurlCurl(Array< OneD, Array< OneD, const NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q, const int j)
Definition: Extrapolate.cpp:792
Nektar::Extrapolate::m_curl_dim
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:226
Nektar::Extrapolate::m_HalfMode
bool m_HalfMode
Flag to determine if half homogeneous mode is used.
Definition: Extrapolate.h:254
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::Extrapolate::CalcOutflowBCs
void CalcOutflowBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, NekDouble kinvis)
Definition: Extrapolate.cpp:270
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873
Nektar::Extrapolate::m_PBndConds
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:232
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
Nektar::Extrapolate::m_timestep
NekDouble m_timestep
Definition: Extrapolate.h:249
Nektar::Extrapolate::~Extrapolate
virtual ~Extrapolate()
Definition: Extrapolate.cpp:75
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Nektar::SolverUtils::Advection
Defines a callback function which evaluates the flux vector.
Definition: Advection.h:69
Nektar::Extrapolate::StifflyStable_Betaq_Coeffs
static NekDouble StifflyStable_Betaq_Coeffs[3][3]
total number of expansion for each plane (if homogeneous)
Definition: Extrapolate.h:312
Nektar::LibUtilities::NekFactory
Provides a generic Factory class.
Definition: NekFactory.hpp:116
Generated on Fri Dec 2 2016 13:10:47 for Nektar++ by   doxygen 1.8.8