Nektar++
Public Member Functions | Protected Member Functions | Protected Attributes | Static Private Attributes | List of all members
Nektar::Extrapolate Class Referenceabstract

#include <Extrapolate.h>

Inheritance diagram for Nektar::Extrapolate:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Extrapolate:
Collaboration graph
[legend]

Public Member Functions

 Extrapolate (const LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields, MultiRegions::ExpListSharedPtr pPressure, const Array< OneD, int > pVel, const SolverUtils::AdvectionSharedPtr advObject)
 
virtual ~Extrapolate ()
 
void GenerateHOPBCMap ()
 
void SubSteppingTimeIntegration (const int intMethod, const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
 
void SubStepSaveFields (const int nstep)
 
void SubStepSetPressureBCs (const Array< OneD, const Array< OneD, NekDouble > > &inarray, const NekDouble Aii_DT, NekDouble kinvis)
 
void SubStepAdvance (const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln, const int nstep, NekDouble time)
 
void MountHOPBCs (int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
 
void EvaluatePressureBCs (const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
 
Array< OneD, NekDoubleGetMaxStdVelocity (const Array< OneD, Array< OneD, NekDouble > > inarray)
 

Protected Member Functions

virtual void v_SubSteppingTimeIntegration (int intMethod, const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)=0
 
virtual void v_SubStepSaveFields (int nstep)=0
 
virtual void v_SubStepSetPressureBCs (const Array< OneD, const Array< OneD, NekDouble > > &inarray, NekDouble Aii_DT, NekDouble kinvis)=0
 
virtual void v_SubStepAdvance (const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln, int nstep, NekDouble time)=0
 
virtual void v_MountHOPBCs (int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)=0
 
void CalcNeumannPressureBCs (const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
 
void CalcOutflowBCs (const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
 
void RollOver (Array< OneD, Array< OneD, NekDouble > > &input)
 
void CurlCurl (Array< OneD, Array< OneD, const NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q, const int j)
 

Protected Attributes

LibUtilities::SessionReaderSharedPtr m_session
 
LibUtilities::CommSharedPtr m_comm
 
Array< OneD, MultiRegions::ExpListSharedPtrm_fields
 
MultiRegions::ExpListSharedPtr m_pressure
 Pointer to field holding pressure field. More...
 
Array< OneD, int > m_velocity
 int which identifies which components of m_fields contains the velocity (u,v,w); More...
 
SolverUtils::AdvectionSharedPtr m_advObject
 
Array< OneD, Array< OneD, NekDouble > > m_previousVelFields
 
int m_curl_dim
 Curl-curl dimensionality. More...
 
int m_bnd_dim
 bounday dimensionality More...
 
Array< OneD, const SpatialDomains::BoundaryConditionShPtrm_PBndConds
 pressure boundary conditions container More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_PBndExp
 pressure boundary conditions expansion container More...
 
int m_pressureCalls
 number of times the high-order pressure BCs have been called More...
 
int m_pressureBCsMaxPts
 Maximum points used in pressure BC evaluation. More...
 
int m_pressureBCsElmtMaxPts
 Maximum points used in Element adjacent to pressure BC evaluation. More...
 
int m_intSteps
 Maximum points used in pressure BC evaluation. More...
 
NekDouble m_timestep
 
bool m_SingleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_HalfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_MultipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
Array< OneD, int > m_pressureBCtoElmtID
 Id of element to which pressure boundary condition belongs. More...
 
Array< OneD, int > m_pressureBCtoTraceID
 Id of edge (2D) or face (3D) to which pressure boundary condition belongs. More...
 
Array< OneD, Array< OneD, NekDouble > > m_pressureHBCs
 Storage for current and previous levels of high order pressure boundary conditions. More...
 
Array< OneD, Array< OneD, NekDouble > > m_acceleration
 Storage for current and previous levels of the acceleration term. More...
 
Array< OneD, HBCInfom_HBCdata
 data structure to old all the information regarding High order pressure BCs More...
 
Array< OneD, NekDoublem_wavenumber
 wave number 2 pi k /Lz More...
 
Array< OneD, NekDoublem_negWavenumberSq
 minus Square of wavenumber More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_outflowVel
 Storage for current and previous velocity fields at the otuflow for high order outflow BCs. More...
 
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 outflow BCs. More...
 
Array< OneD, NekDoublem_nonlinearterm_phys
 (if homogeneous) More...
 
Array< OneD, NekDoublem_nonlinearterm_coeffs
 (if homogeneous) More...
 
Array< OneD, unsigned int > m_expsize_per_plane
 (if homogeneous) More...
 
Array< OneD, NekDoublem_PBndCoeffs
 (if homogeneous) More...
 
Array< OneD, Array< OneD, NekDouble > > m_UBndCoeffs
 (if homogeneous) More...
 
int m_totexps_per_plane
 (if homogeneous) More...
 

Static Private Attributes

static std::string def
 total number of expansion for each plane (if homogeneous) More...
 
static NekDouble StifflyStable_Betaq_Coeffs [3][3]
 
static NekDouble StifflyStable_Alpha_Coeffs [3][3]
 
static NekDouble StifflyStable_Gamma0_Coeffs [3]
 

Detailed Description

Definition at line 83 of file Extrapolate.h.

Constructor & Destructor Documentation

Nektar::Extrapolate::Extrapolate ( const LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr pFields,
MultiRegions::ExpListSharedPtr  pPressure,
const Array< OneD, int >  pVel,
const SolverUtils::AdvectionSharedPtr  advObject 
)

Definition at line 57 of file Extrapolate.cpp.

References m_comm, m_session, and m_timestep.

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  }
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:177
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:170
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition: Extrapolate.h:181
LibUtilities::CommSharedPtr m_comm
Definition: Extrapolate.h:172
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Extrapolate.h:174
SolverUtils::AdvectionSharedPtr m_advObject
Definition: Extrapolate.h:183
NekDouble m_timestep
Definition: Extrapolate.h:211
Nektar::Extrapolate::~Extrapolate ( )
virtual

Definition at line 73 of file Extrapolate.cpp.

74  {
75  }

Member Function Documentation

void Nektar::Extrapolate::CalcNeumannPressureBCs ( const Array< OneD, const Array< OneD, NekDouble > > &  fields,
const Array< OneD, const Array< OneD, NekDouble > > &  N,
NekDouble  kinvis 
)
protected

Unified routine for calculation high-oder terms

Casting the boundary expansion to the specific case

Picking up the element where the HOPBc is located

Assigning

Calculating the curl-curl and storing it in Q

Definition at line 191 of file Extrapolate.cpp.

References ASSERTL0, CurlCurl(), Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::MultiRegions::e3DH1D, Nektar::MultiRegions::e3DH2D, m_acceleration, m_bnd_dim, m_curl_dim, m_HBCdata, m_PBndExp, m_pressure, m_pressureBCsElmtMaxPts, m_pressureBCsMaxPts, and MountHOPBCs().

Referenced by EvaluatePressureBCs().

195  {
200 
203 
206 
207  for(int i = 0; i < m_bnd_dim; i++)
208  {
209  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,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:
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;
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:177
int m_pressureBCsMaxPts
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:203
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:197
int m_pressureBCsElmtMaxPts
Maximum points used in Element adjacent to pressure BC evaluation.
Definition: Extrapolate.h:206
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:191
void MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
Definition: Extrapolate.h:324
The base class for all shapes.
Definition: StdExpansion.h:69
Array< OneD, Array< OneD, NekDouble > > m_acceleration
Storage for current and previous levels of the acceleration term.
Definition: Extrapolate.h:236
Array< OneD, HBCInfo > m_HBCdata
data structure to old all the information regarding High order pressure BCs
Definition: Extrapolate.h:239
void CurlCurl(Array< OneD, Array< OneD, const NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q, const int j)
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:188
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void Nektar::Extrapolate::CalcOutflowBCs ( const Array< OneD, const Array< OneD, NekDouble > > &  fields,
const Array< OneD, const Array< OneD, NekDouble > > &  N,
NekDouble  kinvis 
)
protected

Definition at line 329 of file Extrapolate.cpp.

References Nektar::MultiRegions::DirCartesianMap, Nektar::MultiRegions::e3DH1D, m_bnd_dim, m_curl_dim, m_expsize_per_plane, m_fields, m_intSteps, m_nonlinearterm_coeffs, m_nonlinearterm_phys, m_outflowVel, m_PBndCoeffs, m_PBndConds, m_PBndExp, m_PhyoutfVel, m_pressure, m_pressureBCsElmtMaxPts, m_pressureBCsMaxPts, m_pressureBCtoElmtID, m_pressureBCtoTraceID, m_pressureCalls, m_session, m_totexps_per_plane, m_UBndCoeffs, m_velocity, RollOver(), Vmath::Smul(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vvtvp(), and Vmath::Zero().

Referenced by EvaluatePressureBCs().

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 
361  for(int i = 0; i < m_curl_dim; ++i)
362  {
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  {
374 
375  for(int i = 0; i < m_curl_dim; ++i)
376  {
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 
387 
388  m_PBndCoeffs = Array<OneD, NekDouble> (totbndpts,0.0);
390  for(int i = 0; i < m_curl_dim; ++i)
391  {
392  m_UBndCoeffs[i] = Array<OneD, NekDouble> (totbndpts);
393  }
395  planes = m_pressure->GetZIDs();
396  int num_planes = planes.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  }
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 
413  UBndConds(m_curl_dim);
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);
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  {
433  0.0);
434  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,0.0);
435  nGradu[i] = nGradu[0] + i*m_pressureBCsMaxPts;
437 
438  if (m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
439  {
441  }
442  }
443 
444  int nbc,cnt,cnt_start;
445  int veloffset = 0;
446  int nint = min(m_pressureCalls,m_intSteps);
447 
449  Array<OneD, NekDouble> PBCvals, UBCvals;
450  Array<OneD, Array<OneD, NekDouble> > ubc(m_curl_dim);
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(
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 
845  (m_PBndExp[n]->GetExp(i));
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  }
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:177
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:170
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition: Extrapolate.h:181
int m_pressureBCsMaxPts
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:203
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:197
int m_pressureBCsElmtMaxPts
Maximum points used in Element adjacent to pressure BC evaluation.
Definition: Extrapolate.h:206
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
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
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
Array< OneD, NekDouble > m_PBndCoeffs
(if homogeneous)
Definition: Extrapolate.h:264
Array< OneD, Array< OneD, NekDouble > > m_UBndCoeffs
(if homogeneous)
Definition: Extrapolate.h:267
Array< OneD, NekDouble > m_nonlinearterm_phys
(if homogeneous)
Definition: Extrapolate.h:254
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Extrapolate.h:174
Array< OneD, NekDouble > m_nonlinearterm_coeffs
(if homogeneous)
Definition: Extrapolate.h:257
int m_totexps_per_plane
(if homogeneous)
Definition: Extrapolate.h:269
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
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:191
Array< OneD, int > m_pressureBCtoElmtID
Id of element to which pressure boundary condition belongs.
Definition: Extrapolate.h:227
The base class for all shapes.
Definition: StdExpansion.h:69
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
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:209
Array< OneD, unsigned int > m_expsize_per_plane
(if homogeneous)
Definition: Extrapolate.h:261
double NekDouble
Array< OneD, int > m_pressureBCtoTraceID
Id of edge (2D) or face (3D) to which pressure boundary condition belongs.
Definition: Extrapolate.h:230
int m_pressureCalls
number of times the high-order pressure BCs have been called
Definition: Extrapolate.h:200
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
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
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:188
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:194
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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
void Nektar::Extrapolate::CurlCurl ( Array< OneD, Array< OneD, const NekDouble > > &  Vel,
Array< OneD, Array< OneD, NekDouble > > &  Q,
const int  j 
)
protected

Curl Curl routine - dimension dependent

Definition at line 894 of file Extrapolate.cpp.

References ASSERTL0, Nektar::MultiRegions::DirCartesianMap, Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::MultiRegions::e3DH1D, Nektar::MultiRegions::e3DH2D, m_fields, m_HBCdata, m_negWavenumberSq, m_pressureBCsElmtMaxPts, m_wavenumber, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vsub().

Referenced by CalcNeumannPressureBCs().

898  {
900  = m_fields[0]->GetExp(m_HBCdata[j].m_globalElmtID);
901 
904 
905  switch(m_fields[0]->GetExpType())
906  {
907  case MultiRegions::e2D:
908  {
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 
923  {
925 
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 
961  {
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  {
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  }
Array< OneD, NekDouble > m_wavenumber
wave number 2 pi k /Lz
Definition: Extrapolate.h:242
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
int m_pressureBCsElmtMaxPts
Maximum points used in Element adjacent to pressure BC evaluation.
Definition: Extrapolate.h:206
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Extrapolate.h:174
Array< OneD, NekDouble > m_negWavenumberSq
minus Square of wavenumber
Definition: Extrapolate.h:245
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
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
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Array< OneD, HBCInfo > m_HBCdata
data structure to old all the information regarding High order pressure BCs
Definition: Extrapolate.h:239
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
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
void Nektar::Extrapolate::EvaluatePressureBCs ( const Array< OneD, const Array< OneD, NekDouble > > &  fields,
const Array< OneD, const Array< OneD, NekDouble > > &  N,
NekDouble  kinvis 
)

Function to extrapolate the new pressure boundary condition. Based on the velocity field and on the advection term. Acceleration term is also computed. This routine is a general one for 2d and 3D application and it can be called directly from velocity correction scheme. Specialisation on dimensionality is redirected to the CalcPressureBCs method.

Definition at line 90 of file Extrapolate.cpp.

References CalcNeumannPressureBCs(), CalcOutflowBCs(), m_acceleration, m_HBCdata, m_intSteps, m_PBndConds, m_PBndExp, m_pressureCalls, m_pressureHBCs, m_timestep, RollOver(), Vmath::Smul(), StifflyStable_Alpha_Coeffs, StifflyStable_Betaq_Coeffs, StifflyStable_Gamma0_Coeffs, Vmath::Svtvp(), and Vmath::Vcopy().

Referenced by Nektar::SubSteppingExtrapolate::v_SubStepSetPressureBCs().

94  {
95  if(m_HBCdata.num_elements()>0)
96  {
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
112 
113  // Rotate acceleration term
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  {
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  }
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:197
void CalcOutflowBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
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
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
static NekDouble StifflyStable_Gamma0_Coeffs[3]
Definition: Extrapolate.h:277
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
Array< OneD, Array< OneD, NekDouble > > m_acceleration
Storage for current and previous levels of the acceleration term.
Definition: Extrapolate.h:236
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:209
int m_pressureCalls
number of times the high-order pressure BCs have been called
Definition: Extrapolate.h:200
Array< OneD, Array< OneD, NekDouble > > m_pressureHBCs
Storage for current and previous levels of high order pressure boundary conditions.
Definition: Extrapolate.h:233
Array< OneD, HBCInfo > m_HBCdata
data structure to old all the information regarding High order pressure BCs
Definition: Extrapolate.h:239
static NekDouble StifflyStable_Alpha_Coeffs[3][3]
Definition: Extrapolate.h:276
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:194
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
NekDouble m_timestep
Definition: Extrapolate.h:211
static NekDouble StifflyStable_Betaq_Coeffs[3][3]
Definition: Extrapolate.h:275
void CalcNeumannPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
void Nektar::Extrapolate::GenerateHOPBCMap ( )

Map to directly locate HOPBCs position and offsets in all scenarios

Definition at line 1051 of file Extrapolate.cpp.

References ASSERTL0, Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::MultiRegions::e3DH1D, Nektar::MultiRegions::e3DH2D, m_acceleration, m_bnd_dim, m_comm, m_curl_dim, m_HalfMode, m_HBCdata, m_intSteps, m_LhomZ, m_MultipleModes, m_negWavenumberSq, m_npointsY, m_npointsZ, m_PBndConds, m_PBndExp, m_pressure, m_pressureBCsElmtMaxPts, m_pressureBCsMaxPts, m_pressureBCtoElmtID, m_pressureBCtoTraceID, m_pressureCalls, m_pressureHBCs, m_session, m_SingleMode, m_wavenumber, Nektar::LibUtilities::ReduceSum, and sign.

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,
1060 
1061  // find the maximum values of points for pressure BC evaluation
1062  m_pressureBCsMaxPts = 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  {
1070  m_PBndExp[n]->GetExp(i)->GetTotPoints());
1072  m_pressure->GetExp(m_pressureBCtoElmtID[cnt++])
1073  ->GetTotPoints());
1074  }
1075  }
1076 
1077  // Storage array for high order pressure BCs
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);
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  {
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);
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  }
Array< OneD, NekDouble > m_wavenumber
wave number 2 pi k /Lz
Definition: Extrapolate.h:242
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:177
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:170
int m_pressureBCsMaxPts
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:203
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:197
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
bool m_SingleMode
Flag to determine if single homogeneous mode is used.
Definition: Extrapolate.h:214
int m_pressureBCsElmtMaxPts
Maximum points used in Element adjacent to pressure BC evaluation.
Definition: Extrapolate.h:206
int m_npointsY
number of points in Y direction (if homogeneous)
Definition: Extrapolate.h:223
LibUtilities::CommSharedPtr m_comm
Definition: Extrapolate.h:172
Array< OneD, NekDouble > m_negWavenumberSq
minus Square of wavenumber
Definition: Extrapolate.h:245
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:191
Array< OneD, int > m_pressureBCtoElmtID
Id of element to which pressure boundary condition belongs.
Definition: Extrapolate.h:227
Array< OneD, Array< OneD, NekDouble > > m_acceleration
Storage for current and previous levels of the acceleration term.
Definition: Extrapolate.h:236
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:209
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition: Extrapolate.h:220
bool m_MultipleModes
Flag to determine if use multiple homogenenous modes are used.
Definition: Extrapolate.h:218
double NekDouble
Array< OneD, int > m_pressureBCtoTraceID
Id of edge (2D) or face (3D) to which pressure boundary condition belongs.
Definition: Extrapolate.h:230
int m_pressureCalls
number of times the high-order pressure BCs have been called
Definition: Extrapolate.h:200
Array< OneD, Array< OneD, NekDouble > > m_pressureHBCs
Storage for current and previous levels of high order pressure boundary conditions.
Definition: Extrapolate.h:233
Array< OneD, HBCInfo > m_HBCdata
data structure to old all the information regarding High order pressure BCs
Definition: Extrapolate.h:239
int m_npointsZ
number of points in Z direction (if homogeneous)
Definition: Extrapolate.h:224
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:188
bool m_HalfMode
Flag to determine if half homogeneous mode is used.
Definition: Extrapolate.h:216
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:194
Array< OneD, NekDouble > Nektar::Extrapolate::GetMaxStdVelocity ( const Array< OneD, Array< OneD, NekDouble > >  inarray)

Definition at line 1369 of file Extrapolate.cpp.

References ASSERTL0, Nektar::SpatialDomains::eDeformed, m_curl_dim, and m_fields.

Referenced by Nektar::SubSteppingExtrapolate::GetSubstepTimeStep().

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);
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 
1411  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1412 
1413  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
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 
1473  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1474 
1475  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Extrapolate.h:174
double NekDouble
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:188
Geometry is curved or has non-constant factors.
void Nektar::Extrapolate::MountHOPBCs ( int  HBCdata,
NekDouble  kinvis,
Array< OneD, NekDouble > &  Q,
Array< OneD, const NekDouble > &  Advection 
)
inline

Definition at line 324 of file Extrapolate.h.

References v_MountHOPBCs().

Referenced by CalcNeumannPressureBCs().

329  {
330  v_MountHOPBCs(HBCdata,kinvis,Q,Advection);
331  }
virtual void v_MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)=0
void Nektar::Extrapolate::RollOver ( Array< OneD, Array< OneD, NekDouble > > &  input)
protected

Function to roll time-level storages to the next step layout. The stored data associated with the oldest time-level (not required anymore) are moved to the top, where they will be overwritten as the solution process progresses.

Definition at line 1031 of file Extrapolate.cpp.

Referenced by CalcOutflowBCs(), and EvaluatePressureBCs().

1032  {
1033  int nlevels = input.num_elements();
1034 
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  }
void Nektar::Extrapolate::SubStepAdvance ( const LibUtilities::TimeIntegrationSolutionSharedPtr integrationSoln,
const int  nstep,
NekDouble  time 
)
inline

Definition at line 313 of file Extrapolate.h.

References v_SubStepAdvance().

317  {
318  v_SubStepAdvance(integrationSoln,nstep, time);
319  }
virtual void v_SubStepAdvance(const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln, int nstep, NekDouble time)=0
void Nektar::Extrapolate::SubSteppingTimeIntegration ( const int  intMethod,
const LibUtilities::TimeIntegrationWrapperSharedPtr IntegrationScheme 
)
inline

Definition at line 283 of file Extrapolate.h.

References v_SubSteppingTimeIntegration().

286  {
287  v_SubSteppingTimeIntegration(intMethod, IntegrationScheme);
288  }
virtual void v_SubSteppingTimeIntegration(int intMethod, const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)=0
void Nektar::Extrapolate::SubStepSaveFields ( const int  nstep)
inline

Definition at line 293 of file Extrapolate.h.

References v_SubStepSaveFields().

295  {
296  v_SubStepSaveFields(nstep);
297  }
virtual void v_SubStepSaveFields(int nstep)=0
void Nektar::Extrapolate::SubStepSetPressureBCs ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
const NekDouble  Aii_DT,
NekDouble  kinvis 
)
inline

Definition at line 302 of file Extrapolate.h.

References v_SubStepSetPressureBCs().

306  {
307  v_SubStepSetPressureBCs(inarray,Aii_DT,kinvis);
308  }
virtual void v_SubStepSetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, NekDouble Aii_DT, NekDouble kinvis)=0
virtual void Nektar::Extrapolate::v_MountHOPBCs ( int  HBCdata,
NekDouble  kinvis,
Array< OneD, NekDouble > &  Q,
Array< OneD, const NekDouble > &  Advection 
)
protectedpure virtual
virtual void Nektar::Extrapolate::v_SubStepAdvance ( const LibUtilities::TimeIntegrationSolutionSharedPtr integrationSoln,
int  nstep,
NekDouble  time 
)
protectedpure virtual
virtual void Nektar::Extrapolate::v_SubSteppingTimeIntegration ( int  intMethod,
const LibUtilities::TimeIntegrationWrapperSharedPtr IntegrationScheme 
)
protectedpure virtual
virtual void Nektar::Extrapolate::v_SubStepSaveFields ( int  nstep)
protectedpure virtual
virtual void Nektar::Extrapolate::v_SubStepSetPressureBCs ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
NekDouble  Aii_DT,
NekDouble  kinvis 
)
protectedpure virtual

Member Data Documentation

std::string Nektar::Extrapolate::def
staticprivate
Initial value:
=
"StandardExtrapolate", "StandardExtrapolate")

total number of expansion for each plane (if homogeneous)

Definition at line 272 of file Extrapolate.h.

Array<OneD, Array<OneD, NekDouble> > Nektar::Extrapolate::m_acceleration
protected

Storage for current and previous levels of the acceleration term.

Definition at line 236 of file Extrapolate.h.

Referenced by CalcNeumannPressureBCs(), EvaluatePressureBCs(), and GenerateHOPBCMap().

SolverUtils::AdvectionSharedPtr Nektar::Extrapolate::m_advObject
protected

Definition at line 183 of file Extrapolate.h.

Referenced by Nektar::SubSteppingExtrapolate::SubStepAdvection().

int Nektar::Extrapolate::m_bnd_dim
protected

bounday dimensionality

Definition at line 191 of file Extrapolate.h.

Referenced by CalcNeumannPressureBCs(), CalcOutflowBCs(), and GenerateHOPBCMap().

LibUtilities::CommSharedPtr Nektar::Extrapolate::m_comm
protected
int Nektar::Extrapolate::m_curl_dim
protected
Array<OneD, unsigned int> Nektar::Extrapolate::m_expsize_per_plane
protected

(if homogeneous)

expansion sizes of pressure boundary conditions in each plane at the outflow for high order outflow BCs

Definition at line 261 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::Extrapolate::m_fields
protected
bool Nektar::Extrapolate::m_HalfMode
protected

Flag to determine if half homogeneous mode is used.

Definition at line 216 of file Extrapolate.h.

Referenced by GenerateHOPBCMap().

Array<OneD, HBCInfo > Nektar::Extrapolate::m_HBCdata
protected

data structure to old all the information regarding High order pressure BCs

Definition at line 239 of file Extrapolate.h.

Referenced by CalcNeumannPressureBCs(), CurlCurl(), EvaluatePressureBCs(), and GenerateHOPBCMap().

int Nektar::Extrapolate::m_intSteps
protected
NekDouble Nektar::Extrapolate::m_LhomZ
protected

physical length in Z direction (if homogeneous)

Definition at line 220 of file Extrapolate.h.

Referenced by GenerateHOPBCMap().

bool Nektar::Extrapolate::m_MultipleModes
protected

Flag to determine if use multiple homogenenous modes are used.

Definition at line 218 of file Extrapolate.h.

Referenced by GenerateHOPBCMap().

Array<OneD, NekDouble> Nektar::Extrapolate::m_negWavenumberSq
protected

minus Square of wavenumber

Definition at line 245 of file Extrapolate.h.

Referenced by CurlCurl(), and GenerateHOPBCMap().

Array<OneD, NekDouble> Nektar::Extrapolate::m_nonlinearterm_coeffs
protected

(if homogeneous)

Storage for nonlinear term in wave space at the outflow for high order outflow BCs

Definition at line 257 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

Array<OneD, NekDouble> Nektar::Extrapolate::m_nonlinearterm_phys
protected

(if homogeneous)

Storage for nonlinear term in physical space at the outflow for high order outflow BCs

Definition at line 254 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

int Nektar::Extrapolate::m_npointsX
protected

number of points in X direction (if homogeneous)

Definition at line 222 of file Extrapolate.h.

int Nektar::Extrapolate::m_npointsY
protected

number of points in Y direction (if homogeneous)

Definition at line 223 of file Extrapolate.h.

Referenced by GenerateHOPBCMap().

int Nektar::Extrapolate::m_npointsZ
protected

number of points in Z direction (if homogeneous)

Definition at line 224 of file Extrapolate.h.

Referenced by GenerateHOPBCMap().

Array<OneD, Array<OneD, Array<OneD, NekDouble > > > Nektar::Extrapolate::m_outflowVel
protected

Storage for current and previous velocity fields at the otuflow for high order outflow BCs.

Definition at line 248 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

Array<OneD, NekDouble> Nektar::Extrapolate::m_PBndCoeffs
protected

(if homogeneous)

Storage for Fourier Coeffs of Dirichlet pressure condition from the input file

Definition at line 264 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

Array<OneD, const SpatialDomains::BoundaryConditionShPtr> Nektar::Extrapolate::m_PBndConds
protected

pressure boundary conditions container

Definition at line 194 of file Extrapolate.h.

Referenced by CalcOutflowBCs(), EvaluatePressureBCs(), and GenerateHOPBCMap().

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::Extrapolate::m_PBndExp
protected

pressure boundary conditions expansion container

Definition at line 197 of file Extrapolate.h.

Referenced by CalcNeumannPressureBCs(), CalcOutflowBCs(), EvaluatePressureBCs(), and GenerateHOPBCMap().

Array<OneD, Array<OneD, Array<OneD, NekDouble > > > Nektar::Extrapolate::m_PhyoutfVel
protected

Storage for current and previous velocity fields in physical space at the otuflow for high order outflow BCs.

Definition at line 251 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

MultiRegions::ExpListSharedPtr Nektar::Extrapolate::m_pressure
protected

Pointer to field holding pressure field.

Definition at line 177 of file Extrapolate.h.

Referenced by CalcNeumannPressureBCs(), CalcOutflowBCs(), and GenerateHOPBCMap().

int Nektar::Extrapolate::m_pressureBCsElmtMaxPts
protected

Maximum points used in Element adjacent to pressure BC evaluation.

Definition at line 206 of file Extrapolate.h.

Referenced by CalcNeumannPressureBCs(), CalcOutflowBCs(), CurlCurl(), and GenerateHOPBCMap().

int Nektar::Extrapolate::m_pressureBCsMaxPts
protected
Array<OneD, int> Nektar::Extrapolate::m_pressureBCtoElmtID
protected

Id of element to which pressure boundary condition belongs.

Definition at line 227 of file Extrapolate.h.

Referenced by Nektar::SubSteppingExtrapolate::AddDuDt2D(), Nektar::SubSteppingExtrapolate::AddDuDt3D(), CalcOutflowBCs(), and GenerateHOPBCMap().

Array<OneD, int> Nektar::Extrapolate::m_pressureBCtoTraceID
protected

Id of edge (2D) or face (3D) to which pressure boundary condition belongs.

Definition at line 230 of file Extrapolate.h.

Referenced by Nektar::SubSteppingExtrapolate::AddDuDt2D(), Nektar::SubSteppingExtrapolate::AddDuDt3D(), CalcOutflowBCs(), and GenerateHOPBCMap().

int Nektar::Extrapolate::m_pressureCalls
protected

number of times the high-order pressure BCs have been called

Definition at line 200 of file Extrapolate.h.

Referenced by CalcOutflowBCs(), EvaluatePressureBCs(), and GenerateHOPBCMap().

Array<OneD, Array<OneD, NekDouble> > Nektar::Extrapolate::m_pressureHBCs
protected

Storage for current and previous levels of high order pressure boundary conditions.

Definition at line 233 of file Extrapolate.h.

Referenced by EvaluatePressureBCs(), and GenerateHOPBCMap().

Array<OneD, Array<OneD, NekDouble> > Nektar::Extrapolate::m_previousVelFields
protected

Definition at line 185 of file Extrapolate.h.

LibUtilities::SessionReaderSharedPtr Nektar::Extrapolate::m_session
protected
bool Nektar::Extrapolate::m_SingleMode
protected

Flag to determine if single homogeneous mode is used.

Definition at line 214 of file Extrapolate.h.

Referenced by GenerateHOPBCMap().

NekDouble Nektar::Extrapolate::m_timestep
protected
int Nektar::Extrapolate::m_totexps_per_plane
protected

(if homogeneous)

Definition at line 269 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

Array<OneD, Array<OneD, NekDouble> > Nektar::Extrapolate::m_UBndCoeffs
protected

(if homogeneous)

Storage for Fourier Coeffs of Neumann velocity condition from the input file

Definition at line 267 of file Extrapolate.h.

Referenced by CalcOutflowBCs().

Array<OneD, int> Nektar::Extrapolate::m_velocity
protected
Array<OneD, NekDouble> Nektar::Extrapolate::m_wavenumber
protected

wave number 2 pi k /Lz

Definition at line 242 of file Extrapolate.h.

Referenced by CurlCurl(), and GenerateHOPBCMap().

NekDouble Nektar::Extrapolate::StifflyStable_Alpha_Coeffs
staticprivate
Initial value:
= {
{ 1.0, 0.0, 0.0},{ 2.0, -0.5, 0.0},{ 3.0, -1.5, 1.0/3.0}}

Definition at line 276 of file Extrapolate.h.

Referenced by EvaluatePressureBCs().

NekDouble Nektar::Extrapolate::StifflyStable_Betaq_Coeffs
staticprivate
Initial value:
= {
{ 1.0, 0.0, 0.0},{ 2.0, -1.0, 0.0},{ 3.0, -3.0, 1.0}}

Definition at line 275 of file Extrapolate.h.

Referenced by EvaluatePressureBCs().

NekDouble Nektar::Extrapolate::StifflyStable_Gamma0_Coeffs
staticprivate
Initial value:
= {
1.0, 1.5, 11.0/6.0}

Definition at line 277 of file Extrapolate.h.

Referenced by EvaluatePressureBCs().