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
58namespace Nektar
59{
60
61std::string IncNavierStokes::eqTypeLookupIds[6] = {
67 "EqType", "SteadyNavierStokes", eSteadyNavierStokes),
69 "EqType", "SteadyLinearisedNS", eSteadyLinearisedNS),
71 "EqType", "UnsteadyNavierStokes", eUnsteadyNavierStokes),
72 LibUtilities::SessionReader::RegisterEnumValue("EqType", "UnsteadyStokes",
74
75/**
76 * Constructor. Creates ...
77 *
78 * \param
79 * \param
80 */
84 : UnsteadySystem(pSession, pGraph), AdvectionSystem(pSession, pGraph),
85 m_SmoothAdvection(false)
86{
87}
88
89void IncNavierStokes::v_InitObject(bool DeclareField)
90{
92
93 int i, j;
94 int numfields = m_fields.size();
95 std::string velids[] = {"u", "v", "w"};
96
97 // Set up Velocity field to point to the first m_expdim of m_fields;
99
100 for (i = 0; i < m_spacedim; ++i)
101 {
102 for (j = 0; j < numfields; ++j)
103 {
104 std::string var = m_boundaryConditions->GetVariable(j);
105 if (boost::iequals(velids[i], var))
106 {
107 m_velocity[i] = j;
108 break;
109 }
110
111 ASSERTL0(j != numfields, "Failed to find field: " + var);
112 }
113 }
114
115 // Set up equation type enum using kEquationTypeStr
116 for (i = 0; i < (int)eEquationTypeSize; ++i)
117 {
118 bool match;
119 m_session->MatchSolverInfo("EQTYPE", kEquationTypeStr[i], match, false);
120 if (match)
121 {
123 break;
124 }
125 }
126 ASSERTL0(i != eEquationTypeSize, "EQTYPE not found in SOLVERINFO section");
127
128 m_session->LoadParameter("Kinvis", m_kinvis);
129
130 // Default advection type per solver
131 std::string vConvectiveType;
132 switch (m_equationType)
133 {
134 case eUnsteadyStokes:
136 vConvectiveType = "NoAdvection";
137 break;
140 vConvectiveType = "Convective";
141 break;
143 vConvectiveType = "Linearised";
144 break;
145 default:
146 break;
147 }
148
149 // Check if advection type overridden
150 if (m_session->DefinesTag("AdvectiveType") &&
153 {
154 vConvectiveType = m_session->GetTag("AdvectiveType");
155 }
156
157 // Initialise advection
159 vConvectiveType, vConvectiveType);
160 m_advObject->InitObject(m_session, m_fields);
161
162 // Set up arrays for moving reference frame
163 // Note: this must be done before the forcing
164 if (DefinedForcing("MovingReferenceFrame"))
165 {
166 // 0-5(inertial disp), 6-11(inertial vel), 12-17(inertial acce) current
167 // 18-21(body pivot)
168 // 21-26(inertial disp), 27-32(body vel), 33-38(body acce) next step
169 // 39-41(body pivot)
171 "X", "Y", "Z", "Theta_x", "Theta_y", "Theta_z",
172 "U", "V", "W", "Omega_x", "Omega_y", "Omega_z",
173 "A_x", "A_y", "A_z", "DOmega_x", "DOmega_y", "DOmega_z",
174 "X0", "Y0", "Z0"};
177 }
178
179 // Forcing terms
180 m_forcing = SolverUtils::Forcing::Load(m_session, shared_from_this(),
182
183 // check to see if any Robin boundary conditions and if so set
184 // up m_field to boundary condition maps;
188
189 for (size_t i = 0; i < m_fields.size(); ++i)
190 {
191 bool Set = false;
192
195 int radpts = 0;
196
197 BndConds = m_fields[i]->GetBndConditions();
198 BndExp = m_fields[i]->GetBndCondExpansions();
199 for (size_t n = 0; n < BndConds.size(); ++n)
200 {
201 if (boost::iequals(BndConds[n]->GetUserDefined(), "Radiation"))
202 {
203 ASSERTL0(
204 BndConds[n]->GetBoundaryConditionType() ==
206 "Radiation boundary condition must be of type Robin <R>");
207
208 if (Set == false)
209 {
210 m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],
212 Set = true;
213 }
214 radpts += BndExp[n]->GetTotPoints();
215 }
216 if (boost::iequals(BndConds[n]->GetUserDefined(),
217 "ZeroNormalComponent"))
218 {
219 ASSERTL0(BndConds[n]->GetBoundaryConditionType() ==
221 "Zero Normal Component boundary condition option must "
222 "be of type Dirichlet <D>");
223
224 if (Set == false)
225 {
226 m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],
228 Set = true;
229 }
230 }
231 }
232
234
235 radpts = 0; // reset to use as a counter
236
237 for (size_t n = 0; n < BndConds.size(); ++n)
238 {
239 if (boost::iequals(BndConds[n]->GetUserDefined(), "Radiation"))
240 {
241
242 int npoints = BndExp[n]->GetNpoints();
243 Array<OneD, NekDouble> x0(npoints, 0.0);
244 Array<OneD, NekDouble> x1(npoints, 0.0);
245 Array<OneD, NekDouble> x2(npoints, 0.0);
246 Array<OneD, NekDouble> tmpArray;
247
248 BndExp[n]->GetCoords(x0, x1, x2);
249
251 std::static_pointer_cast<
253 ->m_robinPrimitiveCoeff;
254
255 coeff.Evaluate(x0, x1, x2, m_time,
256 tmpArray = m_fieldsRadiationFactor[i] + radpts);
257 // Vmath::Neg(npoints,tmpArray = m_fieldsRadiationFactor[i]+
258 // radpts,1);
259 radpts += npoints;
260 }
261 }
262 }
263
264 // Set up maping for womersley BC - and load variables
265 for (size_t i = 0; i < m_fields.size(); ++i)
266 {
267 for (size_t n = 0; n < m_fields[i]->GetBndConditions().size(); ++n)
268 {
269 if (boost::istarts_with(
270 m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
271 "Womersley"))
272 {
273 // assumes that boundary condition is applied in normal
274 // direction and is decomposed for each direction. There could
275 // be a unique file for each direction
276 m_womersleyParams[i][n] =
278 m_spacedim);
279 // Read in fourier coeffs and precompute coefficients
281 i, n, m_fields[i]->GetBndConditions()[n]->GetUserDefined());
282
283 m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],
285 }
286 }
287 }
288
289 // Set up Field Meta Data for output files
290 m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
291 m_fieldMetaDataMap["TimeStep"] =
292 boost::lexical_cast<std::string>(m_timestep);
293}
294
295/**
296 * Evaluation -N(V) for all fields except pressure using m_velocity
297 */
299 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
300 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
301{
302 size_t VelDim = m_velocity.size();
303 Array<OneD, Array<OneD, NekDouble>> velocity(VelDim);
304
305 size_t npoints = m_fields[0]->GetNpoints();
306 for (size_t i = 0; i < VelDim; ++i)
307 {
308 velocity[i] = Array<OneD, NekDouble>(npoints);
309 Vmath::Vcopy(npoints, inarray[m_velocity[i]], 1, velocity[i], 1);
310 }
311 for (auto &x : m_forcing)
312 {
313 x->PreApply(m_fields, velocity, velocity, time);
314 }
315
316 m_advObject->Advect(m_nConvectiveFields, m_fields, velocity, inarray,
317 outarray, time);
318}
319
320/**
321 * Time dependent boundary conditions updating
322 */
324{
325 size_t i, n;
326 std::string varName;
327 size_t nvariables = m_fields.size();
328
329 for (i = 0; i < nvariables; ++i)
330 {
331 for (n = 0; n < m_fields[i]->GetBndConditions().size(); ++n)
332 {
333 if (m_fields[i]->GetBndConditions()[n]->IsTimeDependent())
334 {
335 varName = m_session->GetVariable(i);
336 m_fields[i]->EvaluateBoundaryConditions(time, varName);
337 }
338 else if (boost::istarts_with(
339 m_fields[i]->GetBndConditions()[n]->GetUserDefined(),
340 "Womersley"))
341 {
343 }
344 }
345
346 // Set Radiation conditions if required
348 }
349
350 // Enforcing the boundary conditions (Inlet and wall) for the
351 // Moving reference frame
353}
354
355/**
356 * Probably should be pushed back into ContField?
357 */
359{
360 size_t i, n;
361
364
365 BndConds = m_fields[fieldid]->GetBndConditions();
366 BndExp = m_fields[fieldid]->GetBndCondExpansions();
367
370
371 size_t cnt;
372 size_t elmtid, nq, offset, boundary;
373 Array<OneD, NekDouble> Bvals, U;
374 size_t cnt1 = 0;
375
376 for (cnt = n = 0; n < BndConds.size(); ++n)
377 {
378 std::string type = BndConds[n]->GetUserDefined();
379
380 if ((BndConds[n]->GetBoundaryConditionType() ==
382 (boost::iequals(type, "Radiation")))
383 {
384 size_t nExp = BndExp[n]->GetExpSize();
385 for (i = 0; i < nExp; ++i, cnt++)
386 {
387 elmtid = m_fieldsBCToElmtID[m_velocity[fieldid]][cnt];
388 elmt = m_fields[fieldid]->GetExp(elmtid);
389 offset = m_fields[fieldid]->GetPhys_Offset(elmtid);
390
391 U = m_fields[fieldid]->UpdatePhys() + offset;
392 Bc = BndExp[n]->GetExp(i);
393
394 boundary = m_fieldsBCToTraceID[fieldid][cnt];
395
396 // Get edge values and put into ubc
397 nq = Bc->GetTotPoints();
399 elmt->GetTracePhysVals(boundary, Bc, U, ubc);
400
401 Vmath::Vmul(nq,
403 [fieldid][cnt1 + BndExp[n]->GetPhys_Offset(i)],
404 1, &ubc[0], 1, &ubc[0], 1);
405
406 Bvals =
407 BndExp[n]->UpdateCoeffs() + BndExp[n]->GetCoeff_Offset(i);
408
409 Bc->IProductWRTBase(ubc, Bvals);
410 }
411 cnt1 += BndExp[n]->GetTotPoints();
412 }
413 else
414 {
415 cnt += BndExp[n]->GetExpSize();
416 }
417 }
418}
419
421{
422 // use static trip since cannot use UserDefinedTag for zero
423 // velocity and have time dependent conditions
424 static bool Setup = false;
425
426 if (Setup == true)
427 {
428 return;
429 }
430 Setup = true;
431
432 size_t i, n;
433
435 BndConds(m_spacedim);
437
438 for (int i = 0; i < m_spacedim; ++i)
439 {
440 BndConds[i] = m_fields[m_velocity[i]]->GetBndConditions();
441 BndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
442 }
443
445
446 size_t cnt;
447 size_t elmtid, nq, boundary;
448
450 Array<OneD, NekDouble> Bphys, Bcoeffs;
451
452 size_t fldid = m_velocity[0];
453
454 for (cnt = n = 0; n < BndConds[0].size(); ++n)
455 {
456 if ((BndConds[0][n]->GetBoundaryConditionType() ==
458 (boost::iequals(BndConds[0][n]->GetUserDefined(),
459 "ZeroNormalComponent")))
460 {
461 size_t nExp = BndExp[0][n]->GetExpSize();
462 for (i = 0; i < nExp; ++i, cnt++)
463 {
464 elmtid = m_fieldsBCToElmtID[fldid][cnt];
465 elmt = m_fields[0]->GetExp(elmtid);
466 boundary = m_fieldsBCToTraceID[fldid][cnt];
467
468 normals = elmt->GetTraceNormal(boundary);
469
470 nq = BndExp[0][n]->GetExp(i)->GetTotPoints();
471 Array<OneD, NekDouble> normvel(nq, 0.0);
472
473 for (int k = 0; k < m_spacedim; ++k)
474 {
475 Bphys = BndExp[k][n]->UpdatePhys() +
476 BndExp[k][n]->GetPhys_Offset(i);
477 Bc = BndExp[k][n]->GetExp(i);
478 Vmath::Vvtvp(nq, normals[k], 1, Bphys, 1, normvel, 1,
479 normvel, 1);
480 }
481
482 // negate normvel for next step
483 Vmath::Neg(nq, normvel, 1);
484
485 for (int k = 0; k < m_spacedim; ++k)
486 {
487 Bphys = BndExp[k][n]->UpdatePhys() +
488 BndExp[k][n]->GetPhys_Offset(i);
489 Bcoeffs = BndExp[k][n]->UpdateCoeffs() +
490 BndExp[k][n]->GetCoeff_Offset(i);
491 Bc = BndExp[k][n]->GetExp(i);
492 Vmath::Vvtvp(nq, normvel, 1, normals[k], 1, Bphys, 1, Bphys,
493 1);
494 Bc->FwdTransBndConstrained(Bphys, Bcoeffs);
495 }
496 }
497 }
498 else
499 {
500 cnt += BndExp[0][n]->GetExpSize();
501 }
502 }
503}
504
505/**
506 * Womersley boundary condition defintion
507 */
508void IncNavierStokes::SetWomersleyBoundary(const int fldid, const int bndid)
509{
510 ASSERTL1(m_womersleyParams.count(bndid) == 1,
511 "Womersley parameters for this boundary have not been set up");
512
513 WomersleyParamsSharedPtr WomParam = m_womersleyParams[fldid][bndid];
514 NekComplexDouble zvel;
515 size_t i, j, k;
516
517 size_t M_coeffs = WomParam->m_wom_vel.size();
518
519 NekDouble T = WomParam->m_period;
520 NekDouble axis_normal = WomParam->m_axisnormal[fldid];
521
522 // Womersley Number
523 NekComplexDouble omega_c(2.0 * M_PI / T, 0.0);
524 NekComplexDouble k_c(0.0, 0.0);
525 NekComplexDouble m_time_c(m_time, 0.0);
526 NekComplexDouble zi(0.0, 1.0);
527 NekComplexDouble i_pow_3q2(-1.0 / sqrt(2.0), 1.0 / sqrt(2.0));
528
530 BndCondExp = m_fields[fldid]->GetBndCondExpansions()[bndid];
531
533 size_t cnt = 0;
534 size_t nfq;
536 size_t exp_npts = BndCondExp->GetExpSize();
537 Array<OneD, NekDouble> wbc(exp_npts, 0.0);
538
540
541 // preallocate the exponent
542 for (k = 1; k < M_coeffs; k++)
543 {
544 k_c = NekComplexDouble((NekDouble)k, 0.0);
545 zt[k] = std::exp(zi * omega_c * k_c * m_time_c);
546 }
547
548 // Loop over each element in an expansion
549 for (i = 0; i < exp_npts; ++i, cnt++)
550 {
551 // Get Boundary and trace expansion
552 bc = BndCondExp->GetExp(i);
553 nfq = bc->GetTotPoints();
554 Array<OneD, NekDouble> wbc(nfq, 0.0);
555
556 // Compute womersley solution
557 for (j = 0; j < nfq; j++)
558 {
559 wbc[j] = WomParam->m_poiseuille[i][j];
560 for (k = 1; k < M_coeffs; k++)
561 {
562 zvel = WomParam->m_zvel[i][j][k] * zt[k];
563 wbc[j] = wbc[j] + zvel.real();
564 }
565 }
566
567 // Multiply w by normal to get u,v,w component of velocity
568 Vmath::Smul(nfq, axis_normal, wbc, 1, wbc, 1);
569 // get the offset
570 Bvals = BndCondExp->UpdateCoeffs() + BndCondExp->GetCoeff_Offset(i);
571
572 // Push back to Coeff space
573 bc->FwdTrans(wbc, Bvals);
574 }
575}
576
577void IncNavierStokes::SetUpWomersley(const int fldid, const int bndid,
578 std::string womStr)
579{
580 std::string::size_type indxBeg = womStr.find_first_of(':') + 1;
581 std::string filename = womStr.substr(indxBeg, std::string::npos);
582
583 TiXmlDocument doc(filename);
584
585 bool loadOkay = doc.LoadFile();
586 ASSERTL0(loadOkay,
587 (std::string("Failed to load file: ") + filename).c_str());
588
589 TiXmlHandle docHandle(&doc);
590
591 int err; /// Error value returned by TinyXML.
592
593 TiXmlElement *nektar = doc.FirstChildElement("NEKTAR");
594 ASSERTL0(nektar, "Unable to find NEKTAR tag in file.");
595
596 TiXmlElement *wombc = nektar->FirstChildElement("WOMERSLEYBC");
597 ASSERTL0(wombc, "Unable to find WOMERSLEYBC tag in file.");
598
599 // read womersley parameters
600 TiXmlElement *womparam = wombc->FirstChildElement("WOMPARAMS");
601 ASSERTL0(womparam, "Unable to find WOMPARAMS tag in file.");
602
603 // Input coefficients
604 TiXmlElement *params = womparam->FirstChildElement("W");
605 std::map<std::string, std::string> Wparams;
606
607 // read parameter list
608 while (params)
609 {
610
611 std::string propstr;
612 propstr = params->Attribute("PROPERTY");
613
614 ASSERTL0(!propstr.empty(),
615 "Failed to read PROPERTY value Womersley BC Parameter");
616
617 std::string valstr;
618 valstr = params->Attribute("VALUE");
619
620 ASSERTL0(!valstr.empty(),
621 "Failed to read VALUE value Womersley BC Parameter");
622
623 std::transform(propstr.begin(), propstr.end(), propstr.begin(),
624 ::toupper);
625 Wparams[propstr] = valstr;
626
627 params = params->NextSiblingElement("W");
628 }
629 bool parseGood;
630
631 // Read parameters
632
633 ASSERTL0(
634 Wparams.count("RADIUS") == 1,
635 "Failed to find Radius parameter in Womersley boundary conditions");
636 std::vector<NekDouble> rad;
637 ParseUtils::GenerateVector(Wparams["RADIUS"], rad);
638 m_womersleyParams[fldid][bndid]->m_radius = rad[0];
639
640 ASSERTL0(
641 Wparams.count("PERIOD") == 1,
642 "Failed to find period parameter in Womersley boundary conditions");
643 std::vector<NekDouble> period;
644 parseGood = ParseUtils::GenerateVector(Wparams["PERIOD"], period);
645 m_womersleyParams[fldid][bndid]->m_period = period[0];
646
647 ASSERTL0(
648 Wparams.count("AXISNORMAL") == 1,
649 "Failed to find axisnormal parameter in Womersley boundary conditions");
650 std::vector<NekDouble> anorm;
651 parseGood = ParseUtils::GenerateVector(Wparams["AXISNORMAL"], anorm);
652 m_womersleyParams[fldid][bndid]->m_axisnormal[0] = anorm[0];
653 m_womersleyParams[fldid][bndid]->m_axisnormal[1] = anorm[1];
654 m_womersleyParams[fldid][bndid]->m_axisnormal[2] = anorm[2];
655
656 ASSERTL0(
657 Wparams.count("AXISPOINT") == 1,
658 "Failed to find axispoint parameter in Womersley boundary conditions");
659 std::vector<NekDouble> apt;
660 parseGood = ParseUtils::GenerateVector(Wparams["AXISPOINT"], apt);
661 m_womersleyParams[fldid][bndid]->m_axispoint[0] = apt[0];
662 m_womersleyParams[fldid][bndid]->m_axispoint[1] = apt[1];
663 m_womersleyParams[fldid][bndid]->m_axispoint[2] = apt[2];
664
665 // Read Temporal Fourier Coefficients.
666
667 // Find the FourierCoeff tag
668 TiXmlElement *coeff = wombc->FirstChildElement("FOURIERCOEFFS");
669
670 // Input coefficients
671 TiXmlElement *fval = coeff->FirstChildElement("F");
672
673 int indx;
674
675 while (fval)
676 {
677 TiXmlAttribute *fvalAttr = fval->FirstAttribute();
678 std::string attrName(fvalAttr->Name());
679
680 ASSERTL0(attrName == "ID",
681 (std::string("Unknown attribute name: ") + attrName).c_str());
682
683 err = fvalAttr->QueryIntValue(&indx);
684 ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
685
686 std::string coeffStr = fval->FirstChild()->ToText()->ValueStr();
687 std::vector<NekDouble> coeffvals;
688
689 parseGood = ParseUtils::GenerateVector(coeffStr, coeffvals);
690 ASSERTL0(
691 parseGood,
692 (std::string("Problem reading value of fourier coefficient, ID=") +
693 boost::lexical_cast<std::string>(indx))
694 .c_str());
695 ASSERTL1(
696 coeffvals.size() == 2,
697 (std::string(
698 "Have not read two entries of Fourier coefficicent from ID=" +
699 boost::lexical_cast<std::string>(indx))
700 .c_str()));
701
702 m_womersleyParams[fldid][bndid]->m_wom_vel.push_back(
703 NekComplexDouble(coeffvals[0], coeffvals[1]));
704
705 fval = fval->NextSiblingElement("F");
706 }
707
708 // starting point of precalculation
709 size_t i, j, k;
710 // M fourier coefficients
711 size_t M_coeffs = m_womersleyParams[fldid][bndid]->m_wom_vel.size();
712 NekDouble R = m_womersleyParams[fldid][bndid]->m_radius;
713 NekDouble T = m_womersleyParams[fldid][bndid]->m_period;
714 Array<OneD, NekDouble> x0 = m_womersleyParams[fldid][bndid]->m_axispoint;
715
717 // Womersley Number
718 NekComplexDouble omega_c(2.0 * M_PI / T, 0.0);
719 NekComplexDouble alpha_c(R * sqrt(omega_c.real() / m_kinvis), 0.0);
720 NekComplexDouble z1(1.0, 0.0);
721 NekComplexDouble i_pow_3q2(-1.0 / sqrt(2.0), 1.0 / sqrt(2.0));
722
724 BndCondExp = m_fields[fldid]->GetBndCondExpansions()[bndid];
725
727 size_t cnt = 0;
728 size_t nfq;
730
731 size_t exp_npts = BndCondExp->GetExpSize();
732 Array<OneD, NekDouble> wbc(exp_npts, 0.0);
733
734 // allocate time indepedent variables
735 m_womersleyParams[fldid][bndid]->m_poiseuille =
737 m_womersleyParams[fldid][bndid]->m_zvel =
739 // could use M_coeffs - 1 but need to avoid complicating things
740 Array<OneD, NekComplexDouble> zJ0(M_coeffs);
741 Array<OneD, NekComplexDouble> lamda_n(M_coeffs);
742 Array<OneD, NekComplexDouble> k_c(M_coeffs);
743 NekComplexDouble zJ0r;
744
745 for (k = 1; k < M_coeffs; k++)
746 {
747 k_c[k] = NekComplexDouble((NekDouble)k, 0.0);
748 lamda_n[k] = i_pow_3q2 * alpha_c * sqrt(k_c[k]);
749 zJ0[k] = Polylib::ImagBesselComp(0, lamda_n[k]);
750 }
751
752 // Loop over each element in an expansion
753 for (i = 0; i < exp_npts; ++i, cnt++)
754 {
755 // Get Boundary and trace expansion
756 bc = BndCondExp->GetExp(i);
757 nfq = bc->GetTotPoints();
758
759 Array<OneD, NekDouble> x(nfq, 0.0);
760 Array<OneD, NekDouble> y(nfq, 0.0);
761 Array<OneD, NekDouble> z(nfq, 0.0);
762 bc->GetCoords(x, y, z);
763
764 m_womersleyParams[fldid][bndid]->m_poiseuille[i] =
766 m_womersleyParams[fldid][bndid]->m_zvel[i] =
768
769 // Compute coefficients
770 for (j = 0; j < nfq; j++)
771 {
772 rqR = NekComplexDouble(sqrt((x[j] - x0[0]) * (x[j] - x0[0]) +
773 (y[j] - x0[1]) * (y[j] - x0[1]) +
774 (z[j] - x0[2]) * (z[j] - x0[2])) /
775 R,
776 0.0);
777
778 // Compute Poiseulle Flow
779 m_womersleyParams[fldid][bndid]->m_poiseuille[i][j] =
780 m_womersleyParams[fldid][bndid]->m_wom_vel[0].real() *
781 (1. - rqR.real() * rqR.real());
782
783 m_womersleyParams[fldid][bndid]->m_zvel[i][j] =
785
786 // compute the velocity information
787 for (k = 1; k < M_coeffs; k++)
788 {
789 zJ0r = Polylib::ImagBesselComp(0, rqR * lamda_n[k]);
790 m_womersleyParams[fldid][bndid]->m_zvel[i][j][k] =
791 m_womersleyParams[fldid][bndid]->m_wom_vel[k] *
792 (z1 - (zJ0r / zJ0[k]));
793 }
794 }
795 }
796}
797
798/**
799 * Add an additional forcing term programmatically.
800 */
802{
803 m_forcing.push_back(pForce);
804}
805
806/**
807 *
808 */
810 [[maybe_unused]] const NekDouble SpeedSoundFactor)
811{
812 size_t nvel = m_velocity.size();
813 size_t nelmt = m_fields[0]->GetExpSize();
814
815 Array<OneD, NekDouble> stdVelocity(nelmt, 0.0);
817
818 if (m_HomogeneousType == eHomogeneous1D) // just do check on 2D info
819 {
821
822 for (size_t i = 0; i < 2; ++i)
823 {
824 velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
825 }
826 }
827 else
828 {
829 velfields = Array<OneD, Array<OneD, NekDouble>>(nvel);
830
831 for (size_t i = 0; i < nvel; ++i)
832 {
833 velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
834 }
835 }
836
837 stdVelocity = m_extrapolation->GetMaxStdVelocity(velfields);
838
839 return stdVelocity;
840}
841
842/**
843 *
844 */
846 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
848{
849 if (physfield.size())
850 {
851 pressure = physfield[physfield.size() - 1];
852 }
853}
854
855/**
856 *
857 */
859 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
860 Array<OneD, NekDouble> &density)
861{
862 int nPts = physfield[0].size();
863 Vmath::Fill(nPts, 1.0, density, 1);
864}
865
866/**
867 *
868 */
870 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
872{
873 for (int i = 0; i < m_spacedim; ++i)
874 {
875 velocity[i] = physfield[i];
876 }
877}
878
879/**
880 * Function to set the moving frame velocities calucated in the forcing
881 * this gives access to the moving reference forcing to set the velocities
882 * to be later used in enforcing the boundary condition in IncNavierStokes
883 * class
884 */
886 const Array<OneD, NekDouble> &vFrameVels, const int step)
887{
888 if (m_movingFrameData.size())
889 {
890 ASSERTL0(vFrameVels.size() == 12,
891 "Arrays have different dimensions, cannot set moving frame "
892 "velocities");
893 Array<OneD, NekDouble> temp = m_movingFrameData + 6 + 21 * step;
894 Vmath::Vcopy(vFrameVels.size(), vFrameVels, 1, temp, 1);
895 }
896}
897
899 Array<OneD, NekDouble> &vFrameVels, const int step)
900{
901 if (m_movingFrameData.size())
902 {
903 ASSERTL0(vFrameVels.size() == 12,
904 "Arrays have different dimensions, cannot get moving frame "
905 "velocities");
906 Vmath::Vcopy(vFrameVels.size(), m_movingFrameData + 6 + 21 * step, 1,
907 vFrameVels, 1);
908 return true;
909 }
910 else
911 {
912 return false;
913 }
914}
915
916/**
917 * Function to set the angles between the moving frame of reference and
918 * stationary inertial reference frame
919 **/
921 const Array<OneD, NekDouble> &vFrameDisp, const int step)
922{
923 if (m_movingFrameData.size())
924 {
925 ASSERTL0(
926 vFrameDisp.size() == 6,
927 "Arrays have different size, cannot set moving frame displacement");
928 Array<OneD, NekDouble> temp = m_movingFrameData + 21 * step;
929 Vmath::Vcopy(vFrameDisp.size(), vFrameDisp, 1, temp, 1);
930 }
931}
932
933/**
934 * Function to get the angles between the moving frame of reference and
935 * stationary inertial reference frame
936 **/
938 const int step)
939{
940 if (m_movingFrameData.size())
941 {
942 ASSERTL0(
943 vFrameDisp.size() == 6,
944 "Arrays have different size, cannot get moving frame displacement");
945 Vmath::Vcopy(vFrameDisp.size(), m_movingFrameData + 21 * step, 1,
946 vFrameDisp, 1);
947 return true;
948 }
949 else
950 {
951 return false;
952 }
953}
954
956 const Array<OneD, NekDouble> &vFramePivot)
957{
958 ASSERTL0(vFramePivot.size() == 3,
959 "Arrays have different size, cannot set moving frame pivot");
961 Vmath::Vcopy(vFramePivot.size(), vFramePivot, 1, temp, 1);
962 temp = m_movingFrameData + 39;
963 Vmath::Vcopy(vFramePivot.size(), vFramePivot, 1, temp, 1);
964}
965
967{
968 if (m_aeroForces.size() >= 6)
969 {
970 Vmath::Vcopy(6, forces, 1, m_aeroForces, 1);
971 }
972}
973
975{
976 if (m_aeroForces.size() >= 6)
977 {
978 Vmath::Vcopy(6, m_aeroForces, 1, forces, 1);
979 }
980}
981
982/**
983 * Function to check the type of forcing
984 **/
985bool IncNavierStokes::DefinedForcing(const std::string &sForce)
986{
987 std::vector<std::string> vForceList;
988 bool hasForce{false};
989
990 if (!m_session->DefinesElement("Nektar/Forcing"))
991 {
992 return hasForce;
993 }
994
995 TiXmlElement *vForcing = m_session->GetElement("Nektar/Forcing");
996 if (vForcing)
997 {
998 TiXmlElement *vForce = vForcing->FirstChildElement("FORCE");
999 while (vForce)
1000 {
1001 std::string vType = vForce->Attribute("TYPE");
1002
1003 vForceList.push_back(vType);
1004 vForce = vForce->NextSiblingElement("FORCE");
1005 }
1006 }
1007
1008 for (auto &f : vForceList)
1009 {
1010 if (boost::iequals(f, sForce))
1011 {
1012 hasForce = true;
1013 }
1014 }
1015
1016 return hasForce;
1017}
1018
1019/**
1020 * Perform the extrapolation.
1021 */
1023{
1024 m_extrapolation->SubStepSaveFields(step);
1025 m_extrapolation->SubStepAdvance(step, m_time);
1027 return false;
1028}
1029
1030} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
Mapping from BCs to Elmt Edge IDs.
Array< OneD, NekDouble > v_GetMaxStdVelocity(const NekDouble SpeedSoundFactor) override
std::map< int, std::map< int, WomersleyParamsSharedPtr > > m_womersleyParams
Womersley parameters if required.
void v_SetAeroForce(Array< OneD, NekDouble > forces) override
bool v_GetMovingFrameVelocities(Array< OneD, NekDouble > &vFrameVels, const int step) override
virtual int v_GetForceDimension()=0
void SetWomersleyBoundary(const int fldid, const int bndid)
Set Womersley Profile if specified.
void v_GetAeroForce(Array< OneD, NekDouble > forces) override
void SetZeroNormalVelocity()
Set Normal Velocity Component to Zero.
void v_SetMovingFrameDisp(const Array< OneD, NekDouble > &vFrameDisp, const int step) override
MultiRegions::ExpListSharedPtr v_GetPressure() override
Array< OneD, NekDouble > m_aeroForces
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
bool v_GetMovingFrameDisp(Array< OneD, NekDouble > &vFrameDisp, const int step) override
NekDouble m_kinvis
Kinematic viscosity.
void v_GetDensity(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density) override
void SetRadiationBoundaryForcing(int fieldid)
Set Radiation forcing term.
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
RHS Factor for Radiation Condition.
bool DefinedForcing(const std::string &sForce)
bool v_PreIntegrate(int step) override
void v_SetMovingFrameVelocities(const Array< OneD, NekDouble > &vFrameVels, const int step) override
void SetUpWomersley(const int fldid, const int bndid, std::string womstr)
Set Up Womersley details.
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
static std::string eqTypeLookupIds[]
void v_SetMovingFramePivot(const Array< OneD, NekDouble > &vFramePivot) override
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
Mapping from BCs to Elmt IDs.
void v_GetVelocity(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity) override
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.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
IncNavierStokes(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Constructor.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
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:130
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
int m_spacedim
Spatial dimension (>= expansion dim).
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
std::vector< std::string > m_strFrameData
variable name in m_movingFrameData
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
enum HomogeneousType m_HomogeneousType
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Array< OneD, NekDouble > m_movingFrameData
Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z,...
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:76
Base class for unsteady solvers.
static NekDouble rad(NekDouble x, NekDouble y)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
SOLVER_UTILS_EXPORT typedef std::shared_ptr< Forcing > ForcingSharedPtr
A shared pointer to an EquationSystem object.
Definition: Forcing.h:53
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::vector< double > z(NPUPPER)
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:1559
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.hpp:72
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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.hpp:366
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.hpp:100
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285