Nektar++
IncNavierStokes.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File IncNavierStokes.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Incompressible Navier Stokes class definition built on
32 // ADRBase class
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/algorithm/string.hpp>
37 #include <iomanip>
38 
45 
49 #include <algorithm>
50 #include <fstream>
51 #include <iostream>
52 #include <sstream>
53 
56 #include <tinyxml.h>
57 
58 using namespace std;
59 
60 namespace Nektar
61 {
62 
63 std::string IncNavierStokes::eqTypeLookupIds[6] = {
64  LibUtilities::SessionReader::RegisterEnumValue("EqType", "SteadyStokes",
66  LibUtilities::SessionReader::RegisterEnumValue("EqType", "SteadyOseen",
67  eSteadyOseen),
68  LibUtilities::SessionReader::RegisterEnumValue(
69  "EqType", "SteadyNavierStokes", eSteadyNavierStokes),
70  LibUtilities::SessionReader::RegisterEnumValue(
71  "EqType", "SteadyLinearisedNS", eSteadyLinearisedNS),
72  LibUtilities::SessionReader::RegisterEnumValue(
73  "EqType", "UnsteadyNavierStokes", eUnsteadyNavierStokes),
74  LibUtilities::SessionReader::RegisterEnumValue("EqType", "UnsteadyStokes",
76 
77 /**
78  * Constructor. Creates ...
79  *
80  * \param
81  * \param
82  */
83 IncNavierStokes::IncNavierStokes(
86  : UnsteadySystem(pSession, pGraph), AdvectionSystem(pSession, pGraph),
87  m_SmoothAdvection(false)
88 {
89 }
90 
91 void IncNavierStokes::v_InitObject(bool DeclareField)
92 {
93  AdvectionSystem::v_InitObject(DeclareField);
94 
95  int i, j;
96  int numfields = m_fields.size();
97  std::string velids[] = {"u", "v", "w"};
98 
99  // Set up Velocity field to point to the first m_expdim of m_fields;
101 
102  for (i = 0; i < m_spacedim; ++i)
103  {
104  for (j = 0; j < numfields; ++j)
105  {
106  std::string var = m_boundaryConditions->GetVariable(j);
107  if (boost::iequals(velids[i], var))
108  {
109  m_velocity[i] = j;
110  break;
111  }
112 
113  ASSERTL0(j != numfields, "Failed to find field: " + var);
114  }
115  }
116 
117  // Set up equation type enum using kEquationTypeStr
118  for (i = 0; i < (int)eEquationTypeSize; ++i)
119  {
120  bool match;
121  m_session->MatchSolverInfo("EQTYPE", kEquationTypeStr[i], match, false);
122  if (match)
123  {
125  break;
126  }
127  }
128  ASSERTL0(i != eEquationTypeSize, "EQTYPE not found in SOLVERINFO section");
129 
130  m_session->LoadParameter("Kinvis", m_kinvis);
131 
132  // Default advection type per solver
133  std::string vConvectiveType;
134  switch (m_equationType)
135  {
136  case eUnsteadyStokes:
137  case eSteadyLinearisedNS:
138  vConvectiveType = "NoAdvection";
139  break;
141  case eSteadyNavierStokes:
142  vConvectiveType = "Convective";
143  break;
145  vConvectiveType = "Linearised";
146  break;
147  default:
148  break;
149  }
150 
151  // Check if advection type overridden
152  if (m_session->DefinesTag("AdvectiveType") &&
155  {
156  vConvectiveType = m_session->GetTag("AdvectiveType");
157  }
158 
159  // Initialise advection
161  vConvectiveType, vConvectiveType);
162  m_advObject->InitObject(m_session, m_fields);
163 
164  // Set up arrays for moving reference frame
165  // Note: this must be done before the forcing
166  m_movingFrameProjMat = bnu::identity_matrix<NekDouble>(3, 3);
167  if (DefinedForcing("MovingReferenceFrame"))
168  {
173  }
174 
175  // Forcing terms
176  m_forcing = SolverUtils::Forcing::Load(m_session, shared_from_this(),
178 
179  // check to see if any Robin boundary conditions and if so set
180  // up m_field to boundary condition maps;
184 
185  for (i = 0; i < m_fields.size(); ++i)
186  {
187  bool Set = false;
188 
191  int radpts = 0;
192 
193  BndConds = m_fields[i]->GetBndConditions();
194  BndExp = m_fields[i]->GetBndCondExpansions();
195  for (int n = 0; n < BndConds.size(); ++n)
196  {
197  if (boost::iequals(BndConds[n]->GetUserDefined(), "Radiation"))
198  {
199  ASSERTL0(
200  BndConds[n]->GetBoundaryConditionType() ==
202  "Radiation boundary condition must be of type Robin <R>");
203 
204  if (Set == false)
205  {
206  m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],
208  Set = true;
209  }
210  radpts += BndExp[n]->GetTotPoints();
211  }
212  if (boost::iequals(BndConds[n]->GetUserDefined(),
213  "ZeroNormalComponent"))
214  {
215  ASSERTL0(BndConds[n]->GetBoundaryConditionType() ==
217  "Zero Normal Component boundary condition option must "
218  "be of type Dirichlet <D>");
219 
220  if (Set == false)
221  {
222  m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],
224  Set = true;
225  }
226  }
227  }
228 
230 
231  radpts = 0; // reset to use as a counter
232 
233  for (int n = 0; n < BndConds.size(); ++n)
234  {
235  if (boost::iequals(BndConds[n]->GetUserDefined(), "Radiation"))
236  {
237 
238  int npoints = BndExp[n]->GetNpoints();
239  Array<OneD, NekDouble> x0(npoints, 0.0);
240  Array<OneD, NekDouble> x1(npoints, 0.0);
241  Array<OneD, NekDouble> x2(npoints, 0.0);
242  Array<OneD, NekDouble> tmpArray;
243 
244  BndExp[n]->GetCoords(x0, x1, x2);
245 
246  LibUtilities::Equation coeff =
247  std::static_pointer_cast<
249  ->m_robinPrimitiveCoeff;
250 
251  coeff.Evaluate(x0, x1, x2, m_time,
252  tmpArray = m_fieldsRadiationFactor[i] + radpts);
253  // Vmath::Neg(npoints,tmpArray = m_fieldsRadiationFactor[i]+
254  // radpts,1);
255  radpts += npoints;
256  }
257  }
258  }
259 
260  // Set up maping for womersley BC - and load variables
261  for (int i = 0; i < m_fields.size(); ++i)
262  {
263  for (int n = 0; n < m_fields[i]->GetBndConditions().size(); ++n)
264  {
265  if (boost::istarts_with(
266  m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
267  "Womersley"))
268  {
269  // assumes that boundary condition is applied in normal
270  // direction and is decomposed for each direction. There could
271  // be a unique file for each direction
272  m_womersleyParams[i][n] =
274  m_spacedim);
275  // Read in fourier coeffs and precompute coefficients
277  i, n, m_fields[i]->GetBndConditions()[n]->GetUserDefined());
278 
279  m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],
281  }
282  }
283  }
284 
285  // Check if we have moving reference frame boundary conditions,if so, all
286  // velocity vectors at those boundaries should be Drichlet and have a same
287  // tag loop over the boundaries
288  for (int n = 0; n < m_fields[0]->GetBndConditions().size(); ++n)
289  {
290  // loop over the velocities
291  for (int i = 0; i < m_velocity.size(); ++i)
292  {
293  if (boost::iequals(
294  m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
295  "MovingFrameWall"))
296  {
297  // all velocities must have the same userdefined tag and all
298  // must be Dirichlet
299  for (int j = 0; j < m_velocity.size(); ++j)
300  {
301  ASSERTL0(m_fields[j]
302  ->GetBndConditions()[n]
303  ->GetBoundaryConditionType() ==
305  "All velocities with userdefined tag: "
306  "\"MovingFrameWall\" must be Dirichlet boundary "
307  "condition");
308  ASSERTL0(boost::iequals(m_fields[j]
309  ->GetBndConditions()[n]
310  ->GetUserDefined(),
311  "MovingFrameWall"),
312  "If any of the velocity components at a "
313  "boundary is defined as \"MovingFrameWall\", "
314  "the rest of velocity components are also "
315  "must be defined as \"MovingFrameWall\" ");
316  }
317  }
318 
319  if (boost::iequals(
320  m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
321  "MovingFrameDomainVel"))
322  {
323  // all velocities must have the same userdefined tag and all
324  // must be Dirichlet
325  for (int j = 0; j < m_velocity.size(); ++j)
326  {
327  ASSERTL0(
328  m_fields[j]
329  ->GetBndConditions()[n]
330  ->GetBoundaryConditionType() ==
332  "All velocities with userdefined tag: "
333  "\"MovingFrameDomainVel\" must be Dirichlet boundary "
334  "condition");
335  ASSERTL0(boost::iequals(m_fields[j]
336  ->GetBndConditions()[n]
337  ->GetUserDefined(),
338  "MovingFrameDomainVel"),
339  "If any of the velocity components at a "
340  "boundary is defined as \"MovingFrameDomainVel\", "
341  "the rest of velocity components are also "
342  "must be defined as \"MovingFrameDomainVel\" ");
343  }
344  }
345  }
346  }
347 
348  // Set up Field Meta Data for output files
349  m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
350  m_fieldMetaDataMap["TimeStep"] =
351  boost::lexical_cast<std::string>(m_timestep);
352 
353  // Check to see if the metadata for moving reference frame is defined
354  if (DefinedForcing("MovingReferenceFrame"))
355  {
357  {
358  std::vector<std::string> vSuffix = {"_x", "_y", "_z"};
359  for (int i = 0; i < 3; ++i)
360  {
361  NekDouble dTheta{0};
362  std::string sTheta = "Theta" + vSuffix[i];
363  auto it = m_fieldMetaDataMap.find(sTheta);
364  if (it != m_fieldMetaDataMap.end())
365  {
366  dTheta = boost::lexical_cast<NekDouble>(it->second);
367  m_movingFrameTheta[i] = dTheta;
368  }
369  else
370  {
371  m_fieldMetaDataMap[sTheta] =
372  boost::lexical_cast<std::string>(m_movingFrameTheta[i]);
373  }
374  }
375  }
376  else
377  {
378  std::vector<std::string> vSuffix = {"_x", "_y", "_z"};
379  for (int i = 0; i < 3; ++i)
380  {
381  std::string sTheta = "Theta" + vSuffix[i];
382  m_fieldMetaDataMap[sTheta] =
383  boost::lexical_cast<std::string>(m_movingFrameTheta[i]);
384  }
385  }
386  }
387 }
388 
389 /**
390  * Destructor
391  */
393 {
394 }
395 
396 /**
397  * Evaluation -N(V) for all fields except pressure using m_velocity
398  */
400  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
401  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
402 {
403  int i;
404  int VelDim = m_velocity.size();
405  Array<OneD, Array<OneD, NekDouble>> velocity(VelDim);
406 
407  int npoints = m_fields[0]->GetNpoints();
408  for (i = 0; i < VelDim; ++i)
409  {
410  velocity[i] = Array<OneD, NekDouble>(npoints);
411  Vmath::Vcopy(npoints, inarray[m_velocity[i]], 1, velocity[i], 1);
412  }
413  for (auto &x : m_forcing)
414  {
415  x->PreApply(m_fields, velocity, velocity, time);
416  }
417 
418  m_advObject->Advect(m_nConvectiveFields, m_fields, velocity, inarray,
419  outarray, time);
420 }
421 
422 /**
423  * Time dependent boundary conditions updating
424  */
426 {
427  int i, n;
428  std::string varName;
429  int nvariables = m_fields.size();
430 
431  for (i = 0; i < nvariables; ++i)
432  {
433  for (n = 0; n < m_fields[i]->GetBndConditions().size(); ++n)
434  {
435  if (m_fields[i]->GetBndConditions()[n]->IsTimeDependent())
436  {
437  varName = m_session->GetVariable(i);
438  m_fields[i]->EvaluateBoundaryConditions(time, varName);
439  }
440  else if (boost::istarts_with(
441  m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
442  "Womersley"))
443  {
444  SetWomersleyBoundary(i, n);
445  }
446  }
447 
448  // Set Radiation conditions if required
450  }
451 
452  // Enforcing the boundary conditions (Inlet and wall) for the
453  // Moving reference frame
456 }
457 
458 /**
459  * Probably should be pushed back into ContField?
460  */
462 {
463  int i, n;
464 
467 
468  BndConds = m_fields[fieldid]->GetBndConditions();
469  BndExp = m_fields[fieldid]->GetBndCondExpansions();
470 
473 
474  int cnt;
475  int elmtid, nq, offset, boundary;
476  Array<OneD, NekDouble> Bvals, U;
477  int cnt1 = 0;
478 
479  for (cnt = n = 0; n < BndConds.size(); ++n)
480  {
481  std::string type = BndConds[n]->GetUserDefined();
482 
483  if ((BndConds[n]->GetBoundaryConditionType() ==
485  (boost::iequals(type, "Radiation")))
486  {
487  for (i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
488  {
489  elmtid = m_fieldsBCToElmtID[m_velocity[fieldid]][cnt];
490  elmt = m_fields[fieldid]->GetExp(elmtid);
491  offset = m_fields[fieldid]->GetPhys_Offset(elmtid);
492 
493  U = m_fields[fieldid]->UpdatePhys() + offset;
494  Bc = BndExp[n]->GetExp(i);
495 
496  boundary = m_fieldsBCToTraceID[fieldid][cnt];
497 
498  // Get edge values and put into ubc
499  nq = Bc->GetTotPoints();
500  Array<OneD, NekDouble> ubc(nq);
501  elmt->GetTracePhysVals(boundary, Bc, U, ubc);
502 
503  Vmath::Vmul(nq,
505  [fieldid][cnt1 + BndExp[n]->GetPhys_Offset(i)],
506  1, &ubc[0], 1, &ubc[0], 1);
507 
508  Bvals =
509  BndExp[n]->UpdateCoeffs() + BndExp[n]->GetCoeff_Offset(i);
510 
511  Bc->IProductWRTBase(ubc, Bvals);
512  }
513  cnt1 += BndExp[n]->GetTotPoints();
514  }
515  else
516  {
517  cnt += BndExp[n]->GetExpSize();
518  }
519  }
520 }
521 
523 {
524  // use static trip since cannot use UserDefinedTag for zero
525  // velocity and have time dependent conditions
526  static bool Setup = false;
527 
528  if (Setup == true)
529  {
530  return;
531  }
532  Setup = true;
533 
534  int i, n;
535 
537  BndConds(m_spacedim);
539 
540  for (i = 0; i < m_spacedim; ++i)
541  {
542  BndConds[i] = m_fields[m_velocity[i]]->GetBndConditions();
543  BndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
544  }
545 
547 
548  int cnt;
549  int elmtid, nq, boundary;
550 
552  Array<OneD, NekDouble> Bphys, Bcoeffs;
553 
554  int fldid = m_velocity[0];
555 
556  for (cnt = n = 0; n < BndConds[0].size(); ++n)
557  {
558  if ((BndConds[0][n]->GetBoundaryConditionType() ==
560  (boost::iequals(BndConds[0][n]->GetUserDefined(),
561  "ZeroNormalComponent")))
562  {
563  for (i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
564  {
565  elmtid = m_fieldsBCToElmtID[fldid][cnt];
566  elmt = m_fields[0]->GetExp(elmtid);
567  boundary = m_fieldsBCToTraceID[fldid][cnt];
568 
569  normals = elmt->GetTraceNormal(boundary);
570 
571  nq = BndExp[0][n]->GetExp(i)->GetTotPoints();
572  Array<OneD, NekDouble> normvel(nq, 0.0);
573 
574  for (int k = 0; k < m_spacedim; ++k)
575  {
576  Bphys = BndExp[k][n]->UpdatePhys() +
577  BndExp[k][n]->GetPhys_Offset(i);
578  Bc = BndExp[k][n]->GetExp(i);
579  Vmath::Vvtvp(nq, normals[k], 1, Bphys, 1, normvel, 1,
580  normvel, 1);
581  }
582 
583  // negate normvel for next step
584  Vmath::Neg(nq, normvel, 1);
585 
586  for (int k = 0; k < m_spacedim; ++k)
587  {
588  Bphys = BndExp[k][n]->UpdatePhys() +
589  BndExp[k][n]->GetPhys_Offset(i);
590  Bcoeffs = BndExp[k][n]->UpdateCoeffs() +
591  BndExp[k][n]->GetCoeff_Offset(i);
592  Bc = BndExp[k][n]->GetExp(i);
593  Vmath::Vvtvp(nq, normvel, 1, normals[k], 1, Bphys, 1, Bphys,
594  1);
595  Bc->FwdTransBndConstrained(Bphys, Bcoeffs);
596  }
597  }
598  }
599  else
600  {
601  cnt += BndExp[0][n]->GetExpSize();
602  }
603  }
604 }
605 
606 /**
607  * Womersley boundary condition defintion
608  */
609 void IncNavierStokes::SetWomersleyBoundary(const int fldid, const int bndid)
610 {
611  ASSERTL1(m_womersleyParams.count(bndid) == 1,
612  "Womersley parameters for this boundary have not been set up");
613 
614  WomersleyParamsSharedPtr WomParam = m_womersleyParams[fldid][bndid];
615  NekComplexDouble zvel;
616  int i, j, k;
617 
618  int M_coeffs = WomParam->m_wom_vel.size();
619 
620  NekDouble T = WomParam->m_period;
621  NekDouble axis_normal = WomParam->m_axisnormal[fldid];
622 
623  // Womersley Number
624  NekComplexDouble omega_c(2.0 * M_PI / T, 0.0);
625  NekComplexDouble k_c(0.0, 0.0);
626  NekComplexDouble m_time_c(m_time, 0.0);
627  NekComplexDouble zi(0.0, 1.0);
628  NekComplexDouble i_pow_3q2(-1.0 / sqrt(2.0), 1.0 / sqrt(2.0));
629 
631  BndCondExp = m_fields[fldid]->GetBndCondExpansions()[bndid];
632 
634  int cnt = 0;
635  int nfq;
637  int exp_npts = BndCondExp->GetExpSize();
638  Array<OneD, NekDouble> wbc(exp_npts, 0.0);
639 
640  Array<OneD, NekComplexDouble> zt(M_coeffs);
641 
642  // preallocate the exponent
643  for (k = 1; k < M_coeffs; k++)
644  {
645  k_c = NekComplexDouble((NekDouble)k, 0.0);
646  zt[k] = std::exp(zi * omega_c * k_c * m_time_c);
647  }
648 
649  // Loop over each element in an expansion
650  for (i = 0; i < exp_npts; ++i, cnt++)
651  {
652  // Get Boundary and trace expansion
653  bc = BndCondExp->GetExp(i);
654  nfq = bc->GetTotPoints();
655  Array<OneD, NekDouble> wbc(nfq, 0.0);
656 
657  // Compute womersley solution
658  for (j = 0; j < nfq; j++)
659  {
660  wbc[j] = WomParam->m_poiseuille[i][j];
661  for (k = 1; k < M_coeffs; k++)
662  {
663  zvel = WomParam->m_zvel[i][j][k] * zt[k];
664  wbc[j] = wbc[j] + zvel.real();
665  }
666  }
667 
668  // Multiply w by normal to get u,v,w component of velocity
669  Vmath::Smul(nfq, axis_normal, wbc, 1, wbc, 1);
670  // get the offset
671  Bvals = BndCondExp->UpdateCoeffs() + BndCondExp->GetCoeff_Offset(i);
672 
673  // Push back to Coeff space
674  bc->FwdTrans(wbc, Bvals);
675  }
676 }
677 
678 void IncNavierStokes::SetUpWomersley(const int fldid, const int bndid,
679  std::string womStr)
680 {
681  std::string::size_type indxBeg = womStr.find_first_of(':') + 1;
682  string filename = womStr.substr(indxBeg, string::npos);
683 
684  TiXmlDocument doc(filename);
685 
686  bool loadOkay = doc.LoadFile();
687  ASSERTL0(loadOkay,
688  (std::string("Failed to load file: ") + filename).c_str());
689 
690  TiXmlHandle docHandle(&doc);
691 
692  int err; /// Error value returned by TinyXML.
693 
694  TiXmlElement *nektar = doc.FirstChildElement("NEKTAR");
695  ASSERTL0(nektar, "Unable to find NEKTAR tag in file.");
696 
697  TiXmlElement *wombc = nektar->FirstChildElement("WOMERSLEYBC");
698  ASSERTL0(wombc, "Unable to find WOMERSLEYBC tag in file.");
699 
700  // read womersley parameters
701  TiXmlElement *womparam = wombc->FirstChildElement("WOMPARAMS");
702  ASSERTL0(womparam, "Unable to find WOMPARAMS tag in file.");
703 
704  // Input coefficients
705  TiXmlElement *params = womparam->FirstChildElement("W");
706  map<std::string, std::string> Wparams;
707 
708  // read parameter list
709  while (params)
710  {
711 
712  std::string propstr;
713  propstr = params->Attribute("PROPERTY");
714 
715  ASSERTL0(!propstr.empty(),
716  "Failed to read PROPERTY value Womersley BC Parameter");
717 
718  std::string valstr;
719  valstr = params->Attribute("VALUE");
720 
721  ASSERTL0(!valstr.empty(),
722  "Failed to read VALUE value Womersley BC Parameter");
723 
724  std::transform(propstr.begin(), propstr.end(), propstr.begin(),
725  ::toupper);
726  Wparams[propstr] = valstr;
727 
728  params = params->NextSiblingElement("W");
729  }
730  bool parseGood;
731 
732  // Read parameters
733 
734  ASSERTL0(
735  Wparams.count("RADIUS") == 1,
736  "Failed to find Radius parameter in Womersley boundary conditions");
737  std::vector<NekDouble> rad;
738  ParseUtils::GenerateVector(Wparams["RADIUS"], rad);
739  m_womersleyParams[fldid][bndid]->m_radius = rad[0];
740 
741  ASSERTL0(
742  Wparams.count("PERIOD") == 1,
743  "Failed to find period parameter in Womersley boundary conditions");
744  std::vector<NekDouble> period;
745  parseGood = ParseUtils::GenerateVector(Wparams["PERIOD"], period);
746  m_womersleyParams[fldid][bndid]->m_period = period[0];
747 
748  ASSERTL0(
749  Wparams.count("AXISNORMAL") == 1,
750  "Failed to find axisnormal parameter in Womersley boundary conditions");
751  std::vector<NekDouble> anorm;
752  parseGood = ParseUtils::GenerateVector(Wparams["AXISNORMAL"], anorm);
753  m_womersleyParams[fldid][bndid]->m_axisnormal[0] = anorm[0];
754  m_womersleyParams[fldid][bndid]->m_axisnormal[1] = anorm[1];
755  m_womersleyParams[fldid][bndid]->m_axisnormal[2] = anorm[2];
756 
757  ASSERTL0(
758  Wparams.count("AXISPOINT") == 1,
759  "Failed to find axispoint parameter in Womersley boundary conditions");
760  std::vector<NekDouble> apt;
761  parseGood = ParseUtils::GenerateVector(Wparams["AXISPOINT"], apt);
762  m_womersleyParams[fldid][bndid]->m_axispoint[0] = apt[0];
763  m_womersleyParams[fldid][bndid]->m_axispoint[1] = apt[1];
764  m_womersleyParams[fldid][bndid]->m_axispoint[2] = apt[2];
765 
766  // Read Temporal Fourier Coefficients.
767 
768  // Find the FourierCoeff tag
769  TiXmlElement *coeff = wombc->FirstChildElement("FOURIERCOEFFS");
770 
771  // Input coefficients
772  TiXmlElement *fval = coeff->FirstChildElement("F");
773 
774  int indx;
775  int nextFourierCoeff = -1;
776 
777  while (fval)
778  {
779  nextFourierCoeff++;
780 
781  TiXmlAttribute *fvalAttr = fval->FirstAttribute();
782  std::string attrName(fvalAttr->Name());
783 
784  ASSERTL0(attrName == "ID",
785  (std::string("Unknown attribute name: ") + attrName).c_str());
786 
787  err = fvalAttr->QueryIntValue(&indx);
788  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
789 
790  std::string coeffStr = fval->FirstChild()->ToText()->ValueStr();
791  vector<NekDouble> coeffvals;
792 
793  parseGood = ParseUtils::GenerateVector(coeffStr, coeffvals);
794  ASSERTL0(
795  parseGood,
796  (std::string("Problem reading value of fourier coefficient, ID=") +
797  boost::lexical_cast<string>(indx))
798  .c_str());
799  ASSERTL1(
800  coeffvals.size() == 2,
801  (std::string(
802  "Have not read two entries of Fourier coefficicent from ID=" +
803  boost::lexical_cast<string>(indx))
804  .c_str()));
805 
806  m_womersleyParams[fldid][bndid]->m_wom_vel.push_back(
807  NekComplexDouble(coeffvals[0], coeffvals[1]));
808 
809  fval = fval->NextSiblingElement("F");
810  }
811 
812  // starting point of precalculation
813  int i, j, k;
814  // M fourier coefficients
815  int M_coeffs = m_womersleyParams[fldid][bndid]->m_wom_vel.size();
816  NekDouble R = m_womersleyParams[fldid][bndid]->m_radius;
817  NekDouble T = m_womersleyParams[fldid][bndid]->m_period;
818  Array<OneD, NekDouble> x0 = m_womersleyParams[fldid][bndid]->m_axispoint;
819 
820  NekComplexDouble rqR;
821  // Womersley Number
822  NekComplexDouble omega_c(2.0 * M_PI / T, 0.0);
823  NekComplexDouble alpha_c(R * sqrt(omega_c.real() / m_kinvis), 0.0);
824  NekComplexDouble z1(1.0, 0.0);
825  NekComplexDouble i_pow_3q2(-1.0 / sqrt(2.0), 1.0 / sqrt(2.0));
826 
828  BndCondExp = m_fields[fldid]->GetBndCondExpansions()[bndid];
829 
831  int cnt = 0;
832  int nfq;
834 
835  int exp_npts = BndCondExp->GetExpSize();
836  Array<OneD, NekDouble> wbc(exp_npts, 0.0);
837 
838  // allocate time indepedent variables
839  m_womersleyParams[fldid][bndid]->m_poiseuille =
841  m_womersleyParams[fldid][bndid]->m_zvel =
843  // could use M_coeffs - 1 but need to avoid complicating things
844  Array<OneD, NekComplexDouble> zJ0(M_coeffs);
845  Array<OneD, NekComplexDouble> lamda_n(M_coeffs);
846  Array<OneD, NekComplexDouble> k_c(M_coeffs);
847  NekComplexDouble zJ0r;
848 
849  for (k = 1; k < M_coeffs; k++)
850  {
851  k_c[k] = NekComplexDouble((NekDouble)k, 0.0);
852  lamda_n[k] = i_pow_3q2 * alpha_c * sqrt(k_c[k]);
853  zJ0[k] = Polylib::ImagBesselComp(0, lamda_n[k]);
854  }
855 
856  // Loop over each element in an expansion
857  for (i = 0; i < exp_npts; ++i, cnt++)
858  {
859  // Get Boundary and trace expansion
860  bc = BndCondExp->GetExp(i);
861  nfq = bc->GetTotPoints();
862 
863  Array<OneD, NekDouble> x(nfq, 0.0);
864  Array<OneD, NekDouble> y(nfq, 0.0);
865  Array<OneD, NekDouble> z(nfq, 0.0);
866  bc->GetCoords(x, y, z);
867 
868  m_womersleyParams[fldid][bndid]->m_poiseuille[i] =
870  m_womersleyParams[fldid][bndid]->m_zvel[i] =
872 
873  // Compute coefficients
874  for (j = 0; j < nfq; j++)
875  {
876  rqR = NekComplexDouble(sqrt((x[j] - x0[0]) * (x[j] - x0[0]) +
877  (y[j] - x0[1]) * (y[j] - x0[1]) +
878  (z[j] - x0[2]) * (z[j] - x0[2])) /
879  R,
880  0.0);
881 
882  // Compute Poiseulle Flow
883  m_womersleyParams[fldid][bndid]->m_poiseuille[i][j] =
884  m_womersleyParams[fldid][bndid]->m_wom_vel[0].real() *
885  (1. - rqR.real() * rqR.real());
886 
887  m_womersleyParams[fldid][bndid]->m_zvel[i][j] =
889 
890  // compute the velocity information
891  for (k = 1; k < M_coeffs; k++)
892  {
893  zJ0r = Polylib::ImagBesselComp(0, rqR * lamda_n[k]);
894  m_womersleyParams[fldid][bndid]->m_zvel[i][j][k] =
895  m_womersleyParams[fldid][bndid]->m_wom_vel[k] *
896  (z1 - (zJ0r / zJ0[k]));
897  }
898  }
899  }
900 }
901 
902 /**
903  * Set boundary conditions for moving frame of reference
904  */
906 {
907  SetMRFWallBCs(time);
908  SetMRFDomainVelBCs(time);
909 }
910 
911 /**
912  * Set Wall boundary conditions for moving frame of reference
913  */
915 {
916  // for the wall we need to calculate:
917  // [V_wall]_xyz = [V_frame]_xyz + [Omega X r]_xyz
918  // Note all vectors must be in moving frame coordinates xyz
919  // not in inertial frame XYZ
920 
922  BndConds(m_spacedim);
924 
925  for (int i = 0; i < m_spacedim; ++i)
926  {
927  BndConds[i] = m_fields[m_velocity[i]]->GetBndConditions();
928  BndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
929  }
930 
931  int npoints;
932  int fldid = m_velocity[0];
933 
934  Array<OneD, NekDouble> Bphys, Bcoeffs;
935 
936  // loop over the boundary regions
937  for (int n = 0; n < BndExp[0].size(); ++n)
938  {
939  if (BndConds[0][n]->GetBoundaryConditionType() ==
941  (boost::iequals(BndConds[0][n]->GetUserDefined(),
942  "MovingFrameWall")))
943  {
944  npoints = BndExp[0][n]->GetNpoints();
945  //
946  // define arrays to calculate the prescribed velocities and modifed
947  // ones
951  for (int k = 0; k < m_velocity.size(); ++k)
952  {
953  velocities[k] = Array<OneD, NekDouble>(npoints, 0.0);
954  unitArray[k] = Array<OneD, NekDouble>(npoints, 1.0);
955  }
956  Array<OneD, NekDouble> tmp(npoints, 0.0);
957  for (int k = 0; k < 3; ++k)
958  {
959  coords[k] = Array<OneD, NekDouble>(npoints, 0.0);
960  }
961  BndExp[0][n]->GetCoords(coords[0], coords[1], coords[2]);
962 
963  // move the centre to the location of pivot
964  for (int i = 0; i < m_spacedim; ++i)
965  {
966  Vmath::Sadd(npoints, -m_pivotPoint[i], coords[i], 1, coords[i],
967  1);
968  }
969  /////////////////////////////////////
970  // Note that both Omega and the absolute velocities have been
971  // expressed in the moving reference frame unit vectors
972  // therefore the result is in moving ref frame and no furhter
973  // transformaton is required
974  //
975  // compute Omega X r = vx ex + vy ey + vz ez
976  // Note OmegaX : movingFrameVelsxyz[3]
977  // Note OmegaY : movingFrameVelsxyz[4]
978  // Note OmegaZ : movingFrameVelsxyz[5]
979  //
980  // vx = OmegaY*z-OmegaZ*y
981  Vmath::Smul(npoints, -1 * m_movingFrameVelsxyz[5], coords[1], 1,
982  velocities[0], 1);
983  // vy = OmegaZ*x-OmegaX*z
984  Vmath::Smul(npoints, m_movingFrameVelsxyz[5], coords[0], 1,
985  velocities[1], 1);
986  if (m_spacedim == 3)
987  {
988  // add the OmegaY*z to vx
989  Vmath::Svtvp(npoints, m_movingFrameVelsxyz[4], coords[2], 1,
990  velocities[0], 1, velocities[0], 1);
991  // add the -OmegaX*z to vy
992  Vmath::Svtvp(npoints, -1 * m_movingFrameVelsxyz[3], coords[2],
993  1, velocities[1], 1, velocities[1], 1);
994 
995  // vz = OmegaX*y-OmegaY*x
996  Vmath::Svtsvtp(npoints, m_movingFrameVelsxyz[3], coords[1], 1,
997  -1.0 * m_movingFrameVelsxyz[4], coords[0], 1,
998  velocities[2], 1);
999  }
1000 
1001  // add the translation velocity
1002  for (int k = 0; k < m_spacedim; ++k)
1003  {
1004  Vmath::Sadd(npoints, m_movingFrameVelsxyz[k], velocities[k], 1,
1005  velocities[k], 1);
1006  }
1007 
1008  // update the boundary values
1009  for (int k = 0; k < m_spacedim; ++k)
1010  {
1011  Vmath::Vmul(npoints, unitArray[k], 1, velocities[k], 1,
1012  BndExp[k][n]->UpdatePhys(), 1);
1013  }
1014  // update the coefficients
1015  for (int k = 0; k < m_spacedim; ++k)
1016  {
1017  if (m_fields[k]->GetExpType() == MultiRegions::e3DH1D)
1018  {
1019  BndExp[k][n]->SetWaveSpace(false);
1020  }
1021  BndExp[k][n]->FwdTransBndConstrained(
1022  BndExp[k][n]->GetPhys(), BndExp[k][n]->UpdateCoeffs());
1023  }
1024  }
1025  }
1026 }
1027 
1028 /**
1029  * Set inlet boundary conditions for moving frame of reference
1030  */
1032 {
1033  // The inlet conditions for the velocities given in the session xml file
1034  // however, those are in the inertial reference coordinate XYZ and
1035  // needed to be converted to the moving reference coordinate xyz
1036 
1037  // define arrays to calculate the prescribed velocities and modifed ones
1038  Array<OneD, Array<OneD, NekDouble>> definedVels(m_velocity.size());
1039  Array<OneD, Array<OneD, NekDouble>> velocities(m_velocity.size());
1042 
1044  BndConds(m_spacedim);
1046 
1047  for (int i = 0; i < m_spacedim; ++i)
1048  {
1049  BndConds[i] = m_fields[m_velocity[i]]->GetBndConditions();
1050  BndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
1051  }
1052 
1053  int npoints;
1054  int fldid = m_velocity[0];
1055 
1056  Array<OneD, NekDouble> Bphys, Bcoeffs;
1057 
1058  // loop over the boundary regions
1059  for (int n = 0; n < BndExp[0].size(); ++n)
1060  {
1061  if (BndConds[0][n]->GetBoundaryConditionType() ==
1063  (boost::iequals(BndConds[0][n]->GetUserDefined(),
1064  "MovingFrameDomainVel")))
1065  {
1066  npoints = BndExp[0][n]->GetNpoints();
1067  for (int k = 0; k < m_velocity.size(); ++k)
1068  {
1069  definedVels[k] = Array<OneD, NekDouble>(npoints, 0.0);
1070  velocities[k] = Array<OneD, NekDouble>(npoints, 0.0);
1071  unitArray[k] = Array<OneD, NekDouble>(npoints, 1.0);
1072  }
1073  Array<OneD, NekDouble> tmp(npoints, 0.0);
1074  for (int k = 0; k < 3; ++k)
1075  {
1076  coords[k] = Array<OneD, NekDouble>(npoints, 0.0);
1077  }
1078  BndExp[0][n]->GetCoords(coords[0], coords[1], coords[2]);
1079 
1080  // loop over the velocity fields and compute the boundary
1081  // condition
1082  for (int k = 0; k < m_velocity.size(); ++k)
1083  {
1084  LibUtilities::Equation condition =
1085  std::static_pointer_cast<
1087  BndConds[k][n])
1088  ->m_dirichletCondition;
1089  // Evaluate
1090  condition.Evaluate(coords[0], coords[1], coords[2], time,
1091  definedVels[k]);
1092  }
1093 
1094  // We have all velocity components
1095  // transform them to the moving refernce frame
1096  for (int l = 0; l < m_spacedim; ++l)
1097  {
1098  Array<OneD, NekDouble> tmp0(npoints, 0.0);
1099  Array<OneD, NekDouble> tmp1(npoints, 0.0);
1100  Array<OneD, NekDouble> tmp2(npoints, 0.0);
1101  for (int m = 0; m < m_spacedim; ++m)
1102  {
1103  Vmath::Svtvp(npoints, m_movingFrameProjMat(l, m),
1104  tmp0 = definedVels[m], 1, tmp1 = velocities[l],
1105  1, tmp2 = velocities[l], 1);
1106  }
1107  }
1108 
1109  // update the boundary values
1110  for (int k = 0; k < m_spacedim; ++k)
1111  {
1112  Vmath::Vmul(npoints, unitArray[k], 1, velocities[k], 1,
1113  BndExp[k][n]->UpdatePhys(), 1);
1114  }
1115  // update the coefficients
1116  for (int k = 0; k < m_spacedim; ++k)
1117  {
1118  if (m_fields[k]->GetExpType() == MultiRegions::e3DH1D)
1119  {
1120  BndExp[k][n]->SetWaveSpace(false);
1121  }
1122  BndExp[k][n]->FwdTransBndConstrained(
1123  BndExp[k][n]->GetPhys(), BndExp[k][n]->UpdateCoeffs());
1124  }
1125  }
1126  }
1127 }
1128 
1129 /**
1130  * Add an additional forcing term programmatically.
1131  */
1133 {
1134  m_forcing.push_back(pForce);
1135 }
1136 
1137 /**
1138  *
1139  */
1141  const NekDouble SpeedSoundFactor)
1142 {
1143  boost::ignore_unused(SpeedSoundFactor);
1144  int nvel = m_velocity.size();
1145  int nelmt = m_fields[0]->GetExpSize();
1146 
1147  Array<OneD, NekDouble> stdVelocity(nelmt, 0.0);
1149 
1150  if (m_HomogeneousType == eHomogeneous1D) // just do check on 2D info
1151  {
1152  velfields = Array<OneD, Array<OneD, NekDouble>>(2);
1153 
1154  for (int i = 0; i < 2; ++i)
1155  {
1156  velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
1157  }
1158  }
1159  else
1160  {
1161  velfields = Array<OneD, Array<OneD, NekDouble>>(nvel);
1162 
1163  for (int i = 0; i < nvel; ++i)
1164  {
1165  velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
1166  }
1167  }
1168 
1169  stdVelocity = m_extrapolation->GetMaxStdVelocity(velfields);
1170 
1171  return stdVelocity;
1172 }
1173 
1174 /**
1175  *
1176  */
1178  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1180 {
1181  pressure = physfield[m_nConvectiveFields];
1182 }
1183 
1184 /**
1185  *
1186  */
1188  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1189  Array<OneD, NekDouble> &density)
1190 {
1191  int nPts = physfield[0].size();
1192  Vmath::Fill(nPts, 1.0, density, 1);
1193 }
1194 
1195 /**
1196  *
1197  */
1199  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1200  Array<OneD, Array<OneD, NekDouble>> &velocity)
1201 {
1202  for (int i = 0; i < m_spacedim; ++i)
1203  {
1204  velocity[i] = physfield[i];
1205  }
1206 }
1207 
1208 /**
1209  * Function to set the moving frame velocities calucated in the forcing
1210  * this gives access to the moving reference forcing to set the velocities
1211  * to be later used in enforcing the boundary condition in IncNavierStokes
1212  * class
1213  */
1215  const Array<OneD, NekDouble> &vFrameVels)
1216 {
1217  ASSERTL0(vFrameVels.size() == m_movingFrameVelsxyz.size(),
1218  "Arrays have different dimensions, cannot set moving frame "
1219  "velocities");
1220  Vmath::Vcopy(vFrameVels.size(), vFrameVels, 1, m_movingFrameVelsxyz, 1);
1221 }
1222 
1224  Array<OneD, NekDouble> &vFrameVels)
1225 {
1226  ASSERTL0(vFrameVels.size() == m_movingFrameVelsxyz.size(),
1227  "Arrays have different dimensions, cannot get moving frame "
1228  "velocities");
1229  unsigned int size = m_movingFrameVelsxyz.size();
1230  Vmath::Vcopy(size, m_movingFrameVelsxyz, 1, vFrameVels, 1);
1231 }
1232 
1233 /**
1234  * Function to set the rotation matrix computed in the forcing
1235  * this gives access to the moving reference forcing to set the Projection
1236  * matrix to be used later in IncNavierStokes calss for enforcing the
1237  * boundary conditions
1238  */
1240  const bnu::matrix<NekDouble> &vProjMat)
1241 {
1242  ASSERTL0(vProjMat.size1() == m_movingFrameProjMat.size1(),
1243  "Matrices have different numbers of rows, cannot Set the "
1244  "moving frame projection matrix");
1245  ASSERTL0(vProjMat.size2() == m_movingFrameProjMat.size2(),
1246  "Matrices have different numbers of columns, cannot Set the "
1247  "moving frame projection matrix");
1248  for (int i = 0; i < vProjMat.size1(); ++i)
1249  {
1250  for (int j = 0; j < vProjMat.size2(); ++j)
1251  {
1252  m_movingFrameProjMat(i, j) = vProjMat(i, j);
1253  }
1254  }
1255 }
1256 
1258  bnu::matrix<NekDouble> &vProjMat)
1259 {
1260  ASSERTL0(vProjMat.size1() == m_movingFrameProjMat.size1(),
1261  "Matrices have different numbers of rows, cannot Get the "
1262  "moving frame projection matrix");
1263  ASSERTL0(vProjMat.size2() == m_movingFrameProjMat.size2(),
1264  "Matrices have different numbers of columns, cannot Get the "
1265  "moving frame projection matrix");
1266 
1267  for (int i = 0; i < vProjMat.size1(); ++i)
1268  {
1269  for (int j = 0; j < vProjMat.size2(); ++j)
1270  {
1271  vProjMat(i, j) = m_movingFrameProjMat(i, j);
1272  }
1273  }
1274 }
1275 
1276 /**
1277  * Function to set the angles between the moving frame of reference and
1278  * stationary inertial reference frame
1279  **/
1281  const Array<OneD, NekDouble> &vFrameTheta)
1282 {
1283  ASSERTL0(vFrameTheta.size() == m_movingFrameTheta.size(),
1284  "Arrays have different size, cannot set moving frame angles");
1285  for (int i = 0; i < vFrameTheta.size(); ++i)
1286  {
1287  m_movingFrameTheta[i] = vFrameTheta[i];
1288  }
1289 }
1290 
1291 /**
1292  * Function to get the angles between the moving frame of reference and
1293  * stationary inertial reference frame
1294  **/
1296 {
1297  ASSERTL0(vFrameTheta.size() == m_movingFrameTheta.size(),
1298  "Arrays have different size, cannot get moving frame angles");
1299  for (int i = 0; i < m_movingFrameTheta.size(); ++i)
1300  {
1301  vFrameTheta[i] = m_movingFrameTheta[i];
1302  }
1303 }
1304 
1305 /**
1306  * Function to check the type of forcing
1307  **/
1308 bool IncNavierStokes::DefinedForcing(const std::string &sForce)
1309 {
1310  vector<std::string> vForceList;
1311  bool hasForce{false};
1312 
1313  if (!m_session->DefinesElement("Nektar/Forcing"))
1314  {
1315  return hasForce;
1316  }
1317 
1318  TiXmlElement *vForcing = m_session->GetElement("Nektar/Forcing");
1319  if (vForcing)
1320  {
1321  TiXmlElement *vForce = vForcing->FirstChildElement("FORCE");
1322  while (vForce)
1323  {
1324  string vType = vForce->Attribute("TYPE");
1325 
1326  vForceList.push_back(vType);
1327  vForce = vForce->NextSiblingElement("FORCE");
1328  }
1329  }
1330 
1331  for (auto &f : vForceList)
1332  {
1333  if (boost::iequals(f, sForce))
1334  {
1335  hasForce = true;
1336  }
1337  }
1338 
1339  return hasForce;
1340 }
1341 
1342 /**
1343  * Get the pivot point for moving reference
1344  **/
1346 {
1347  std::string sMRFForcingType = "MovingReferenceFrame";
1348 
1349  // only if we use moving reference frame formulation
1350  if (!DefinedForcing(sMRFForcingType))
1351  {
1352  return;
1353  }
1354 
1355  TiXmlElement *vMRFForcing;
1356 
1357  TiXmlElement *vForcing = m_session->GetElement("Nektar/Forcing");
1358  if (vForcing)
1359  {
1360  TiXmlElement *vForce = vForcing->FirstChildElement("FORCE");
1361  while (vForce)
1362  {
1363  string vType = vForce->Attribute("TYPE");
1364 
1365  // if it is moving reference frame
1366  if (boost::iequals(vType, sMRFForcingType))
1367  {
1368  TiXmlElement *pivotElmt =
1369  vForce->FirstChildElement("PivotPoint");
1370 
1371  // if not defined, zero would be the default
1372  if (pivotElmt)
1373  {
1374  std::stringstream pivotPointStream;
1375  std::string pivotPointStr = pivotElmt->GetText();
1376  pivotPointStream.str(pivotPointStr);
1377 
1378  for (int j = 0; j < m_spacedim; ++j)
1379  {
1380  pivotPointStream >> pivotPointStr;
1381  if (!pivotPointStr.empty())
1382  {
1384  m_session->GetInterpreter(), pivotPointStr);
1385  vPivotPoint[j] = equ.Evaluate();
1386  }
1387  }
1388  }
1389 
1390  break;
1391  }
1392  }
1393  }
1394 }
1395 
1396 /**
1397  * Perform the extrapolation.
1398  */
1400 {
1401  m_extrapolation->SubStepSaveFields(step);
1402  m_extrapolation->SubStepAdvance(step, m_time);
1404  return false;
1405 }
1406 
1407 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
Array< OneD, int > & GetVelocity(void)
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
Mapping from BCs to Elmt Edge IDs.
Array< OneD, NekDouble > m_pivotPoint
virtual Array< OneD, NekDouble > v_GetMaxStdVelocity(const NekDouble SpeedSoundFactor) override
virtual void SetMovingFrameAngles(const Array< OneD, NekDouble > &vFrameTheta) override
std::map< int, std::map< int, WomersleyParamsSharedPtr > > m_womersleyParams
Womersley parameters if required.
virtual void GetMovingFrameAngles(Array< OneD, NekDouble > &vFrameTheta) override
virtual int v_GetForceDimension()=0
void SetWomersleyBoundary(const int fldid, const int bndid)
Set Womersley Profile if specified.
void SetZeroNormalVelocity()
Set Normal Velocity Component to Zero.
virtual void SetMovingFrameVelocities(const Array< OneD, NekDouble > &vFrameVels) override
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
virtual void GetMovingFrameVelocities(Array< OneD, NekDouble > &vFrameVels) override
virtual void GetMovingFrameProjectionMat(bnu::matrix< NekDouble > &vProjMat) override
void SetMRFDomainVelBCs(const NekDouble &time)
NekDouble m_kinvis
Kinematic viscosity.
void SetRadiationBoundaryForcing(int fieldid)
Set Radiation forcing term.
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
RHS Factor for Radiation Condition.
virtual bool DefinedForcing(const std::string &sForce)
void SetMRFWallBCs(const NekDouble &time)
virtual void GetPivotPoint(Array< OneD, NekDouble > &vPivotPoint)
virtual bool v_PreIntegrate(int step) override
void SetUpWomersley(const int fldid, const int bndid, std::string womstr)
Set Up Womersley details.
virtual void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
virtual void GetDensity(const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, NekDouble > &density) override
Extract array with density from physfield.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
Mapping from BCs to Elmt IDs.
EquationType m_equationType
equation type;
int m_nConvectiveFields
Number of fields to be convected;.
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
virtual void SetMovingFrameProjectionMat(const bnu::matrix< NekDouble > &vProjMat) override
void SetMovingReferenceFrameBCs(const NekDouble &time)
Set the moving reference frame boundary conditions.
NekDouble Evaluate() const
Definition: Equation.cpp:95
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
Definition: ParseUtils.cpp:129
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
int m_spacedim
Spatial dimension (>= expansion dim).
boost::numeric::ublas::matrix< NekDouble > m_movingFrameProjMat
Projection matrix for transformation between inertial and moving.
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Array< OneD, NekDouble > m_movingFrameVelsxyz
Moving frame of reference velocities.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Array< OneD, NekDouble > m_movingFrameTheta
Moving frame of reference angles with respect to the.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:120
Base class for unsteady solvers.
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter.cpp:86
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:53
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47
SOLVER_UTILS_EXPORT typedef std::shared_ptr< Forcing > ForcingSharedPtr
A shared pointer to an EquationSystem object.
Definition: Forcing.h:52
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::complex< double > NekComplexDouble
@ eSteadyNavierStokes
@ eUnsteadyStokes
@ eUnsteadyNavierStokes
@ eSteadyLinearisedNS
@ eUnsteadyLinearisedNS
@ eEquationTypeSize
std::shared_ptr< WomersleyParams > WomersleyParamsSharedPtr
const std::string kEquationTypeStr[]
double NekDouble
std::complex< Nektar::NekDouble > ImagBesselComp(int n, std::complex< Nektar::NekDouble > y)
Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Rec...
Definition: Polylib.cpp:1731
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:751
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:209
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291