Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Incompressible Navier Stokes class definition built on
33 // ADRBase class
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iomanip>
38 #include <boost/algorithm/string.hpp>
39 
46 
50 #include <algorithm>
51 #include <complex>
52 #include <iostream>
53 #include <fstream>
54 #include <sstream>
55 
56 #include <tinyxml.h>
58 
59 using namespace std;
60 
61 namespace Nektar
62 {
63 
64  /**
65  * Constructor. Creates ...
66  *
67  * \param
68  * \param
69  */
70  IncNavierStokes::IncNavierStokes(const LibUtilities::SessionReaderSharedPtr& pSession):
71  UnsteadySystem(pSession),
72  AdvectionSystem(pSession),
73  m_SmoothAdvection(false),
74  m_steadyStateSteps(0)
75  {
76  }
77 
79  {
81 
82  int i,j;
83  int numfields = m_fields.num_elements();
84  std::string velids[] = {"u","v","w"};
85 
86  // Set up Velocity field to point to the first m_expdim of m_fields;
88 
89  for(i = 0; i < m_spacedim; ++i)
90  {
91  for(j = 0; j < numfields; ++j)
92  {
93  std::string var = m_boundaryConditions->GetVariable(j);
94  if(boost::iequals(velids[i], var))
95  {
96  m_velocity[i] = j;
97  break;
98  }
99 
100  ASSERTL0(j != numfields, "Failed to find field: " + var);
101  }
102  }
103 
104  // Set up equation type enum using kEquationTypeStr
105  for(i = 0; i < (int) eEquationTypeSize; ++i)
106  {
107  bool match;
108  m_session->MatchSolverInfo("EQTYPE",kEquationTypeStr[i],match,false);
109  if(match)
110  {
112  break;
113  }
114  }
115  ASSERTL0(i != eEquationTypeSize,"EQTYPE not found in SOLVERINFO section");
116 
117  // This probably should to into specific implementations
118  // Equation specific Setups
119  switch(m_equationType)
120  {
121  case eSteadyStokes:
122  case eSteadyOseen:
123  case eSteadyNavierStokes:
124  case eSteadyLinearisedNS:
125  break;
127  case eUnsteadyStokes:
128  {
129  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
130  m_session->LoadParameter("IO_CFLSteps", m_cflsteps, 0);
131  m_session->LoadParameter("SteadyStateSteps", m_steadyStateSteps, 0);
132  m_session->LoadParameter("SteadyStateTol", m_steadyStateTol, 1e-6);
133  }
134  break;
135  case eNoEquationType:
136  default:
137  ASSERTL0(false,"Unknown or undefined equation type");
138  }
139 
140  m_session->LoadParameter("Kinvis", m_kinvis);
141 
142  // Default advection type per solver
143  std::string vConvectiveType;
144  switch(m_equationType)
145  {
146  case eUnsteadyStokes:
147  case eSteadyLinearisedNS:
148  vConvectiveType = "NoAdvection";
149  break;
151  case eSteadyNavierStokes:
152  vConvectiveType = "Convective";
153  break;
155  vConvectiveType = "Linearised";
156  break;
157  default:
158  break;
159  }
160 
161  // Check if advection type overridden
162  if (m_session->DefinesTag("AdvectiveType") && m_equationType != eUnsteadyStokes &&
164  {
165  vConvectiveType = m_session->GetTag("AdvectiveType");
166  }
167 
168  // Initialise advection
169  m_advObject = SolverUtils::GetAdvectionFactory().CreateInstance(vConvectiveType, vConvectiveType);
170  m_advObject->InitObject( m_session, m_fields);
171 
172  // Forcing terms
175 
176  // check to see if any Robin boundary conditions and if so set
177  // up m_field to boundary condition maps;
181 
182  for (i = 0; i < m_fields.num_elements(); ++i)
183  {
184  bool Set = false;
185 
188  int radpts = 0;
189 
190  BndConds = m_fields[i]->GetBndConditions();
191  BndExp = m_fields[i]->GetBndCondExpansions();
192  for(int n = 0; n < BndConds.num_elements(); ++n)
193  {
194  if(boost::iequals(BndConds[n]->GetUserDefined(),"Radiation"))
195  {
196  ASSERTL0(BndConds[n]->GetBoundaryConditionType() == SpatialDomains::eRobin,
197  "Radiation boundary condition must be of type Robin <R>");
198 
199  if(Set == false)
200  {
201  m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],m_fieldsBCToTraceID[i]);
202  Set = true;
203  }
204  radpts += BndExp[n]->GetTotPoints();
205  }
206  if(boost::iequals(BndConds[n]->GetUserDefined(),"ZeroNormalComponent"))
207  {
208  ASSERTL0(BndConds[n]->GetBoundaryConditionType() == SpatialDomains::eDirichlet,
209  "Zero Normal Component boundary condition option must be of type Dirichlet <D>");
210 
211  if(Set == false)
212  {
213  m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],m_fieldsBCToTraceID[i]);
214  Set = true;
215  }
216  }
217  }
218 
219  m_fieldsRadiationFactor[i] = Array<OneD, NekDouble>(radpts);
220 
221  radpts = 0; // reset to use as a counter
222 
223  for(int n = 0; n < BndConds.num_elements(); ++n)
224  {
225  if(boost::iequals(BndConds[n]->GetUserDefined(),"Radiation"))
226  {
227 
228  int npoints = BndExp[n]->GetNpoints();
229  Array<OneD, NekDouble> x0(npoints,0.0);
230  Array<OneD, NekDouble> x1(npoints,0.0);
231  Array<OneD, NekDouble> x2(npoints,0.0);
232  Array<OneD, NekDouble> tmpArray;
233 
234  BndExp[n]->GetCoords(x0,x1,x2);
235 
236  LibUtilities::Equation coeff =
237  boost::static_pointer_cast<
239  >(BndConds[n])->m_robinPrimitiveCoeff;
240 
241  coeff.Evaluate(x0,x1,x2,m_time,
242  tmpArray = m_fieldsRadiationFactor[i]+ radpts);
243  //Vmath::Neg(npoints,tmpArray = m_fieldsRadiationFactor[i]+ radpts,1);
244  radpts += npoints;
245  }
246  }
247  }
248 
249  // Set up maping for womersley BC - and load variables
250  for (int i = 0; i < m_fields.num_elements(); ++i)
251  {
252  for(int n = 0; n < m_fields[i]->GetBndConditions().num_elements(); ++n)
253  {
254  if(boost::istarts_with(m_fields[i]->GetBndConditions()[n]->GetUserDefined(),"Womersley"))
255  {
256 
258 
259 
260 #if 0
261  m_session->LoadParameter("Period",m_womersleyParams[n]->m_period);
262  m_session->LoadParameter("Radius",m_womersleyParams[n]->m_radius);
263 
264  NekDouble n0,n1,n2;
265  m_session->LoadParameter("n0",n0);
266  m_session->LoadParameter("n1",n1);
267  m_session->LoadParameter("n2",n2);
268  m_womersleyParams[n]->m_axisnormal[0] = n0;
269  m_womersleyParams[n]->m_axisnormal[1] = n1;
270  m_womersleyParams[n]->m_axisnormal[2] = n2;
271 
272  NekDouble x0,x1,x2;
273  m_session->LoadParameter("x0",x0);
274  m_session->LoadParameter("x1",x1);
275  m_session->LoadParameter("x2",x2);
276  m_womersleyParams[n]->m_axispoint[0] = x0;
277  m_womersleyParams[n]->m_axispoint[1] = x1;
278  m_womersleyParams[n]->m_axispoint[2] = x2;
279 #endif
280 
281  // Read in fourier coeffs
282  SetUpWomersley(n,
283  m_fields[i]->GetBndConditions()[n]->GetUserDefined());
284 
285  m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],m_fieldsBCToTraceID[i]);
286 
287  }
288  }
289  }
290 
291  // Set up Field Meta Data for output files
292  m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
293  m_fieldMetaDataMap["TimeStep"] = boost::lexical_cast<std::string>(m_timestep);
294  }
295 
296  /**
297  * Destructor
298  */
300  {
301  }
302 
303 
304  /**
305  *
306  */
308  Array<OneD, Array<OneD, NekDouble> > &physfield,
310  {
311  ASSERTL1(flux.num_elements() == m_velocity.num_elements(),"Dimension of flux array and velocity array do not match");
312 
313  for(int j = 0; j < flux.num_elements(); ++j)
314  {
315  Vmath::Vmul(GetNpoints(), physfield[i], 1, m_fields[m_velocity[j]]->GetPhys(), 1, flux[j], 1);
316  }
317  }
318 
319  /**
320  * Calcualate numerical fluxes
321  */
323  Array<OneD, Array<OneD, NekDouble> > &numflux)
324  {
325  /// Counter variable
326  int i;
327 
328  /// Number of trace points
329  int nTracePts = GetTraceNpoints();
330 
331  /// Number of spatial dimensions
332  int nDimensions = m_spacedim;
333 
334  /// Forward state array
335  Array<OneD, NekDouble> Fwd(2*nTracePts);
336 
337  /// Backward state array
338  Array<OneD, NekDouble> Bwd = Fwd + nTracePts;
339 
340  /// Normal velocity array
341  Array<OneD, NekDouble> Vn (nTracePts, 0.0);
342 
343  // Extract velocity field along the trace space and multiply by trace normals
344  for(i = 0; i < nDimensions; ++i)
345  {
346  m_fields[0]->ExtractTracePhys(m_fields[m_velocity[i]]->GetPhys(), Fwd);
347  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Fwd, 1, Vn, 1, Vn, 1);
348  }
349 
350  /// Compute the numerical fluxes at the trace points
351  for(i = 0; i < numflux.num_elements(); ++i)
352  {
353  /// Extract forwards/backwards trace spaces
354  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
355 
356  /// Upwind between elements
357  m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
358 
359  /// Calculate the numerical fluxes multipling Fwd or Bwd
360  /// by the normal advection velocity
361  Vmath::Vmul(nTracePts, numflux[i], 1, Vn, 1, numflux[i], 1);
362  }
363  }
364 
365  /**
366  * Evaluation -N(V) for all fields except pressure using m_velocity
367  */
369  Array<OneD, Array<OneD, NekDouble> > &outarray)
370  {
371  int i;
372  int VelDim = m_velocity.num_elements();
373  Array<OneD, Array<OneD, NekDouble> > velocity(VelDim);
374 
375  for(i = 0; i < VelDim; ++i)
376  {
377  velocity[i] = inarray[m_velocity[i]];
378  }
379 
381  velocity, inarray, outarray, m_time);
382  }
383 
384  /**
385  * Time dependent boundary conditions updating
386  */
388  {
389  int i, n;
390  std::string varName;
391  int nvariables = m_fields.num_elements();
392 
393  for (i = 0; i < nvariables; ++i)
394  {
395  for(n = 0; n < m_fields[i]->GetBndConditions().num_elements(); ++n)
396  {
397  if(m_fields[i]->GetBndConditions()[n]->IsTimeDependent())
398  {
399  varName = m_session->GetVariable(i);
400  m_fields[i]->EvaluateBoundaryConditions(time, varName);
401  }
402  else if(boost::istarts_with(m_fields[i]->GetBndConditions()[n]->GetUserDefined(),"Womersley"))
403  {
405  }
406  }
407 
408  // Set Radiation conditions if required
410  }
411 
413  }
414 
415  /**
416  * Probably should be pushed back into ContField?
417  */
419  {
420  int i,n;
421 
424 
425  BndConds = m_fields[fieldid]->GetBndConditions();
426  BndExp = m_fields[fieldid]->GetBndCondExpansions();
427 
430 
431  int cnt;
432  int elmtid,nq,offset, boundary;
433  Array<OneD, NekDouble> Bvals, U;
434  int cnt1 = 0;
435 
436  for(cnt = n = 0; n < BndConds.num_elements(); ++n)
437  {
438  std::string type = BndConds[n]->GetUserDefined();
439 
440  if((BndConds[n]->GetBoundaryConditionType() == SpatialDomains::eRobin)&&(boost::iequals(type,"Radiation")))
441  {
442  for(i = 0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
443  {
444  elmtid = m_fieldsBCToElmtID[m_velocity[fieldid]][cnt];
445  elmt = m_fields[fieldid]->GetExp(elmtid);
446  offset = m_fields[fieldid]->GetPhys_Offset(elmtid);
447 
448  U = m_fields[fieldid]->UpdatePhys() + offset;
449  Bc = BndExp[n]->GetExp(i);
450 
451  boundary = m_fieldsBCToTraceID[fieldid][cnt];
452 
453  // Get edge values and put into ubc
454  nq = Bc->GetTotPoints();
455  Array<OneD, NekDouble> ubc(nq);
456  elmt->GetTracePhysVals(boundary,Bc,U,ubc);
457 
458  Vmath::Vmul(nq,&m_fieldsRadiationFactor[fieldid][cnt1 +
459  BndExp[n]->GetPhys_Offset(i)],1,&ubc[0],1,&ubc[0],1);
460 
461  Bvals = BndExp[n]->UpdateCoeffs()+BndExp[n]->GetCoeff_Offset(i);
462 
463  Bc->IProductWRTBase(ubc,Bvals);
464  }
465  cnt1 += BndExp[n]->GetTotPoints();
466  }
467  else
468  {
469  cnt += BndExp[n]->GetExpSize();
470  }
471  }
472  }
473 
474 
476  {
477  // use static trip since cannot use UserDefinedTag for zero
478  // velocity and have time dependent conditions
479  static bool Setup = false;
480 
481  if(Setup == true)
482  {
483  return;
484  }
485  Setup = true;
486 
487  int i,n;
488 
490  BndConds(m_spacedim);
492  BndExp(m_spacedim);
493 
494 
495  for(i = 0; i < m_spacedim; ++i)
496  {
497  BndConds[i] = m_fields[m_velocity[i]]->GetBndConditions();
498  BndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
499  }
500 
502 
503  int cnt;
504  int elmtid,nq, boundary;
505 
507  Array<OneD, NekDouble> Bphys,Bcoeffs;
508 
509  int fldid = m_velocity[0];
510 
511  for(cnt = n = 0; n < BndConds[0].num_elements(); ++n)
512  {
513  if((BndConds[0][n]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)&& (boost::iequals(BndConds[0][n]->GetUserDefined(),"ZeroNormalComponent")))
514  {
515  for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i,cnt++)
516  {
517  elmtid = m_fieldsBCToElmtID[fldid][cnt];
518  elmt = m_fields[0]->GetExp(elmtid);
519  boundary = m_fieldsBCToTraceID[fldid][cnt];
520 
521  normals = elmt->GetSurfaceNormal(boundary);
522 
523  nq = BndExp[0][n]->GetExp(i)->GetTotPoints();
524  Array<OneD, NekDouble> normvel(nq,0.0);
525 
526  for(int k = 0; k < m_spacedim; ++k)
527  {
528  Bphys = BndExp[k][n]->UpdatePhys()+
529  BndExp[k][n]->GetPhys_Offset(i);
530  Bc = BndExp[k][n]->GetExp(i);
531  Vmath::Vvtvp(nq,normals[k],1,Bphys,1,normvel,1,
532  normvel,1);
533  }
534 
535  // negate normvel for next step
536  Vmath::Neg(nq,normvel,1);
537 
538  for(int k = 0; k < m_spacedim; ++k)
539  {
540  Bphys = BndExp[k][n]->UpdatePhys()+
541  BndExp[k][n]->GetPhys_Offset(i);
542  Bcoeffs = BndExp[k][n]->UpdateCoeffs()+
543  BndExp[k][n]->GetCoeff_Offset(i);
544  Bc = BndExp[k][n]->GetExp(i);
545  Vmath::Vvtvp(nq,normvel,1,normals[k],1,Bphys,1,
546  Bphys,1);
547  Bc->FwdTrans_BndConstrained(Bphys,Bcoeffs);
548  }
549  }
550  }
551  else
552  {
553  cnt += BndExp[0][n]->GetExpSize();
554  }
555  }
556  }
557 
558 
559  /**
560  * Womersley boundary condition defintion
561  */
562  void IncNavierStokes::SetWomersleyBoundary(const int fldid, const int bndid)
563  {
564  ASSERTL1(m_womersleyParams.count(bndid) == 1, "Womersley parameters for this boundary have not been set up");
565 
567  std::complex<NekDouble> za, zar, zJ0, zJ0r, zq, zvel, zJ0rJ0;
568  int i,j,k;
569 
570  int M = WomParam->m_wom_vel_r.size();
571 
572  NekDouble R = WomParam->m_radius;
573  NekDouble T = WomParam->m_period;
574 
575  Array<OneD, NekDouble > normals = WomParam->m_axisnormal;
576  Array<OneD, NekDouble > x0 = WomParam->m_axispoint;
577 
578  // Womersley Number
579  NekDouble alpha = R*sqrt(2*M_PI/T/m_kinvis);
580  NekDouble r,kt;
581 
582  std::complex<NekDouble> z1 (1.0,0.0);
583  std::complex<NekDouble> zi (0.0,1.0);
584  std::complex<NekDouble> comp_conj (-1.0,1.0); //complex conjugate
585 
587 
588  BndExp = m_fields[fldid]->GetBndCondExpansions();
589 
592  int cnt=0;
593  int elmtid,offset, boundary,nfq;
594 
595  Array<OneD, NekDouble> Bvals,w;
596 
597  //Loop over all expansions
598  for(i = 0; i < BndExp[bndid]->GetExpSize(); ++i,cnt++)
599  {
600  // Get element id and offset
601  elmtid = m_fieldsBCToElmtID[fldid][cnt];
602  elmt = m_fields[fldid]->GetExp(elmtid);
603  offset = m_fields[fldid]->GetPhys_Offset(elmtid);
604 
605  // Get Boundary and trace expansion
606  bc = BndExp[bndid]->GetExp(i);
607  boundary = m_fieldsBCToTraceID[fldid][cnt];
608 
609  nfq=bc->GetTotPoints();
610  w = m_fields[fldid]->UpdatePhys() + offset;
611 
612  Array<OneD, NekDouble> x(nfq,0.0);
613  Array<OneD, NekDouble> y(nfq,0.0);
614  Array<OneD, NekDouble> z(nfq,0.0);
615  Array<OneD, NekDouble> wbc(nfq,0.0);
616  bc->GetCoords(x,y,z);
617 
618  // Add edge values (trace) into the wbc
619  elmt->GetTracePhysVals(boundary,bc,w,wbc);
620 
621  //Compute womersley solution
622  for (j=0;j<nfq;j++)
623  {
624  //NOTE: only need to calculate these two once, could
625  //be stored or precomputed?
626  r = sqrt((x[j]-x0[0])*(x[j]-x0[0]) +
627  (y[j]-x0[1])*(y[j]-x0[1]) +
628  (z[j]-x0[2])*(z[j]-x0[2]))/R;
629 
630  wbc[j] = WomParam->m_wom_vel_r[0]*(1. - r*r); // Compute Poiseulle Flow
631 
632  for (k=1; k<M; k++)
633  {
634  kt = 2.0 * M_PI * k * m_time / T;
635  za = alpha * sqrt((NekDouble)k/2.0) * comp_conj;
636  zar = r * za;
637  zJ0 = Polylib::ImagBesselComp(0,za);
638  zJ0r = Polylib::ImagBesselComp(0,zar);
639  zJ0rJ0 = zJ0r / zJ0;
640  zq = std::exp(zi * kt) * std::complex<NekDouble>(
641  WomParam->m_wom_vel_r[k],
642  WomParam->m_wom_vel_i[k]);
643  zvel = zq * (z1 - zJ0rJ0);
644  wbc[j] = wbc[j] + zvel.real();
645  }
646  }
647 
648  // Multiply w by normal to get u,v,w component of velocity
649  Vmath::Smul(nfq,normals[fldid],wbc,1,wbc,1);
650 
651  Bvals = BndExp[bndid]->UpdateCoeffs()+
652  BndExp[bndid]->GetCoeff_Offset(i);
653  // Push back to Coeff space
654  bc->FwdTrans(wbc,Bvals);
655  }
656  }
657 
658 
659  void IncNavierStokes::SetUpWomersley(const int bndid, std::string womStr)
660  {
661  std::string::size_type indxBeg = womStr.find_first_of(':') + 1;
662  string filename = womStr.substr(indxBeg,string::npos);
663 
664  std::complex<NekDouble> coef;
665 
666 #if 1
667  TiXmlDocument doc(filename);
668 
669  bool loadOkay = doc.LoadFile();
670  ASSERTL0(loadOkay,(std::string("Failed to load file: ") +
671  filename).c_str());
672 
673  TiXmlHandle docHandle(&doc);
674 
675  int err; /// Error value returned by TinyXML.
676 
677  TiXmlElement *nektar = doc.FirstChildElement("NEKTAR");
678  ASSERTL0(nektar, "Unable to find NEKTAR tag in file.");
679 
680  TiXmlElement *wombc = nektar->FirstChildElement("WOMERSLEYBC");
681  ASSERTL0(wombc, "Unable to find WOMERSLEYBC tag in file.");
682 
683  // read womersley parameters
684  TiXmlElement *womparam = wombc->FirstChildElement("WOMPARAMS");
685  ASSERTL0(womparam, "Unable to find WOMPARAMS tag in file.");
686 
687  // Input coefficients
688  TiXmlElement *params = womparam->FirstChildElement("W");
689  map<std::string,std::string> Wparams;
690 
691  // read parameter list
692  while (params)
693  {
694 
695  std::string propstr;
696  propstr = params->Attribute("PROPERTY");
697 
698  ASSERTL0(!propstr.empty(),"Failed to read PROPERTY value Womersley BC Parameter");
699 
700 
701  std::string valstr;
702  valstr = params->Attribute("VALUE");
703 
704  ASSERTL0(!valstr.empty(),"Failed to read VALUE value Womersley BC Parameter");
705 
706  std::transform(propstr.begin(),propstr.end(),propstr.begin(),
707  ::toupper);
708  Wparams[propstr] = valstr;
709 
710  params = params->NextSiblingElement("W");
711  }
712 
713  // Read parameters
714 
715  ASSERTL0(Wparams.count("RADIUS") == 1,
716  "Failed to find Radius parameter in Womersley boundary conditions");
717  std::vector<NekDouble> rad;
719  Wparams["RADIUS"].c_str(),rad);
720  m_womersleyParams[bndid]->m_radius = rad[0];
721 
722  ASSERTL0(Wparams.count("PERIOD") == 1,
723  "Failed to find period parameter in Womersley boundary conditions");
724  std::vector<NekDouble> period;
726  Wparams["PERIOD"].c_str(),period);
727  m_womersleyParams[bndid]->m_period = period[0];
728 
729 
730  ASSERTL0(Wparams.count("AXISNORMAL") == 1,
731  "Failed to find axisnormal parameter in Womersley boundary conditions");
732  std::vector<NekDouble> anorm;
734  Wparams["AXISNORMAL"].c_str(),anorm);
735  m_womersleyParams[bndid]->m_axisnormal[0] = anorm[0];
736  m_womersleyParams[bndid]->m_axisnormal[1] = anorm[1];
737  m_womersleyParams[bndid]->m_axisnormal[2] = anorm[2];
738 
739 
740  ASSERTL0(Wparams.count("AXISPOINT") == 1,
741  "Failed to find axispoint parameter in Womersley boundary conditions");
742  std::vector<NekDouble> apt;
744  Wparams["AXISPOINT"].c_str(),apt);
745  m_womersleyParams[bndid]->m_axispoint[0] = apt[0];
746  m_womersleyParams[bndid]->m_axispoint[1] = apt[1];
747  m_womersleyParams[bndid]->m_axispoint[2] = apt[2];
748 
749  // Read Temporal Foruier Coefficients.
750 
751  // Find the FourierCoeff tag
752  TiXmlElement *coeff = wombc->FirstChildElement("FOURIERCOEFFS");
753 
754  // Input coefficients
755  TiXmlElement *fval = coeff->FirstChildElement("F");
756 
757  int indx;
758  int nextFourierCoeff = -1;
759 
760  while (fval)
761  {
762  nextFourierCoeff++;
763 
764  TiXmlAttribute *fvalAttr = fval->FirstAttribute();
765  std::string attrName(fvalAttr->Name());
766 
767  ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
768 
769  err = fvalAttr->QueryIntValue(&indx);
770  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
771 
772  std::string coeffStr = fval->FirstChild()->ToText()->ValueStr();
773  vector<NekDouble> coeffvals;
774  bool parseGood = ParseUtils::GenerateUnOrderedVector(coeffStr.c_str(),
775  coeffvals);
776  ASSERTL0(parseGood,(std::string("Problem reading value of fourier coefficient, ID=") + boost::lexical_cast<string>(indx)).c_str());
777  ASSERTL1(coeffvals.size() == 2,(std::string("Have not read two entries of Fourier coefficicent from ID="+ boost::lexical_cast<string>(indx)).c_str()));
778  m_womersleyParams[bndid]->m_wom_vel_r.push_back(coeffvals[0]);
779  m_womersleyParams[bndid]->m_wom_vel_i.push_back(coeffvals[1]);
780 
781  fval = fval->NextSiblingElement("F");
782  }
783 
784 #else
785  std::ifstream file(filename);
786  std::string line;
787 
788  ASSERTL1(file.is_open(),(std::string("Missing file ") + filename).c_str());
789  int count = 0;
790  while(std::getline(file,line))
791  {
792  std::stringstream stream(line);
793  while(stream>>coef)
794  {
795  m_womersleyParams[bndid]->m_wom_vel_r.push_back(coef.real());
796  m_womersleyParams[bndid]->m_wom_vel_i.push_back(coef.imag());
797  count++;
798  }
799  }
800 #endif
801  }
802 
803  /**
804  * Add an additional forcing term programmatically.
805  */
807  {
808  m_forcing.push_back(pForce);
809  }
810 
811 
812  /**
813  * Decide if at a steady state if the discrerte L2 sum of the
814  * coefficients is the same as the previous step to within the
815  * tolerance m_steadyStateTol;
816  */
818  {
819  static NekDouble previousL2 = 0.0;
820  bool returnval = false;
821 
822  NekDouble L2 = 0.0;
823 
824  // calculate L2 discrete summation
825  int ncoeffs = m_fields[0]->GetNcoeffs();
826 
827  for(int i = 0; i < m_fields.num_elements(); ++i)
828  {
829  L2 += Vmath::Dot(ncoeffs,m_fields[i]->GetCoeffs(),1,m_fields[i]->GetCoeffs(),1);
830  }
831 
832  if(fabs(L2-previousL2) < ncoeffs*m_steadyStateTol)
833  {
834  returnval = true;
835  }
836 
837  previousL2 = L2;
838 
839  return returnval;
840  }
841 
842  /**
843  *
844  */
846  {
847  int n_vel = m_velocity.num_elements();
848  int n_element = m_fields[0]->GetExpSize();
849 
850  const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
851  Array<OneD, int> ExpOrderList (n_element, ExpOrder);
852 
853  const NekDouble cLambda = 0.2; // Spencer book pag. 317
854 
855  Array<OneD, NekDouble> cfl (n_element, 0.0);
856  Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
858 
859  if(m_HomogeneousType == eHomogeneous1D) // just do check on 2D info
860  {
861  velfields = Array<OneD, Array<OneD, NekDouble> >(2);
862 
863  for(int i = 0; i < 2; ++i)
864  {
865  velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
866  }
867  }
868  else
869  {
870  velfields = Array<OneD, Array<OneD, NekDouble> >(n_vel);
871 
872  for(int i = 0; i < n_vel; ++i)
873  {
874  velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
875  }
876  }
877 
878  stdVelocity = m_extrapolation->GetMaxStdVelocity(velfields);
879 
880  for(int el = 0; el < n_element; ++el)
881  {
882  cfl[el] = m_timestep*(stdVelocity[el] * cLambda *
883  (ExpOrder[el]-1) * (ExpOrder[el]-1));
884  }
885 
886  return cfl;
887  }
888 
889  /**
890  *
891  */
893  {
894  int n_element = m_fields[0]->GetExpSize();
895 
897 
898  elmtid = Vmath::Imax(n_element,cfl,1);
899  NekDouble CFL,CFL_loc;
900 
901  CFL = CFL_loc = cfl[elmtid];
902  m_comm->AllReduce(CFL,LibUtilities::ReduceMax);
903 
904  // unshuffle elmt id if data is not stored in consecutive order.
905  elmtid = m_fields[0]->GetExp(elmtid)->GetGeom()->GetGlobalID();
906  if(CFL != CFL_loc)
907  {
908  elmtid = -1;
909  }
910 
911  m_comm->AllReduce(elmtid,LibUtilities::ReduceMax);
912 
913  // express element id with respect to plane
915  {
916  elmtid = elmtid%m_fields[0]->GetPlane(0)->GetExpSize();
917  }
918  return CFL;
919  }
920 
921 
922  /**
923  * Perform the extrapolation.
924  */
926  {
927  m_extrapolation->SubStepSaveFields(step);
928  m_extrapolation->SubStepAdvance(m_intSoln,step,m_time);
930  return false;
931  }
932 
933 
934  /**
935  * Estimate CFL and perform steady-state check
936  */
938  {
939  if(m_cflsteps && !((step+1)%m_cflsteps))
940  {
941  int elmtid;
942  NekDouble cfl = GetCFLEstimate(elmtid);
943 
944  if(m_comm->GetRank() == 0)
945  {
946  cout << "CFL (zero plane): "<< cfl << " (in elmt "
947  << elmtid << ")" << endl;
948  }
949  }
950 
951  if(m_steadyStateSteps && step && (!((step+1)%m_steadyStateSteps)))
952  {
953  if(CalcSteadyState() == true)
954  {
955  cout << "Reached Steady State to tolerance "
956  << m_steadyStateTol << endl;
957  return true;
958  }
959  }
960 
961  return false;
962  }
963 } //end of namespace
964 
EquationType m_equationType
equation type;
void SetZeroNormalVelocity()
Set Normal Velocity Component to Zero.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
Mapping from BCs to Elmt Edge IDs.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
virtual bool v_PreIntegrate(int step)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekDouble m_time
Current time of simulation.
void SetRadiationBoundaryForcing(int fieldid)
Set Radiation forcing term.
NekDouble m_kinvis
Kinematic viscosity.
bool CalcSteadyState(void)
evaluate steady state
NekDouble m_timestep
Time step size.
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp()
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
ExtrapolateSharedPtr m_extrapolation
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:442
SOLVER_UTILS_EXPORT typedef boost::shared_ptr< Forcing > ForcingSharedPtr
A shared pointer to an EquationSystem object.
Definition: Forcing.h:51
STL namespace.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:86
virtual bool v_PostIntegrate(int step)
int m_nConvectiveFields
Number of fields to be convected;.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:755
LibUtilities::CommSharedPtr m_comm
Communicator.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
RHS Factor for Radiation Condition.
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:213
void SetUpWomersley(const int bndid, std::string womstr)
Set Up Womersley details.
int m_cflsteps
dump cfl estimate
std::map< int, WomersleyParamsSharedPtr > m_womersleyParams
Womersley parameters if required.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Base class for unsteady solvers.
boost::shared_ptr< WomersleyParams > WomersleyParamsSharedPtr
NekDouble Evaluate() const
Definition: Equation.h:102
int m_steadyStateSteps
Check for steady state at step interval.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
int m_spacedim
Spatial dimension (>= expansion dim).
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
virtual void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:914
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
static NekDouble rad(NekDouble x, NekDouble y)
SOLVER_UTILS_EXPORT int GetNpoints()
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
Mapping from BCs to Elmt IDs.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
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:2999
const std::string kEquationTypeStr[]
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_steadyStateTol
Tolerance to which steady state should be evaluated at.
void SetWomersleyBoundary(const int fldid, const int bndid)
Set Womersley Profile if specified.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
Definition: ParseUtils.hpp:128
A base class for PDEs which include an advection component.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
virtual void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
virtual void v_InitObject()
Init object for UnsteadySystem class.
virtual int v_GetForceDimension()=0
enum HomogeneousType m_HomogeneousType
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:183
NekDouble GetCFLEstimate(int &elmtid)
Array< OneD, NekDouble > GetElmtCFLVals(void)