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