Nektar++
Extrapolate.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Extrapolate.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// 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: Abstract base class for Extrapolate.
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38namespace Nektar
39{
40
42 {1.0, 0.0, 0.0}, {2.0, -1.0, 0.0}, {3.0, -3.0, 1.0}};
44 {1.0, 0.0, 0.0}, {2.0, -0.5, 0.0}, {3.0, -1.5, 1.0 / 3.0}};
45NekDouble Extrapolate::StifflyStable_Gamma0_Coeffs[3] = {1.0, 1.5, 11.0 / 6.0};
46
48{
49 static ExtrapolateFactory instance;
50 return instance;
51}
52
56 const Array<OneD, int> pVel,
57 const SolverUtils::AdvectionSharedPtr advObject)
58 : m_session(pSession), m_fields(pFields), m_pressure(pPressure),
59 m_velocity(pVel), m_advObject(advObject)
60{
61 m_session->LoadParameter("TimeStep", m_timestep, 0.01);
62 m_comm = m_session->GetComm();
63}
64
65std::string Extrapolate::def =
67 "StandardExtrapolate", "StandardExtrapolate");
68
69/**
70 *
71 */
73{
74 if (m_numHBCDof)
75 {
76 // Update velocity BF at n+1 (actually only needs doing if
77 // velocity is time dependent on HBCs)
79
80 // Calculate acceleration term at level n based on previous steps
82
83 // Subtract acceleration term off m_pressureHBCs[nlevels-1]
85 1, m_pressureHBCs[m_intSteps - 1], 1,
87 }
88}
89
90/**
91 *
92 */
94{
95 if (m_numHBCDof)
96 {
97 int order = std::min(m_pressureCalls, m_intSteps);
98
99 // Update velocity BF at n+1 (actually only needs doing if
100 // velocity is time dependent on HBCs)
102
103 // Subtract acceleration term off m_pressureHBCs[nlevels-1]
105 -1.0 * StifflyStable_Gamma0_Coeffs[order - 1] / m_timestep,
107 m_pressureHBCs[m_intSteps - 1], 1);
108 }
109}
110
111/**
112 * Unified routine for calculation high-oder terms
113 */
115 const Array<OneD, const Array<OneD, NekDouble>> &fields,
116 const Array<OneD, const Array<OneD, NekDouble>> &N, NekDouble kinvis)
117{
118 size_t n, cnt;
119
121
124
127
128 // Loop all boundary conditions
129 for (n = cnt = 0; n < m_PBndConds.size(); ++n)
130 {
131 // Detect higher order boundary conditions
132 if ((m_hbcType[n] == eHBCNeumann) || (m_hbcType[n] == eConvectiveOBC))
133 {
134 m_bndElmtExps[n]->SetWaveSpace(m_fields[0]->GetWaveSpace());
135 int nqb = m_PBndExp[n]->GetTotPoints();
136 int nq = m_bndElmtExps[n]->GetTotPoints();
137 int ncoeffs = m_PBndExp[n]->GetNcoeffs();
138
139 for (int i = 0; i < m_bnd_dim; i++)
140 {
141 BndValues[i] = Array<OneD, NekDouble>(nqb, 0.0);
142 }
143
144 for (int i = 0; i < m_curl_dim; i++)
145 {
146 Q[i] = Array<OneD, NekDouble>(nq, 0.0);
147 }
148
149 // Obtaining fields on BndElmtExp
150 for (int i = 0; i < m_curl_dim; i++)
151 {
152 m_fields[0]->ExtractPhysToBndElmt(n, fields[i], Velocity[i]);
153 }
154
155 if (N.size()) // not required for some extrapolation
156 {
157 for (int i = 0; i < m_bnd_dim; i++)
158 {
159 m_fields[0]->ExtractPhysToBndElmt(n, N[i], Advection[i]);
160 }
161 }
162
163 // CurlCurl
164 m_bndElmtExps[n]->CurlCurl(Velocity, Q);
165
166 // Mounting advection component into the high-order condition
167 for (int i = 0; i < m_bnd_dim; i++)
168 {
169 MountHOPBCs(nq, kinvis, Q[i], Advection[i]);
170 }
171
172 Pvals = (m_pressureHBCs[m_intSteps - 1]) + cnt;
173
174 // Getting values on the boundary and filling the pressure bnd
175 // expansion. Multiplication by the normal is required
176 for (int i = 0; i < m_bnd_dim; i++)
177 {
178 m_fields[0]->ExtractElmtToBndPhys(n, Q[i], BndValues[i]);
179 }
180
181 m_PBndExp[n]->NormVectorIProductWRTBase(BndValues, Pvals);
182
183 // Get offset for next terms
184 cnt += ncoeffs;
185 }
186 }
187}
188
189// do nothing unless otherwise defined.
191 [[maybe_unused]] const Array<OneD, NekDouble> &pressure)
192{
193}
194
195// do nothing unless otherwise defined.
197 [[maybe_unused]] const int nbcoeffs, [[maybe_unused]] const int nreg,
198 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &u)
199{
200}
201
203 const Array<OneD, const Array<OneD, NekDouble>> &fields, NekDouble kinvis)
204{
205 if (!m_houtflow.get())
206 {
207 return;
208 }
209
211 size_t cnt = 0;
212
213 // Evaluate robin primitive coefficient here so they can be
214 // updated whem m_int > 1 Currently not using this update
215 // since we only using u^n at outflow instead of BDF rule.
217
218 for (size_t n = 0; n < m_PBndConds.size(); ++n)
219 {
220 if ((m_hbcType[n] == eOBC) || (m_hbcType[n] == eConvectiveOBC))
221 {
222 // Get expansion with element on this boundary
223 m_bndElmtExps[n]->SetWaveSpace(m_fields[0]->GetWaveSpace());
224 int nqb = m_PBndExp[n]->GetTotPoints();
225 int nq = m_bndElmtExps[n]->GetTotPoints();
226
227 // Get velocity and extrapolate
228 for (int i = 0; i < m_curl_dim; i++)
229 {
230 m_fields[0]->ExtractPhysToBndElmt(
231 n, fields[i],
232 m_houtflow->m_outflowVel[cnt][i][m_intSteps - 1]);
233 ExtrapolateArray(m_houtflow->m_outflowVel[cnt][i]);
234 Velocity[i] = m_houtflow->m_outflowVel[cnt][i][m_intSteps - 1];
235 }
236
237 // Homogeneous case needs conversion to physical space
238 if (m_fields[0]->GetWaveSpace())
239 {
240 for (int i = 0; i < m_curl_dim; i++)
241 {
242 m_bndElmtExps[n]->HomogeneousBwdTrans(
243 Velocity[i].size(), Velocity[i], Velocity[i]);
244 }
245 m_bndElmtExps[n]->SetWaveSpace(false);
246 }
247
248 // Get normal vector
250 m_fields[0]->GetBoundaryNormals(n, normals);
251
252 // Calculate n.gradU.n, div(U)
253 Array<OneD, NekDouble> nGradUn(nqb, 0.0);
254 Array<OneD, NekDouble> divU(nqb, 0.0);
256 Array<OneD, NekDouble> bndVal(nqb, 0.0);
257 for (int i = 0; i < m_curl_dim; i++)
258 {
259 grad[i] = Array<OneD, NekDouble>(nq, 0.0);
260 }
261 for (int i = 0; i < m_curl_dim; i++)
262 {
263 if (m_curl_dim == 2)
264 {
265 m_bndElmtExps[n]->PhysDeriv(Velocity[i], grad[0], grad[1]);
266 }
267 else
268 {
269 m_bndElmtExps[n]->PhysDeriv(Velocity[i], grad[0], grad[1],
270 grad[2]);
271 }
272
273 for (int j = 0; j < m_curl_dim; j++)
274 {
275 m_fields[0]->ExtractElmtToBndPhys(n, grad[j], bndVal);
276 // div(U) = gradU_ii
277 if (i == j)
278 {
279 Vmath::Vadd(nqb, divU, 1, bndVal, 1, divU, 1);
280 }
281 // n.gradU.n = gradU_ij n_i n_j
282 Vmath::Vmul(nqb, normals[i], 1, bndVal, 1, bndVal, 1);
283 Vmath::Vvtvp(nqb, normals[j], 1, bndVal, 1, nGradUn, 1,
284 nGradUn, 1);
285 }
286 }
287
288 // Reset WaveSpace in m_bndElmtExp[n] for next time step
289 if (m_fields[0]->GetWaveSpace())
290 {
291 m_bndElmtExps[n]->SetWaveSpace(true);
292 }
293
294 // Obtain u at the boundary
296 for (int i = 0; i < m_curl_dim; i++)
297 {
298 u[i] = Array<OneD, NekDouble>(nqb, 0.0);
299 m_fields[0]->ExtractElmtToBndPhys(n, Velocity[i], u[i]);
300 }
301
302 // Calculate u.n and u^2
303 Array<OneD, NekDouble> un(nqb, 0.0);
304 Array<OneD, NekDouble> u2(nqb, 0.0);
305 for (int i = 0; i < m_curl_dim; i++)
306 {
307 Vmath::Vvtvp(nqb, normals[i], 1, u[i], 1, un, 1, un, 1);
308 Vmath::Vvtvp(nqb, u[i], 1, u[i], 1, u2, 1, u2, 1);
309 }
310
311 // Calculate S_0(u.n) = 0.5*(1-tanh(u.n/*U0*delta))
312 Array<OneD, NekDouble> S0(nqb, 0.0);
313 for (int i = 0; i < nqb; i++)
314 {
315 S0[i] = 0.5 * (1.0 - tanh(un[i] / (m_houtflow->m_U0 *
316 m_houtflow->m_delta)));
317 }
318
319 // Calculate E(n,u) = ((theta+alpha2)*0.5*(u^2)n +
320 // (1-theta+alpha1)*0.5*(n.u)u ) * S_0(u.n)
321 NekDouble k1 =
322 0.5 * (m_houtflow->m_obcTheta + m_houtflow->m_obcAlpha2);
323 NekDouble k2 =
324 0.5 * (1 - m_houtflow->m_obcTheta + m_houtflow->m_obcAlpha1);
325
327 for (int i = 0; i < m_curl_dim; i++)
328 {
329 E[i] = Array<OneD, NekDouble>(nqb, 0.0);
330
331 Vmath::Smul(nqb, k1, u2, 1, E[i], 1);
332 Vmath::Vmul(nqb, E[i], 1, normals[i], 1, E[i], 1);
333 // Use bndVal as a temporary storage
334 Vmath::Smul(nqb, k2, un, 1, bndVal, 1);
335 Vmath::Vvtvp(nqb, u[i], 1, bndVal, 1, E[i], 1, E[i], 1);
336 Vmath::Vmul(nqb, E[i], 1, S0, 1, E[i], 1);
337 }
338
339 // if non-zero forcing is provided we want to subtract
340 // value if we want to interpret values as being the
341 // desired pressure value. This is now precribed from
342 // the velocity forcing to be consistent with the
343 // paper except f_b = -f_b
344
345 // Calculate (E(n,u) + f_b).n
346 Array<OneD, NekDouble> En(nqb, 0.0);
347 for (int i = 0; i < m_bnd_dim; i++)
348 {
349 // Use bndVal as temporary
350 Vmath::Vsub(nqb, E[i], 1,
351 m_houtflow->m_UBndExp[i][n]->GetPhys(), 1, bndVal,
352 1);
353
354 Vmath::Vvtvp(nqb, normals[i], 1, bndVal, 1, En, 1, En, 1);
355 }
356
357 // Calculate pressure bc = kinvis*n.gradU.n - E.n + f_b.n
358 Array<OneD, NekDouble> pbc(nqb, 0.0);
359 Vmath::Svtvm(nqb, kinvis, nGradUn, 1, En, 1, pbc, 1);
360
361 if (m_hbcType[n] == eOBC)
362 {
363
364 if (m_PBndExp[n]->GetWaveSpace())
365 {
366 m_PBndExp[n]->HomogeneousFwdTrans(nqb, pbc, bndVal);
367 m_PBndExp[n]->FwdTrans(bndVal,
368 m_PBndExp[n]->UpdateCoeffs());
369 }
370 else
371 {
372 m_PBndExp[n]->FwdTrans(pbc, m_PBndExp[n]->UpdateCoeffs());
373 }
374 }
375 else if (m_hbcType[n] == eConvectiveOBC) // add outflow values to
376 // calculation from HBC
377 {
378 int nbcoeffs = m_PBndExp[n]->GetNcoeffs();
379 Array<OneD, NekDouble> bndCoeffs(nbcoeffs, 0.0);
380 if (m_PBndExp[n]->GetWaveSpace())
381 {
382 m_PBndExp[n]->HomogeneousFwdTrans(nqb, pbc, bndVal);
383 m_PBndExp[n]->IProductWRTBase(bndVal, bndCoeffs);
384 }
385 else
386 {
387 m_PBndExp[n]->IProductWRTBase(pbc, bndCoeffs);
388 }
389 // Note we have the negative of what is in the Dong paper in
390 // bndVal
391 Vmath::Svtvp(nbcoeffs, m_houtflow->m_pressurePrimCoeff[n],
392 bndCoeffs, 1, m_PBndExp[n]->UpdateCoeffs(), 1,
393 m_PBndExp[n]->UpdateCoeffs(), 1);
394
395 // evaluate u^n at outflow boundary for velocity BC
396 for (int i = 0; i < m_curl_dim; i++)
397 {
398 m_fields[0]->ExtractElmtToBndPhys(
399 n, m_houtflow->m_outflowVel[cnt][i][0],
400 m_houtflow->m_outflowVelBnd[cnt][i][m_intSteps - 1]);
401
402 EvaluateBDFArray(m_houtflow->m_outflowVelBnd[cnt][i]);
403
404 // point u[i] to BDF evalauted value \hat{u}
405 u[i] = m_houtflow->m_outflowVelBnd[cnt][i][m_intSteps - 1];
406 }
407
408 // Add normal velocity if weak pressure
409 // formulation. In this case there is an
410 // additional \int \hat{u}.n ds on the outflow
411 // boundary since we use the inner product wrt
412 // deriv of basis in pressure solve.
413 AddNormVelOnOBC(cnt, n, u);
414 }
415
416 // Calculate velocity boundary conditions
417 if (m_hbcType[n] == eOBC)
418 {
419 // = (pbc n - kinvis divU n)
420 Vmath::Smul(nqb, kinvis, divU, 1, divU, 1);
421 Vmath::Vsub(nqb, pbc, 1, divU, 1, bndVal, 1);
422 }
423 else if (m_hbcType[n] == eConvectiveOBC)
424 {
425 // = (-kinvis divU n)
426 Vmath::Smul(nqb, -1.0 * kinvis, divU, 1, bndVal, 1);
427
428 // pbc needs to be added after pressure solve
429 }
430
431 for (int i = 0; i < m_curl_dim; ++i)
432 {
433 // Reuse divU -> En
434 Vmath::Vvtvp(nqb, normals[i], 1, bndVal, 1, E[i], 1, divU, 1);
435 // - f_b
436 Vmath::Vsub(nqb, divU, 1,
437 m_houtflow->m_UBndExp[i][n]->GetPhys(), 1, divU, 1);
438 // * 1/kinvis
439 Vmath::Smul(nqb, 1.0 / kinvis, divU, 1, divU, 1);
440
441 if (m_hbcType[n] == eConvectiveOBC)
442 {
443 Vmath::Svtvp(nqb, m_houtflow->m_velocityPrimCoeff[i][n],
444 u[i], 1, divU, 1, divU, 1);
445 }
446
447 if (m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
448 {
449 m_houtflow->m_UBndExp[i][n]->HomogeneousFwdTrans(nqb, divU,
450 divU);
451 }
452
453 m_houtflow->m_UBndExp[i][n]->IProductWRTBase(
454 divU, m_houtflow->m_UBndExp[i][n]->UpdateCoeffs());
455 }
456
457 // Get offset for next terms
458 cnt++;
459 }
460 }
461}
462
464{
465 if (!m_houtflow.get())
466 {
467 return;
468 }
469
470 for (size_t n = 0; n < m_PBndConds.size(); ++n)
471 {
472 if (m_hbcType[n] == eConvectiveOBC)
473 {
474 int nqb = m_PBndExp[n]->GetTotPoints();
475 int ncb = m_PBndExp[n]->GetNcoeffs();
476
477 m_pressure->FillBndCondFromField(n, m_pressure->GetCoeffs());
478 Array<OneD, NekDouble> pbc(nqb);
479
480 m_PBndExp[n]->BwdTrans(m_PBndExp[n]->GetCoeffs(), pbc);
481
482 if (m_PBndExp[n]->GetWaveSpace())
483 {
484 m_PBndExp[n]->HomogeneousBwdTrans(nqb, pbc, pbc);
485 }
486
488 Array<OneD, NekDouble> wk1(ncb);
489
490 // Get normal vector
492 m_fields[0]->GetBoundaryNormals(n, normals);
493
494 // Add 1/kinvis * (pbc n )
495 for (int i = 0; i < m_curl_dim; ++i)
496 {
497 Vmath::Vmul(nqb, normals[i], 1, pbc, 1, wk, 1);
498
499 Vmath::Smul(nqb, 1.0 / kinvis, wk, 1, wk, 1);
500
501 if (m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
502 {
503 m_houtflow->m_UBndExp[i][n]->HomogeneousFwdTrans(nqb, wk,
504 wk);
505 }
506 m_houtflow->m_UBndExp[i][n]->IProductWRTBase(wk, wk1);
507
508 Vmath::Vadd(ncb, wk1, 1,
509 m_houtflow->m_UBndExp[i][n]->GetCoeffs(), 1,
510 m_houtflow->m_UBndExp[i][n]->UpdateCoeffs(), 1);
511 }
512 }
513 }
514}
515
517 const Array<OneD, const Array<OneD, NekDouble>> &Vel,
518 Array<OneD, NekDouble> &IProdVn)
519{
520 int i;
521 size_t n, cnt;
522 Array<OneD, NekDouble> IProdVnTmp;
524
525 for (n = cnt = 0; n < m_PBndConds.size(); ++n)
526 {
527 // High order boundary condition;
528 if (m_hbcType[n] == eHBCNeumann)
529 {
530 for (i = 0; i < m_bnd_dim; ++i)
531 {
532 m_fields[0]->ExtractPhysToBnd(n, Vel[i], velbc[i]);
533 }
534 IProdVnTmp = IProdVn + cnt;
535 m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
536 cnt += m_PBndExp[n]->GetNcoeffs();
537 }
538 else if (m_hbcType[n] == eConvectiveOBC) // skip over conective OBC
539 {
540 cnt += m_PBndExp[n]->GetNcoeffs();
541 }
542 }
543}
544
546{
547
548 if (!m_HBCnumber)
549 {
550 return;
551 }
552 int i;
553 size_t n, cnt;
554 Array<OneD, NekDouble> IProdVnTmp;
557 m_bnd_dim);
558 for (i = 0; i < m_bnd_dim; ++i)
559 {
560 VelBndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
561 }
562
563 for (n = cnt = 0; n < m_PBndConds.size(); ++n)
564 {
565 // High order boundary condition;
566 if (m_hbcType[n] == eHBCNeumann)
567 {
568 for (i = 0; i < m_bnd_dim; ++i)
569 {
570 velbc[i] = Array<OneD, NekDouble>(
571 VelBndExp[i][n]->GetTotPoints(), 0.0);
572 VelBndExp[i][n]->SetWaveSpace(
573 m_fields[m_velocity[i]]->GetWaveSpace());
574 VelBndExp[i][n]->BwdTrans(VelBndExp[i][n]->GetCoeffs(),
575 velbc[i]);
576 }
577 IProdVnTmp = IProdVn + cnt;
578 m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
579 cnt += m_PBndExp[n]->GetNcoeffs();
580 }
581 else if (m_hbcType[n] == eConvectiveOBC)
582 {
583 // skip over convective OBC
584 cnt += m_PBndExp[n]->GetNcoeffs();
585 }
586 }
587}
588
589/**
590 * Function to roll time-level storages to the next step layout.
591 * The stored data associated with the oldest time-level
592 * (not required anymore) are moved to the top, where they will
593 * be overwritten as the solution process progresses.
594 */
596{
597 int nlevels = input.size();
598
600
601 tmp = input[nlevels - 1];
602
603 for (int n = nlevels - 1; n > 0; --n)
604 {
605 input[n] = input[n - 1];
606 }
607
608 input[0] = tmp;
609}
610
611/**
612 * Initialise boundary expansion lists for each domain boundary
613 * Each boundary expansion list contains all elements that touch the boundary.
614 * Construct for every boundary and not only higher-order pressure BCs.
615 */
617{
618 size_t n, nBndElmtExp = m_pressure->GetBndConditions().size();
619
620 // Initialise Array of pointers to BndEltmExpansion(-Lists)
622
623 // Loop n domain boundaries and initialise the boundary expansion list
624 for (n = 0; n < nBndElmtExp; ++n)
625 {
626 m_fields[0]->GetBndElmtExpansion(n, m_bndElmtExps[n], false);
627 }
628}
629
630/**
631 * Initialize HOBCs
632 */
635{
636 m_PBndConds = m_pressure->GetBndConditions();
637 m_PBndExp = m_pressure->GetBndCondExpansions();
638
639 size_t cnt, n;
640
641 // Storage array for high order pressure BCs
644
645 // Get useful values for HOBCs
646 m_HBCnumber = 0;
647 m_numHBCDof = 0;
648
649 int outHBCnumber = 0;
650 int numOutHBCPts = 0;
651
653 for (n = 0; n < m_PBndConds.size(); ++n)
654 {
655 // High order boundary Neumann Condiiton
656 if (boost::iequals(m_PBndConds[n]->GetUserDefined(), "H"))
657 {
659 m_numHBCDof += m_PBndExp[n]->GetNcoeffs();
660 m_HBCnumber += m_PBndExp[n]->GetExpSize();
661 }
662
663 // High order outflow convective condition
664 if (m_PBndConds[n]->GetBoundaryConditionType() ==
666 boost::iequals(m_PBndConds[n]->GetUserDefined(), "HOutflow"))
667 {
669 m_numHBCDof += m_PBndExp[n]->GetNcoeffs();
670 m_HBCnumber += m_PBndExp[n]->GetExpSize();
671 numOutHBCPts += m_PBndExp[n]->GetTotPoints();
672 outHBCnumber++;
673 }
674 // High order outflow boundary condition;
675 else if (boost::iequals(m_PBndConds[n]->GetUserDefined(), "HOutflow"))
676 {
677 m_hbcType[n] = eOBC;
678 numOutHBCPts += m_PBndExp[n]->GetTotPoints();
679 outHBCnumber++;
680 }
681 }
682
684 for (int n = 0; n < m_intSteps; ++n)
685 {
688 }
689
690 m_pressureCalls = 0;
691
692 switch (m_pressure->GetExpType())
693 {
695 {
696 m_curl_dim = 2;
697 m_bnd_dim = 2;
698 }
699 break;
701 {
702 m_curl_dim = 3;
703 m_bnd_dim = 2;
704 }
705 break;
707 {
708 m_curl_dim = 3;
709 m_bnd_dim = 1;
710 }
711 break;
713 {
714 m_curl_dim = 3;
715 m_bnd_dim = 3;
716 }
717 break;
718 default:
719 ASSERTL0(0, "Dimension not supported");
720 break;
721 }
722
723 // Initialise storage for outflow HOBCs
724 if (numOutHBCPts > 0)
725 {
727 numOutHBCPts, outHBCnumber, m_curl_dim, pSession);
728
729 // set up boundary expansions link
730 for (int i = 0; i < m_curl_dim; ++i)
731 {
732 m_houtflow->m_UBndExp[i] =
733 m_fields[m_velocity[i]]->GetBndCondExpansions();
734 }
735
736 for (n = 0, cnt = 0; n < m_PBndConds.size(); ++n)
737 {
738 if (boost::iequals(m_PBndConds[n]->GetUserDefined(), "HOutflow"))
739 {
740 m_houtflow->m_outflowVel[cnt] =
742 m_curl_dim);
743
744 m_houtflow->m_outflowVelBnd[cnt] =
746 m_curl_dim);
747
748 int nqb = m_PBndExp[n]->GetTotPoints();
749 int nq = m_bndElmtExps[n]->GetTotPoints();
750 for (int j = 0; j < m_curl_dim; ++j)
751 {
752 m_houtflow->m_outflowVel[cnt][j] =
754
755 m_houtflow->m_outflowVelBnd[cnt][j] =
757
758 for (int k = 0; k < m_intSteps; ++k)
759 {
760 m_houtflow->m_outflowVel[cnt][j][k] =
761 Array<OneD, NekDouble>(nq, 0.0);
762 m_houtflow->m_outflowVelBnd[cnt][j][k] =
763 Array<OneD, NekDouble>(nqb, 0.0);
764 }
765 }
766 cnt++;
767 }
768
769 // evaluate convective primitive coefficient if
770 // convective OBCs are used
771 if (m_hbcType[n] == eConvectiveOBC)
772 {
773 // initialise convective members of
774 // HighOrderOutflow struct
775 if (m_houtflow->m_pressurePrimCoeff.size() == 0)
776 {
777 m_houtflow->m_pressurePrimCoeff =
779 m_houtflow->m_velocityPrimCoeff =
781
782 for (int i = 0; i < m_curl_dim; ++i)
783 {
784 m_houtflow->m_velocityPrimCoeff[i] =
786 }
787 }
788
790 std::static_pointer_cast<
792 ->m_robinPrimitiveCoeff;
793
794 // checkout equation evaluation options!!
795 m_houtflow->m_pressurePrimCoeff[n] = coeff.Evaluate();
796
797 for (int i = 0; i < m_curl_dim; ++i)
798 {
800 UBndConds = m_fields[m_velocity[i]]->GetBndConditions();
801
803 std::static_pointer_cast<
805 UBndConds[n])
806 ->m_robinPrimitiveCoeff;
807
808 m_houtflow->m_defVelPrimCoeff[i] = coeff1.GetExpression();
809
810 ASSERTL1(UBndConds[n]->GetBoundaryConditionType() ==
812 "Require Velocity "
813 "conditions to be of Robin type when pressure"
814 "outflow is specticied as Robin Boundary type");
815
816 // checkout equation evaluation options!!
817 m_houtflow->m_velocityPrimCoeff[i][n] = coeff1.Evaluate();
818 }
819 }
820 }
821 }
822}
823
825{
826
828 {
829 return;
830 }
831
832 for (size_t n = 0; n < m_PBndConds.size(); ++n)
833 {
834 // Get expansion with element on this boundary
835 if (m_hbcType[n] == eConvectiveOBC)
836 {
837 for (int i = 0; i < m_curl_dim; ++i)
838 {
840 m_fields[m_velocity[i]]->UpdateBndConditions();
841
842 std::string primcoeff =
843 m_houtflow->m_defVelPrimCoeff[i] + "*" +
844 boost::lexical_cast<std::string>(
846
847 SpatialDomains::RobinBCShPtr rcond = std::dynamic_pointer_cast<
849
853 m_session, rcond->m_robinFunction.GetExpression(),
854 primcoeff, rcond->GetUserDefined(),
855 rcond->m_filename);
856
857 UBndConds[n] = bcond;
858 }
859 }
860 }
861}
862
863/**
864 *
865 */
867 const Array<OneD, Array<OneD, NekDouble>> inarray)
868{
869 // Checking if the problem is 2D
870 ASSERTL0(m_curl_dim >= 2, "Method not implemented for 1D");
871
872 size_t n_points_0 = m_fields[0]->GetExp(0)->GetTotPoints();
873 size_t n_element = m_fields[0]->GetExpSize();
874 size_t nvel = inarray.size();
875 size_t cnt;
876
877 NekDouble pntVelocity;
878
879 // Getting the standard velocity vector
880 Array<OneD, Array<OneD, NekDouble>> stdVelocity(nvel);
882 Array<OneD, NekDouble> maxV(n_element, 0.0);
884
885 for (size_t i = 0; i < nvel; ++i)
886 {
887 stdVelocity[i] = Array<OneD, NekDouble>(n_points_0);
888 }
889
890 cnt = 0.0;
891 for (size_t el = 0; el < n_element; ++el)
892 {
893 size_t n_points = m_fields[0]->GetExp(el)->GetTotPoints();
894 ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
895
896 // reset local space
897 if (n_points != n_points_0)
898 {
899 for (size_t j = 0; j < nvel; ++j)
900 {
901 stdVelocity[j] = Array<OneD, NekDouble>(n_points, 0.0);
902 }
903 n_points_0 = n_points;
904 }
905 else
906 {
907 for (size_t j = 0; j < nvel; ++j)
908 {
909 Vmath::Zero(n_points, stdVelocity[j], 1);
910 }
911 }
912
914 ->GetExp(el)
915 ->GetGeom()
916 ->GetMetricInfo()
917 ->GetDerivFactors(ptsKeys);
918
919 if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype() ==
921 {
922 for (size_t j = 0; j < nvel; ++j)
923 {
924 for (size_t k = 0; k < nvel; ++k)
925 {
926 Vmath::Vvtvp(n_points, gmat[k * nvel + j], 1,
927 tmp = inarray[k] + cnt, 1, stdVelocity[j], 1,
928 stdVelocity[j], 1);
929 }
930 }
931 }
932 else
933 {
934 for (size_t j = 0; j < nvel; ++j)
935 {
936 for (size_t k = 0; k < nvel; ++k)
937 {
938 Vmath::Svtvp(n_points, gmat[k * nvel + j][0],
939 tmp = inarray[k] + cnt, 1, stdVelocity[j], 1,
940 stdVelocity[j], 1);
941 }
942 }
943 }
944 cnt += n_points;
945
946 // Calculate total velocity in stdVelocity[0]
947 Vmath::Vmul(n_points, stdVelocity[0], 1, stdVelocity[0], 1,
948 stdVelocity[0], 1);
949 for (size_t k = 1; k < nvel; ++k)
950 {
951 Vmath::Vvtvp(n_points, stdVelocity[k], 1, stdVelocity[k], 1,
952 stdVelocity[0], 1, stdVelocity[0], 1);
953 }
954 pntVelocity = Vmath::Vmax(n_points, stdVelocity[0], 1);
955 maxV[el] = sqrt(pntVelocity);
956 }
957
958 return maxV;
959}
960
962{
963 return "";
964}
965
966/**
967 * At the start, the newest value is stored in array[nlevels-1]
968 * and the previous values in the first positions
969 * At the end, the extrapolated value is stored in array[nlevels-1]
970 * and the storage has been updated to included the new value
971 */
973{
974 int nint = std::min(m_pressureCalls, m_intSteps);
975 int nlevels = array.size();
976 int nPts = array[0].size();
977
978 // Check integer for time levels
979 // Note that ExtrapolateArray assumes m_pressureCalls is >= 1
980 // meaning v_EvaluatePressureBCs has been called previously
981 ASSERTL0(nint > 0, "nint must be > 0 when calling ExtrapolateArray.");
982
983 // Update array
984 RollOver(array);
985
986 // Extrapolate to outarray
987 Vmath::Smul(nPts, StifflyStable_Betaq_Coeffs[nint - 1][nint - 1],
988 array[nint - 1], 1, array[nlevels - 1], 1);
989
990 for (int n = 0; n < nint - 1; ++n)
991 {
992 Vmath::Svtvp(nPts, StifflyStable_Betaq_Coeffs[nint - 1][n], array[n], 1,
993 array[nlevels - 1], 1, array[nlevels - 1], 1);
994 }
995}
996
997/**
998 * At the start, the newest value is stored in array[nlevels-1]
999 * and the previous values in the first positions
1000 * At the end, the value of the bdf explicit part is stored in
1001 * array[nlevels-1] and the storage has been updated to included the new value
1002 */
1004{
1005 int nint = std::min(m_pressureCalls, m_intSteps);
1006 int nlevels = array.size();
1007 int nPts = array[0].size();
1008
1009 // Update array
1010 RollOver(array);
1011
1012 // Extrapolate to outarray
1013 Vmath::Smul(nPts, StifflyStable_Alpha_Coeffs[nint - 1][nint - 1],
1014 array[nint - 1], 1, array[nlevels - 1], 1);
1015
1016 for (int n = 0; n < nint - 1; ++n)
1017 {
1018 Vmath::Svtvp(nPts, StifflyStable_Alpha_Coeffs[nint - 1][n], array[n], 1,
1019 array[nlevels - 1], 1, array[nlevels - 1], 1);
1020 }
1021}
1022
1023/**
1024 * At the start, the newest value is stored in array[nlevels-1]
1025 * and the previous values in the first positions
1026 * At the end, the acceleration from BDF is stored in array[nlevels-1]
1027 * and the storage has been updated to included the new value
1028 */
1030{
1031 int nlevels = array.size();
1032 int nPts = array[0].size();
1033
1034 if (nPts)
1035 {
1036 // Update array
1037 RollOver(array);
1038
1039 // Calculate acceleration using Backward Differentiation Formula
1040 Array<OneD, NekDouble> accelerationTerm(nPts, 0.0);
1041 if (m_pressureCalls > 2)
1042 {
1043 int acc_order = std::min(m_pressureCalls - 2, m_intSteps);
1044 Vmath::Smul(nPts, StifflyStable_Gamma0_Coeffs[acc_order - 1],
1045 array[0], 1, accelerationTerm, 1);
1046
1047 for (int i = 0; i < acc_order; i++)
1048 {
1050 nPts, -1 * StifflyStable_Alpha_Coeffs[acc_order - 1][i],
1051 array[i + 1], 1, accelerationTerm, 1, accelerationTerm, 1);
1052 }
1053 }
1054 array[nlevels - 1] = accelerationTerm;
1055 }
1056}
1057
1059{
1060 size_t n, cnt;
1061 for (cnt = n = 0; n < m_PBndConds.size(); ++n)
1062 {
1063 if ((m_hbcType[n] == eHBCNeumann) || (m_hbcType[n] == eConvectiveOBC))
1064 {
1065 int nq = m_PBndExp[n]->GetNcoeffs();
1066 Vmath::Vcopy(nq, &(m_pressureHBCs[m_intSteps - 1])[cnt], 1,
1067 &(m_PBndExp[n]->UpdateCoeffs()[0]), 1);
1068 cnt += nq;
1069 }
1070 }
1071}
1072
1073} // 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, NekDouble > > m_pressureHBCs
Storage for current and previous levels of high order pressure boundary conditions.
Definition: Extrapolate.h:277
virtual void v_AddNormVelOnOBC(const int nbcoeffs, const int nreg, Array< OneD, Array< OneD, NekDouble > > &u)
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:250
Extrapolate(const LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields, MultiRegions::ExpListSharedPtr pPressure, const Array< OneD, int > pVel, const SolverUtils::AdvectionSharedPtr advObject)
Definition: Extrapolate.cpp:53
int m_curl_dim
Curl-curl dimensionality.
Definition: Extrapolate.h:247
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:234
static std::string def
Definition: Extrapolate.h:294
static NekDouble StifflyStable_Betaq_Coeffs[3][3]
Definition: Extrapolate.h:286
virtual void v_AccelerationBDF(Array< OneD, Array< OneD, NekDouble > > &array)
virtual void v_CalcNeumannPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
void CopyPressureHBCsToPbndExp(void)
virtual void v_CorrectPressureBCs(const Array< OneD, NekDouble > &pressure)
void EvaluateBDFArray(Array< OneD, Array< OneD, NekDouble > > &array)
void IProductNormVelocityBCOnHBC(Array< OneD, NekDouble > &IprodVn)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Velocity fields.
Definition: Extrapolate.h:231
Array< OneD, Array< OneD, NekDouble > > m_iprodnormvel
Storage for current and previous levels of the inner product of normal velocity.
Definition: Extrapolate.h:281
Array< OneD, NekDouble > GetMaxStdVelocity(const Array< OneD, Array< OneD, NekDouble > > inarray)
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:271
void GenerateHOPBCMap(const LibUtilities::SessionReaderSharedPtr &pSsession)
void IProductNormVelocityOnHBC(const Array< OneD, const Array< OneD, NekDouble > > &Vel, Array< OneD, NekDouble > &IprodVn)
Array< OneD, HBCType > m_hbcType
Array of type of high order BCs for splitting shemes.
Definition: Extrapolate.h:228
HighOrderOutflowSharedPtr m_houtflow
Definition: Extrapolate.h:291
NekDouble m_timestep
Definition: Extrapolate.h:273
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:256
void ExtrapolateArray(Array< OneD, Array< OneD, NekDouble > > &array)
void GenerateBndElmtExpansion(void)
void MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
Definition: Extrapolate.h:107
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition: Extrapolate.h:238
static NekDouble StifflyStable_Alpha_Coeffs[3][3]
Definition: Extrapolate.h:287
void AddNormVelOnOBC(const int nbcoeffs, const int nreg, Array< OneD, Array< OneD, NekDouble > > &u)
Definition: Extrapolate.h:121
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndElmtExps
Boundary expansions on each domain boundary.
Definition: Extrapolate.h:259
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:253
int m_pressureCalls
number of times the high-order pressure BCs have been called
Definition: Extrapolate.h:262
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:223
void UpdateRobinPrimCoeff(void)
void CalcOutflowBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, NekDouble kinvis)
void AddPressureToOutflowBCs(NekDouble kinvis)
LibUtilities::CommSharedPtr m_comm
Definition: Extrapolate.h:225
static NekDouble StifflyStable_Gamma0_Coeffs[3]
Definition: Extrapolate.h:288
virtual std::string v_GetSubStepName(void)
std::string GetExpression(void) const
Provides a generic Factory class.
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
An abstract base class encapsulating the concept of advection of a vector field.
Definition: Advection.h:81
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:54
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:213
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:216
@ eNOHBC
Definition: Extrapolate.h:52
@ eConvectiveOBC
Definition: Extrapolate.h:55
@ eHBCNeumann
Definition: Extrapolate.h:53
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:47
double NekDouble
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 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.hpp:396
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 Svtvm(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvm (scalar times vector minus vector): z = alpha*x - y.
Definition: Vmath.hpp:424
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285