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