Nektar++
CoupledLinearNS.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: CoupledLinearNS.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: Coupled Solver for the Linearised Incompressible
32// Navier Stokes equations
33///////////////////////////////////////////////////////////////////////////////
34
35#include <boost/algorithm/string.hpp>
36
42
43namespace Nektar
44{
45
48 "CoupledLinearisedNS", CoupledLinearNS::create);
49
52 "SolverType", "CoupledLinearisedNS", eCoupledLinearisedNS);
53
54/**
55 * @class CoupledLinearNS
56 *
57 * Set up expansion field for velocity and pressure, the local to
58 * global mapping arrays and the basic memory definitions for
59 * coupled matrix system
60 */
64 : UnsteadySystem(pSession, pGraph), IncNavierStokes(pSession, pGraph),
65 m_zeroMode(false)
66{
67}
68
69void CoupledLinearNS::v_InitObject(bool DeclareField)
70{
72
73 size_t i;
74 int expdim = m_graph->GetMeshDimension();
75
76 // Get Expansion list for orthogonal expansion at p-2
77 const SpatialDomains::ExpansionInfoMap &pressure_exp =
78 GenPressureExp(m_graph->GetExpansionInfo("u"));
79
81 if (boost::iequals(
82 m_boundaryConditions->GetVariable(m_nConvectiveFields - 1), "p"))
83 {
84 ASSERTL0(false, "Last field is defined as pressure but this is not "
85 "suitable for this solver, please remove this field as "
86 "it is implicitly defined");
87 }
88 // Decide how to declare explist for pressure.
89 if (expdim == 2)
90 {
91 size_t nz;
92
94 {
95 ASSERTL0(m_fields.size() > 2, "Expect to have three at least three "
96 "components of velocity variables");
97 LibUtilities::BasisKey Homo1DKey =
98 m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
99
102 m_homogen_dealiasing, pressure_exp);
103
104 ASSERTL1(m_npointsZ % 2 == 0,
105 "Non binary number of planes have been specified");
106 nz = m_npointsZ / 2;
107 }
108 else
109 {
110 m_pressure =
112 m_session, pressure_exp);
113 nz = 1;
114 }
115
117 for (i = 0; i < m_velocity.size(); ++i)
118 {
119 velocity[i] = m_fields[m_velocity[i]];
120 }
121
122 // Set up Array of mappings
124
125 if (m_singleMode)
126 {
127
128 ASSERTL0(nz <= 2,
129 "For single mode calculation can only have nz <= 2");
130 if (m_session->DefinesSolverInfo("BetaZero"))
131 {
132 m_zeroMode = true;
133 }
134 int nz_loc = 2;
135 m_locToGloMap[0] =
138 m_pressure, nz_loc, m_zeroMode);
139 }
140 else
141 {
142 // base mode
143 int nz_loc = 1;
144 m_locToGloMap[0] =
147 m_pressure, nz_loc);
148
149 if (nz > 1)
150 {
151 nz_loc = 2;
152 // Assume all higher modes have the same boundary conditions and
153 // re-use mapping
154 m_locToGloMap[1] =
158 m_pressure->GetPlane(2), nz_loc, false);
159 // Note high order modes cannot be singular
160 for (i = 2; i < nz; ++i)
161 {
163 }
164 }
165 }
166 }
167 else if (expdim == 3)
168 {
169 // m_pressure =
170 // MemoryManager<MultiRegions::ExpList3D>::AllocateSharedPtr(pressure_exp);
171 ASSERTL0(false, "Setup mapping aray");
172 }
173 else
174 {
175 ASSERTL0(false, "Exp dimension not recognised");
176 }
177
178 // creation of the extrapolation object
180 {
181 std::string vExtrapolation = "Standard";
182
183 if (m_session->DefinesSolverInfo("Extrapolation"))
184 {
185 vExtrapolation = m_session->GetSolverInfo("Extrapolation");
186 }
187
189 vExtrapolation, m_session, m_fields, m_pressure, m_velocity,
191 }
192}
193
194/**
195 * Set up a coupled linearised Naviers-Stokes solve. Primarily a
196 * wrapper fuction around the full mode by mode version of
197 * SetUpCoupledMatrix.
198 *
199 */
201 const NekDouble lambda, const Array<OneD, Array<OneD, NekDouble>> &Advfield,
202 bool IsLinearNSEquation)
203{
204
205 int nz;
206 if (m_singleMode)
207 {
208
209 NekDouble lambda_imag;
210
211 // load imaginary component of any potential shift
212 // Probably should be called from DriverArpack but not yet
213 // clear how to do this
214 m_session->LoadParameter("imagShift", lambda_imag,
216 nz = 1;
218
219 ASSERTL1(m_npointsZ <= 2, "Only expected a maxmimum of two planes in "
220 "single mode linear NS solver");
221
222 if (m_zeroMode)
223 {
224 SetUpCoupledMatrix(lambda, Advfield, IsLinearNSEquation, 0,
225 m_mat[0], m_locToGloMap[0], lambda_imag);
226 }
227 else
228 {
229 NekDouble beta = 2 * M_PI / m_LhomZ;
230 NekDouble lam = lambda + m_kinvis * beta * beta;
231
232 SetUpCoupledMatrix(lam, Advfield, IsLinearNSEquation, 1, m_mat[0],
233 m_locToGloMap[0], lambda_imag);
234 }
235 }
236 else
237 {
238 int n;
239 if (m_npointsZ > 1)
240 {
241 nz = m_npointsZ / 2;
242 }
243 else
244 {
245 nz = 1;
246 }
247
249
250 // mean mode or 2D mode.
251 SetUpCoupledMatrix(lambda, Advfield, IsLinearNSEquation, 0, m_mat[0],
252 m_locToGloMap[0]);
253
254 for (n = 1; n < nz; ++n)
255 {
256 NekDouble beta = 2 * M_PI * n / m_LhomZ;
257
258 NekDouble lam = lambda + m_kinvis * beta * beta;
259
260 SetUpCoupledMatrix(lam, Advfield, IsLinearNSEquation, n, m_mat[n],
261 m_locToGloMap[n]);
262 }
263 }
264}
265
266/**
267 * Set up a coupled linearised Naviers-Stokes solve in the
268 * following manner:
269 *
270 * Consider the linearised Navier-Stokes element matrix denoted as
271 *
272 * \f[ \left [ \begin{array}{ccc} A
273 * & -D_{bnd}^T & B\\ -D_{bnd}& 0 & -D_{int}\\ \tilde{B}^T &
274 * -D_{int}^T & C \end{array} \right ] \left [ \begin{array}{c} {\bf
275 * v}_{bnd} \\ p\\ {\bf v}_{int} \end{array} \right ] = \left [
276 * \begin{array}{c} {\bf f}_{bnd} \\ 0\\ {\bf f}_{int} \end{array}
277 * \right ] \f]
278 *
279 * where \f${\bf v}_{bnd}, {\bf f}_{bnd}\f$ denote the degrees of
280 * freedom of the elemental velocities on the boundary of the
281 * element, \f${\bf v}_{int}, {\bf f}_{int}\f$ denote the degrees
282 * of freedom of the elemental velocities on the interior of the
283 * element and \f$p\f$ is the piecewise continuous pressure. The
284 * matrices have the interpretation
285 *
286 * \f[ A[n,m] = (\nabla \phi^b_n, \nu \nabla
287 * \phi^b_m) + (\phi^b_n,{\bf U \cdot \nabla} \phi^b_m) +
288 * (\phi^b_n \nabla^T {\bf U} \phi^b_m), \f]
289 *
290 * \f[ B[n,m] = (\nabla \phi^b_n, \nu \nabla \phi^i_m) +
291 * (\phi^b_n,{\bf U \cdot \nabla} \phi^i_m) + (\phi^b_n \nabla^T
292 * {\bf U} \phi^i_m),\f]
293 *
294 * \f[\tilde{B}^T[n,m] = (\nabla \phi^i_n, \nu \nabla \phi^b_m) +
295 * (\phi^i_n,{\bf U \cdot \nabla} \phi^b_m) + (\phi^i_n \nabla^T
296 * {\bf U} \phi^b_m) \f]
297 *
298 * \f[ C[n,m] = (\nabla \phi^i_n, \nu \nabla
299 * \phi^i_m) + (\phi^i_n,{\bf U \cdot \nabla} \phi^i_m) +
300 * (\phi^i_n \nabla^T {\bf U} \phi^i_m),\f]
301 *
302 * \f[ D_{bnd}[n,m] = (\psi_m,\nabla \phi^b),\f]
303 *
304 * \f[ D_{int}[n,m] = (\psi_m,\nabla \phi^i) \f]
305 *
306 * where \f$\psi\f$ is the space of pressures typically at order
307 * \f$P-2\f$ and \f$\phi\f$ is the velocity vector space of
308 * polynomial order \f$P\f$. (Note the above definitions for the
309 * transpose terms shoudl be interpreted with care since we have
310 * used a tensor differential for the \f$\nabla^T \f$ terms)
311 *
312 * Note \f$B = \tilde{B}^T\f$ if just considering the stokes
313 * operator and then \f$C\f$ is also symmetric. Also note that
314 * \f$A,B\f$ and \f$C\f$ are block diagonal in the Oseen equation
315 * when \f$\nabla^T {\bf U}\f$ are zero.
316 *
317 * Since \f$C\f$ is invertible we can premultiply the governing
318 * elemental equation by the following matrix:
319 *
320 * \f[ \left [ \begin{array}{ccc} I & 0 &
321 * -BC^{-1}\\ 0& I & D_{int}C^{-1}\\ 0 & 0 & I \end{array}
322 * \right ] \left \{ \left [ \begin{array}{ccc} A & -D_{bnd}^T &
323 * B\\ -D_{bnd}& 0 & -D_{int}\\ \tilde{B}^T & -D_{int}^T & C
324 * \end{array} \right ] \left [ \begin{array}{c} {\bf v}_{bnd} \\
325 * p\\ {\bf v_{int}} \end{array} \right ] = \left [
326 * \begin{array}{c} {\bf f}_{bnd} \\ 0\\ {\bf f_{int}} \end{array}
327 * \right ] \right \} \f]
328 *
329 * which if we multiply out the matrix equation we obtain
330 *
331 * \f[ \left [ \begin{array}{ccc} A - B
332 * C^{-1}\tilde{B}^T & -D_{bnd}^T +B C^{-1} D_{int}^T& 0\\
333 * -D_{bnd}+D_{int} C^{-1} \tilde{B}^T & -D_{int} C^{-1}
334 * D_{int}^T & 0\\ \tilde{B}^T & -D_{int}^T & C \end{array} \right ]
335 * \left [ \begin{array}{c} {\bf v}_{bnd} \\ p\\ {\bf v_{int}}
336 * \end{array} \right ] = \left [ \begin{array}{c} {\bf f}_{bnd}
337 * -B C^{-1} {\bf f}_{int}\\ f_p = D_{int} C^{-1} {\bf
338 * f}_{int}\\ {\bf f_{int}} \end{array} \right ] \f]
339 *
340 * In the above equation the \f${\bf v}_{int}\f$ degrees of freedom
341 * are decoupled and so we need to solve for the \f${\bf v}_{bnd},
342 * p\f$ degrees of freedom. The final step is to perform a second
343 * level of static condensation but where we will lump the mean
344 * pressure mode (or a pressure degree of freedom containing a
345 * mean component) with the velocity boundary degrees of
346 * freedom. To do we define \f${\bf b} = [{\bf v}_{bnd}, p_0]\f$ where
347 * \f$p_0\f$ is the mean pressure mode and \f$\hat{p}\f$ to be the
348 * remainder of the pressure space. We now repartition the top \f$2
349 * \times 2\f$ block of matrices of previous matrix equation as
350 *
351 * \f[ \left [ \begin{array}{cc} \hat{A} & \hat{B}\\ \hat{C} &
352 * \hat{D} \end{array} \right ] \left [ \begin{array}{c} {\bf b}
353 * \\ \hat{p} \end{array} \right ] = \left [ \begin{array}{c}
354 * \hat{\bf f}_{bnd} \\ \hat{f}_p \end{array} \right ]
355 * \label{eqn.linNS_fac2} \f]
356 *
357 * where
358 *
359 * \f[ \hat{A}_{ij} = \left [ \begin{array}{cc} A - B
360 * C^{-1}\tilde{B}^T & [-D_{bnd}^T +B C^{-1} D_{int}^T]_{i0}\\
361 * {[}-D_{bnd}+D_{int} C^{-1} \tilde{B}^T]_{0j} & -[D_{int}
362 * C^{-1} D_{int}^T ]_{00} \end{array} \right ] \f]
363 *
364 * \f[ \hat{B}_{ij} = \left [ \begin{array}{c} [-D_{bnd}^T +B
365 * C^{-1} D_{int}^T]_{i,j+1} \\ {[} -D_{int} C^{-1} D^T_{int}
366 * ]_{0j}\end{array} \right ] \f]
367 *
368 * \f[ \hat{C}_{ij} = \left [\begin{array}{cc} -D_{bnd} +
369 * D_{int} C^{-1} \tilde{B}^T, & {[} -D_{int} C^{-1} D^T_{int}
370 * ]_{i+1,0}\end{array} \right ] \f]
371 *
372 * \f[ \hat{D}_{ij} = \left [\begin{array}{c} {[} -D_{int}
373 * C^{-1} D^T_{int} ]_{i+1,j+1}\end{array} \right ] \f]
374 *
375 * and
376 *
377 * \f[ fh\_{bnd} = \left [ \begin{array}{c} {\bf
378 * f}_{bnd} -B C^{-1} {\bf f}_{int}\\ {[}D_{int} C^{-1} {\bf
379 * f}_{int}]_0 \end{array}\right ] \hspace{1cm} [fh\_p_{i} =
380 * \left [ \begin{array}{c} {[}D_{int} C^{-1} {\bf f}_{iint}]_{i+1}
381 * \end{array}\right ] \f]
382 *
383 * Since the \f$\hat{D}\f$ is decoupled and invertible we can now
384 * statically condense the previous matrix equationto decouple
385 * \f${\bf b}\f$ from \f$\hat{p}\f$ by solving the following system
386 *
387 * \f[ \left [ \begin{array}{cc} \hat{A} - \hat{B} \hat{D}^{-1}
388 * \hat{C} & 0 \\ \hat{C} & \hat{D} \end{array} \right ] \left
389 * [ \begin{array}{c} {\bf b} \\ \hat{p} \end{array} \right ] =
390 * \left [ \begin{array}{c} \hat{\bf f}_{bnd} - \hat{B}
391 * \hat{D}^{-1} \hat{f}_p\\ \hat{f}_p \end{array} \right ] \f]
392 *
393 * The matrix \f$\hat{A} - \hat{B} \hat{D}^{-1} \hat{C}\f$ has
394 * to be globally assembled and solved iteratively or
395 * directly. One we obtain the solution to \f${\bf b}\f$ we can use
396 * the second row of equation fourth matrix equation to solve for
397 * \f$\hat{p}\f$ and finally the last row of equation second
398 * matrix equation to solve for \f${\bf v}_{int}\f$.
399 *
400 */
401
403 const NekDouble lambda, const Array<OneD, Array<OneD, NekDouble>> &Advfield,
404 bool IsLinearNSEquation, const int HomogeneousMode,
407 const NekDouble lambda_imag)
408{
409 size_t n, i, j, k;
410 size_t nel = m_fields[m_velocity[0]]->GetNumElmts();
411 size_t nvel = m_velocity.size();
412
413 // if Advfield is defined can assume it is an Oseen or LinearNS equation
414 bool AddAdvectionTerms =
415 (Advfield == NullNekDoubleArrayOfArray) ? false : true;
416 // bool AddAdvectionTerms = true; // Temporary debugging trip
417
418 // call velocity space Static condensation and fill block
419 // matrices. Need to set this up differently for Oseen and
420 // Lin NS. Ideally should make block diagonal for Stokes and
421 // Oseen problems.
422 DNekScalMatSharedPtr loc_mat;
424 NekDouble one = 1.0;
425 size_t nint, nbndry;
426 size_t rows, cols;
427 NekDouble zero = 0.0;
428 Array<OneD, unsigned int> bmap, imap;
429
430 Array<OneD, unsigned int> nsize_bndry(nel);
431 Array<OneD, unsigned int> nsize_bndry_p1(nel);
432 Array<OneD, unsigned int> nsize_int(nel);
433 Array<OneD, unsigned int> nsize_p(nel);
434 Array<OneD, unsigned int> nsize_p_m1(nel);
435
436 size_t nz_loc;
437
438 if (HomogeneousMode) // Homogeneous mode flag
439 {
440 nz_loc = 2;
441 }
442 else
443 {
444 if (m_singleMode)
445 {
446 nz_loc = 2;
447 }
448 else
449 {
450 nz_loc = 1;
451 }
452 }
453
454 // Set up block matrix sizes -
455 for (n = 0; n < nel; ++n)
456 {
457 nsize_bndry[n] = nvel *
458 m_fields[m_velocity[0]]->GetExp(n)->NumBndryCoeffs() *
459 nz_loc;
460 nsize_bndry_p1[n] = nsize_bndry[n] + nz_loc;
461 nsize_int[n] =
462 (nvel * m_fields[m_velocity[0]]->GetExp(n)->GetNcoeffs() * nz_loc -
463 nsize_bndry[n]);
464 nsize_p[n] = m_pressure->GetExp(n)->GetNcoeffs() * nz_loc;
465 nsize_p_m1[n] = nsize_p[n] - nz_loc;
466 }
467
468 MatrixStorage blkmatStorage = eDIAGONAL;
471 nsize_bndry_p1, nsize_bndry_p1, blkmatStorage);
473 nsize_bndry, nsize_int, blkmatStorage);
475 nsize_bndry, nsize_int, blkmatStorage);
477 nsize_int, nsize_int, blkmatStorage);
478
480 nsize_p, nsize_bndry, blkmatStorage);
481
483 nsize_p, nsize_int, blkmatStorage);
484
485 // Final level static condensation matrices.
488 nsize_bndry_p1, nsize_p_m1, blkmatStorage);
491 nsize_p_m1, nsize_bndry_p1, blkmatStorage);
494 blkmatStorage);
495
497 timer.Start();
498
499 for (n = 0; n < nel; ++n)
500 {
501 nbndry = nsize_bndry[n];
502 nint = nsize_int[n];
503 k = nsize_bndry_p1[n];
506 Array<OneD, NekDouble> Ah_data = Ah->GetPtr();
507 int AhRows = k;
510 Array<OneD, NekDouble> B_data = B->GetPtr();
513 Array<OneD, NekDouble> C_data = C->GetPtr();
516 Array<OneD, NekDouble> D_data = D->GetPtr();
517
519 nsize_p[n], nsize_bndry[n], zero);
521 nsize_p[n], nsize_int[n], zero);
522
523 locExp = m_fields[m_velocity[0]]->GetExp(n);
524 locExp->GetBoundaryMap(bmap);
525 locExp->GetInteriorMap(imap);
529 StdRegions::eHelmholtz, locExp->DetShapeType(), *locExp, factors);
530
531 size_t ncoeffs = m_fields[m_velocity[0]]->GetExp(n)->GetNcoeffs();
532 size_t nphys = m_fields[m_velocity[0]]->GetExp(n)->GetTotPoints();
533 size_t nbmap = bmap.size();
534 size_t nimap = imap.size();
535
536 Array<OneD, NekDouble> coeffs(ncoeffs);
537 Array<OneD, NekDouble> phys(nphys);
538 size_t psize = m_pressure->GetExp(n)->GetNcoeffs();
539 size_t pqsize = m_pressure->GetExp(n)->GetTotPoints();
540
541 Array<OneD, NekDouble> deriv(pqsize);
542 Array<OneD, NekDouble> pcoeffs(psize);
543
544 if (AddAdvectionTerms == false) // use static condensed managed matrices
545 {
546 // construct velocity matrices using statically
547 // condensed elemental matrices and then construct
548 // pressure matrix systems
550 CondMat = locExp->GetLocStaticCondMatrix(helmkey);
551
552 for (k = 0; k < nvel * nz_loc; ++k)
553 {
554 DNekScalMat &SubBlock = *CondMat->GetBlock(0, 0);
555 rows = SubBlock.GetRows();
556 cols = SubBlock.GetColumns();
557 for (i = 0; i < rows; ++i)
558 {
559 for (j = 0; j < cols; ++j)
560 {
561 (*Ah)(i + k * rows, j + k * cols) =
562 m_kinvis * SubBlock(i, j);
563 }
564 }
565 }
566
567 for (k = 0; k < nvel * nz_loc; ++k)
568 {
569 DNekScalMat &SubBlock = *CondMat->GetBlock(0, 1);
570 DNekScalMat &SubBlock1 = *CondMat->GetBlock(1, 0);
571 rows = SubBlock.GetRows();
572 cols = SubBlock.GetColumns();
573 for (i = 0; i < rows; ++i)
574 {
575 for (j = 0; j < cols; ++j)
576 {
577 (*B)(i + k * rows, j + k * cols) = SubBlock(i, j);
578 (*C)(i + k * rows, j + k * cols) =
579 m_kinvis * SubBlock1(j, i);
580 }
581 }
582 }
583
584 for (k = 0; k < nvel * nz_loc; ++k)
585 {
586 DNekScalMat &SubBlock = *CondMat->GetBlock(1, 1);
587 NekDouble inv_kinvis = 1.0 / m_kinvis;
588 rows = SubBlock.GetRows();
589 cols = SubBlock.GetColumns();
590 for (i = 0; i < rows; ++i)
591 {
592 for (j = 0; j < cols; ++j)
593 {
594 (*D)(i + k * rows, j + k * cols) =
595 inv_kinvis * SubBlock(i, j);
596 }
597 }
598 }
599
600 // Loop over pressure space and construct boundary block matrices.
601 for (i = 0; i < bmap.size(); ++i)
602 {
603 // Fill element with mode
604 Vmath::Zero(ncoeffs, coeffs, 1);
605 coeffs[bmap[i]] = 1.0;
606 m_fields[m_velocity[0]]->GetExp(n)->BwdTrans(coeffs, phys);
607
608 // Differentiation & Inner product wrt base.
609 for (j = 0; j < nvel; ++j)
610 {
611 if ((nz_loc == 2) && (j == 2)) // handle d/dz derivative
612 {
613 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
614
616 m_fields[m_velocity[0]]->GetExp(n)->GetTotPoints(),
617 beta, phys, 1, deriv, 1);
618
619 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
620
621 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
622 Dbnd->GetRawPtr() +
623 ((nz_loc * j + 1) * bmap.size() + i) *
624 nsize_p[n],
625 1);
626
627 Vmath::Neg(psize, pcoeffs, 1);
628 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
629 Dbnd->GetRawPtr() +
630 ((nz_loc * j) * bmap.size() + i) *
631 nsize_p[n] +
632 psize,
633 1);
634 }
635 else
636 {
637 if (j < 2) // required for mean mode of homogeneous
638 // expansion
639 {
640 locExp->PhysDeriv(MultiRegions::DirCartesianMap[j],
641 phys, deriv);
642 m_pressure->GetExp(n)->IProductWRTBase(deriv,
643 pcoeffs);
644 // copy into column major storage.
645 for (k = 0; k < nz_loc; ++k)
646 {
648 psize, &(pcoeffs)[0], 1,
649 Dbnd->GetRawPtr() +
650 ((nz_loc * j + k) * bmap.size() + i) *
651 nsize_p[n] +
652 k * psize,
653 1);
654 }
655 }
656 }
657 }
658 }
659
660 for (i = 0; i < imap.size(); ++i)
661 {
662 // Fill element with mode
663 Vmath::Zero(ncoeffs, coeffs, 1);
664 coeffs[imap[i]] = 1.0;
665 m_fields[m_velocity[0]]->GetExp(n)->BwdTrans(coeffs, phys);
666
667 // Differentiation & Inner product wrt base.
668 for (j = 0; j < nvel; ++j)
669 {
670 if ((nz_loc == 2) && (j == 2)) // handle d/dz derivative
671 {
672 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
673
675 m_fields[m_velocity[0]]->GetExp(n)->GetTotPoints(),
676 beta, phys, 1, deriv, 1);
677
678 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
679
680 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
681 Dint->GetRawPtr() +
682 ((nz_loc * j + 1) * imap.size() + i) *
683 nsize_p[n],
684 1);
685
686 Vmath::Neg(psize, pcoeffs, 1);
687 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
688 Dint->GetRawPtr() +
689 ((nz_loc * j) * imap.size() + i) *
690 nsize_p[n] +
691 psize,
692 1);
693 }
694 else
695 {
696 if (j < 2) // required for mean mode of homogeneous
697 // expansion
698 {
699 // m_fields[m_velocity[0]]->GetExp(n)->PhysDeriv(j,phys,
700 // deriv);
701 locExp->PhysDeriv(MultiRegions::DirCartesianMap[j],
702 phys, deriv);
703
704 m_pressure->GetExp(n)->IProductWRTBase(deriv,
705 pcoeffs);
706
707 // copy into column major storage.
708 for (k = 0; k < nz_loc; ++k)
709 {
711 psize, &(pcoeffs)[0], 1,
712 Dint->GetRawPtr() +
713 ((nz_loc * j + k) * imap.size() + i) *
714 nsize_p[n] +
715 k * psize,
716 1);
717 }
718 }
719 }
720 }
721 }
722 }
723 else
724 {
725 // construct velocity matrices and pressure systems at
726 // the same time resusing differential of velocity
727 // space
728
729 DNekScalMat &HelmMat =
730 *(locExp->as<LocalRegions::Expansion>()->GetLocMatrix(helmkey));
731 DNekScalMatSharedPtr MassMat;
732
733 Array<OneD, const NekDouble> HelmMat_data =
734 HelmMat.GetOwnedMatrix()->GetPtr();
735 NekDouble HelmMatScale = HelmMat.Scale();
736 int HelmMatRows = HelmMat.GetRows();
737
738 if ((lambda_imag != NekConstants::kNekUnsetDouble) && (nz_loc == 2))
739 {
741 StdRegions::eMass, locExp->DetShapeType(), *locExp);
742 MassMat = locExp->as<LocalRegions::Expansion>()->GetLocMatrix(
743 masskey);
744 }
745
747 Array<OneD, Array<OneD, NekDouble>> AdvDeriv(nvel * nvel);
748 // Use ExpList phys array for temporaary storage
749 Array<OneD, NekDouble> tmpphys = m_fields[0]->UpdatePhys();
750 int phys_offset = m_fields[m_velocity[0]]->GetPhys_Offset(n);
751 size_t nv;
752 int npoints = locExp->GetTotPoints();
753
754 // Calculate derivative of base flow
755 if (IsLinearNSEquation)
756 {
757 size_t nv1;
758 size_t cnt = 0;
759 AdvDeriv[0] = Array<OneD, NekDouble>(nvel * nvel * npoints);
760 for (nv = 0; nv < nvel; ++nv)
761 {
762 for (nv1 = 0; nv1 < nvel; ++nv1)
763 {
764 if (cnt < nvel * nvel - 1)
765 {
766 AdvDeriv[cnt + 1] = AdvDeriv[cnt] + npoints;
767 ++cnt;
768 }
769
770 if ((nv1 == 2) && (m_HomogeneousType == eHomogeneous1D))
771 {
772 Vmath::Zero(npoints, AdvDeriv[nv * nvel + nv1],
773 1); // dU/dz = 0
774 }
775 else
776 {
777 locExp->PhysDeriv(
779 Advfield[nv] + phys_offset,
780 AdvDeriv[nv * nvel + nv1]);
781 }
782 }
783 }
784 }
785
786 for (i = 0; i < nbmap; ++i)
787 {
788
789 // Fill element with mode
790 Vmath::Zero(ncoeffs, coeffs, 1);
791 coeffs[bmap[i]] = 1.0;
792 locExp->BwdTrans(coeffs, phys);
793
794 for (k = 0; k < nvel * nz_loc; ++k)
795 {
796 for (j = 0; j < nbmap; ++j)
797 {
798 // Ah_data[i+k*nbmap +
799 // (j+k*nbmap)*AhRows] +=
800 // m_kinvis*HelmMat(bmap[i],bmap[j]);
801 Ah_data[i + k * nbmap + (j + k * nbmap) * AhRows] +=
802 m_kinvis * HelmMatScale *
803 HelmMat_data[bmap[i] + HelmMatRows * bmap[j]];
804 }
805
806 for (j = 0; j < nimap; ++j)
807 {
808 B_data[i + k * nbmap + (j + k * nimap) * nbndry] +=
809 m_kinvis * HelmMatScale *
810 HelmMat_data[bmap[i] + HelmMatRows * imap[j]];
811 }
812 }
813
814 if ((lambda_imag != NekConstants::kNekUnsetDouble) &&
815 (nz_loc == 2))
816 {
817 for (k = 0; k < nvel; ++k)
818 {
819 for (j = 0; j < nbmap; ++j)
820 {
821 Ah_data[i + 2 * k * nbmap +
822 (j + (2 * k + 1) * nbmap) * AhRows] -=
823 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
824 }
825
826 for (j = 0; j < nbmap; ++j)
827 {
828 Ah_data[i + (2 * k + 1) * nbmap +
829 (j + 2 * k * nbmap) * AhRows] +=
830 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
831 }
832
833 for (j = 0; j < nimap; ++j)
834 {
835 B_data[i + 2 * k * nbmap +
836 (j + (2 * k + 1) * nimap) * nbndry] -=
837 lambda_imag * (*MassMat)(bmap[i], imap[j]);
838 }
839
840 for (j = 0; j < nimap; ++j)
841 {
842 B_data[i + (2 * k + 1) * nbmap +
843 (j + 2 * k * nimap) * nbndry] +=
844 lambda_imag * (*MassMat)(bmap[i], imap[j]);
845 }
846 }
847 }
848
849 for (k = 0; k < nvel; ++k)
850 {
851 if ((nz_loc == 2) && (k == 2)) // handle d/dz derivative
852 {
853 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
854
855 // Real Component
856 Vmath::Smul(npoints, beta, phys, 1, deriv, 1);
857
858 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
859 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
860 Dbnd->GetRawPtr() +
861 ((nz_loc * k + 1) * bmap.size() + i) *
862 nsize_p[n],
863 1);
864
865 // Imaginary Component
866 Vmath::Neg(psize, pcoeffs, 1);
867 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
868 Dbnd->GetRawPtr() +
869 ((nz_loc * k) * bmap.size() + i) *
870 nsize_p[n] +
871 psize,
872 1);
873
874 // now do advection terms
875 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
876 1, deriv, 1, tmpphys, 1);
877
878 locExp->IProductWRTBase(tmpphys, coeffs);
879
880 // real contribution
881 for (nv = 0; nv < nvel; ++nv)
882 {
883 for (j = 0; j < nbmap; ++j)
884 {
885 Ah_data[j + 2 * nv * nbmap +
886 (i + (2 * nv + 1) * nbmap) * AhRows] +=
887 coeffs[bmap[j]];
888 }
889
890 for (j = 0; j < nimap; ++j)
891 {
892 C_data[i + (2 * nv + 1) * nbmap +
893 (j + 2 * nv * nimap) * nbndry] +=
894 coeffs[imap[j]];
895 }
896 }
897
898 Vmath::Neg(ncoeffs, coeffs, 1);
899 // imaginary contribution
900 for (nv = 0; nv < nvel; ++nv)
901 {
902 for (j = 0; j < nbmap; ++j)
903 {
904 Ah_data[j + (2 * nv + 1) * nbmap +
905 (i + 2 * nv * nbmap) * AhRows] +=
906 coeffs[bmap[j]];
907 }
908
909 for (j = 0; j < nimap; ++j)
910 {
911 C_data[i + 2 * nv * nbmap +
912 (j + (2 * nv + 1) * nimap) * nbndry] +=
913 coeffs[imap[j]];
914 }
915 }
916 }
917 else
918 {
919 if (k < 2)
920 {
921 locExp->PhysDeriv(MultiRegions::DirCartesianMap[k],
922 phys, deriv);
923 Vmath::Vmul(npoints,
924 Advtmp = Advfield[k] + phys_offset, 1,
925 deriv, 1, tmpphys, 1);
926 locExp->IProductWRTBase(tmpphys, coeffs);
927
928 for (nv = 0; nv < nvel * nz_loc; ++nv)
929 {
930 for (j = 0; j < nbmap; ++j)
931 {
932 Ah_data[j + nv * nbmap +
933 (i + nv * nbmap) * AhRows] +=
934 coeffs[bmap[j]];
935 }
936
937 for (j = 0; j < nimap; ++j)
938 {
939 C_data[i + nv * nbmap +
940 (j + nv * nimap) * nbndry] +=
941 coeffs[imap[j]];
942 }
943 }
944
945 // copy into column major storage.
946 m_pressure->GetExp(n)->IProductWRTBase(deriv,
947 pcoeffs);
948 for (j = 0; j < nz_loc; ++j)
949 {
951 psize, &(pcoeffs)[0], 1,
952 Dbnd->GetRawPtr() +
953 ((nz_loc * k + j) * bmap.size() + i) *
954 nsize_p[n] +
955 j * psize,
956 1);
957 }
958 }
959 }
960
961 if (IsLinearNSEquation)
962 {
963 for (nv = 0; nv < nvel; ++nv)
964 {
965 // u' . Grad U terms
966 Vmath::Vmul(npoints, phys, 1,
967 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
968 locExp->IProductWRTBase(tmpphys, coeffs);
969
970 for (size_t n1 = 0; n1 < nz_loc; ++n1)
971 {
972 for (j = 0; j < nbmap; ++j)
973 {
974 Ah_data[j + (k * nz_loc + n1) * nbmap +
975 (i + (nv * nz_loc + n1) * nbmap) *
976 AhRows] += coeffs[bmap[j]];
977 }
978
979 for (j = 0; j < nimap; ++j)
980 {
981 C_data[i + (nv * nz_loc + n1) * nbmap +
982 (j + (k * nz_loc + n1) * nimap) *
983 nbndry] += coeffs[imap[j]];
984 }
985 }
986 }
987 }
988 }
989 }
990
991 for (i = 0; i < nimap; ++i)
992 {
993 // Fill element with mode
994 Vmath::Zero(ncoeffs, coeffs, 1);
995 coeffs[imap[i]] = 1.0;
996 locExp->BwdTrans(coeffs, phys);
997
998 for (k = 0; k < nvel * nz_loc; ++k)
999 {
1000 for (j = 0; j < nbmap; ++j) // C set up as transpose
1001 {
1002 C_data[j + k * nbmap + (i + k * nimap) * nbndry] +=
1003 m_kinvis * HelmMatScale *
1004 HelmMat_data[imap[i] + HelmMatRows * bmap[j]];
1005 }
1006
1007 for (j = 0; j < nimap; ++j)
1008 {
1009 D_data[i + k * nimap + (j + k * nimap) * nint] +=
1010 m_kinvis * HelmMatScale *
1011 HelmMat_data[imap[i] + HelmMatRows * imap[j]];
1012 }
1013 }
1014
1015 if ((lambda_imag != NekConstants::kNekUnsetDouble) &&
1016 (nz_loc == 2))
1017 {
1018 for (k = 0; k < nvel; ++k)
1019 {
1020 for (j = 0; j < nbmap; ++j) // C set up as transpose
1021 {
1022 C_data[j + 2 * k * nbmap +
1023 (i + (2 * k + 1) * nimap) * nbndry] +=
1024 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1025 }
1026
1027 for (j = 0; j < nbmap; ++j) // C set up as transpose
1028 {
1029 C_data[j + (2 * k + 1) * nbmap +
1030 (i + 2 * k * nimap) * nbndry] -=
1031 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1032 }
1033
1034 for (j = 0; j < nimap; ++j)
1035 {
1036 D_data[i + 2 * k * nimap +
1037 (j + (2 * k + 1) * nimap) * nint] -=
1038 lambda_imag * (*MassMat)(imap[i], imap[j]);
1039 }
1040
1041 for (j = 0; j < nimap; ++j)
1042 {
1043 D_data[i + (2 * k + 1) * nimap +
1044 (j + 2 * k * nimap) * nint] +=
1045 lambda_imag * (*MassMat)(imap[i], imap[j]);
1046 }
1047 }
1048 }
1049
1050 for (k = 0; k < nvel; ++k)
1051 {
1052 if ((nz_loc == 2) && (k == 2)) // handle d/dz derivative
1053 {
1054 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
1055
1056 // Real Component
1057 Vmath::Smul(npoints, beta, phys, 1, deriv, 1);
1058 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
1059 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
1060 Dint->GetRawPtr() +
1061 ((nz_loc * k + 1) * imap.size() + i) *
1062 nsize_p[n],
1063 1);
1064 // Imaginary Component
1065 Vmath::Neg(psize, pcoeffs, 1);
1066 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
1067 Dint->GetRawPtr() +
1068 ((nz_loc * k) * imap.size() + i) *
1069 nsize_p[n] +
1070 psize,
1071 1);
1072
1073 // Advfield[k] *d/dx_k to all velocity
1074 // components on diagonal
1075 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
1076 1, deriv, 1, tmpphys, 1);
1077 locExp->IProductWRTBase(tmpphys, coeffs);
1078
1079 // Real Components
1080 for (nv = 0; nv < nvel; ++nv)
1081 {
1082 for (j = 0; j < nbmap; ++j)
1083 {
1084 B_data[j + 2 * nv * nbmap +
1085 (i + (2 * nv + 1) * nimap) * nbndry] +=
1086 coeffs[bmap[j]];
1087 }
1088
1089 for (j = 0; j < nimap; ++j)
1090 {
1091 D_data[j + 2 * nv * nimap +
1092 (i + (2 * nv + 1) * nimap) * nint] +=
1093 coeffs[imap[j]];
1094 }
1095 }
1096 Vmath::Neg(ncoeffs, coeffs, 1);
1097 // Imaginary
1098 for (nv = 0; nv < nvel; ++nv)
1099 {
1100 for (j = 0; j < nbmap; ++j)
1101 {
1102 B_data[j + (2 * nv + 1) * nbmap +
1103 (i + 2 * nv * nimap) * nbndry] +=
1104 coeffs[bmap[j]];
1105 }
1106
1107 for (j = 0; j < nimap; ++j)
1108 {
1109 D_data[j + (2 * nv + 1) * nimap +
1110 (i + 2 * nv * nimap) * nint] +=
1111 coeffs[imap[j]];
1112 }
1113 }
1114 }
1115 else
1116 {
1117 if (k < 2)
1118 {
1119 // Differentiation & Inner product wrt base.
1120 // locExp->PhysDeriv(k,phys, deriv);
1121 locExp->PhysDeriv(MultiRegions::DirCartesianMap[k],
1122 phys, deriv);
1123 Vmath::Vmul(npoints,
1124 Advtmp = Advfield[k] + phys_offset, 1,
1125 deriv, 1, tmpphys, 1);
1126 locExp->IProductWRTBase(tmpphys, coeffs);
1127
1128 for (nv = 0; nv < nvel * nz_loc; ++nv)
1129 {
1130 for (j = 0; j < nbmap; ++j)
1131 {
1132 B_data[j + nv * nbmap +
1133 (i + nv * nimap) * nbndry] +=
1134 coeffs[bmap[j]];
1135 }
1136
1137 for (j = 0; j < nimap; ++j)
1138 {
1139 D_data[j + nv * nimap +
1140 (i + nv * nimap) * nint] +=
1141 coeffs[imap[j]];
1142 }
1143 }
1144 // copy into column major storage.
1145 m_pressure->GetExp(n)->IProductWRTBase(deriv,
1146 pcoeffs);
1147 for (j = 0; j < nz_loc; ++j)
1148 {
1150 psize, &(pcoeffs)[0], 1,
1151 Dint->GetRawPtr() +
1152 ((nz_loc * k + j) * imap.size() + i) *
1153 nsize_p[n] +
1154 j * psize,
1155 1);
1156 }
1157 }
1158 }
1159
1160 if (IsLinearNSEquation)
1161 {
1162 size_t n1;
1163 for (nv = 0; nv < nvel; ++nv)
1164 {
1165 // u'.Grad U terms
1166 Vmath::Vmul(npoints, phys, 1,
1167 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
1168 locExp->IProductWRTBase(tmpphys, coeffs);
1169
1170 for (n1 = 0; n1 < nz_loc; ++n1)
1171 {
1172 for (j = 0; j < nbmap; ++j)
1173 {
1174 B_data[j + (k * nz_loc + n1) * nbmap +
1175 (i + (nv * nz_loc + n1) * nimap) *
1176 nbndry] += coeffs[bmap[j]];
1177 }
1178
1179 for (j = 0; j < nimap; ++j)
1180 {
1181 D_data[j + (k * nz_loc + n1) * nimap +
1182 (i + (nv * nz_loc + n1) * nimap) *
1183 nint] += coeffs[imap[j]];
1184 }
1185 }
1186 }
1187 }
1188 }
1189 }
1190
1191 D->Invert();
1192 (*B) = (*B) * (*D);
1193
1194 // perform (*Ah) = (*Ah) - (*B)*(*C) but since size of
1195 // Ah is larger than (*B)*(*C) easier to call blas
1196 // directly
1197 Blas::Dgemm('N', 'T', B->GetRows(), C->GetRows(), B->GetColumns(),
1198 -1.0, B->GetRawPtr(), B->GetRows(), C->GetRawPtr(),
1199 C->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1200 }
1201
1202 mat.m_BCinv->SetBlock(
1203 n, n,
1205 mat.m_Btilde->SetBlock(
1206 n, n,
1208 mat.m_Cinv->SetBlock(
1209 n, n,
1211 mat.m_D_bnd->SetBlock(
1212 n, n,
1214 mat.m_D_int->SetBlock(
1215 n, n,
1217
1218 // Do matrix manipulations and get final set of block matries
1219 // reset boundary to put mean mode into boundary system.
1220
1221 DNekMatSharedPtr Cinv, BCinv, Btilde;
1222 DNekMat DintCinvDTint, BCinvDTint_m_DTbnd, DintCinvBTtilde_m_Dbnd;
1223
1224 Cinv = D;
1225 BCinv = B;
1226 Btilde = C;
1227
1228 DintCinvDTint = (*Dint) * (*Cinv) * Transpose(*Dint);
1229 BCinvDTint_m_DTbnd = (*BCinv) * Transpose(*Dint) - Transpose(*Dbnd);
1230
1231 // This could be transpose of BCinvDint in some cases
1232 DintCinvBTtilde_m_Dbnd =
1233 (*Dint) * (*Cinv) * Transpose(*Btilde) - (*Dbnd);
1234
1235 // Set up final set of matrices.
1237 nsize_bndry_p1[n], nsize_p_m1[n], zero);
1239 nsize_p_m1[n], nsize_bndry_p1[n], zero);
1241 nsize_p_m1[n], nsize_p_m1[n], zero);
1242 Array<OneD, NekDouble> Dh_data = Dh->GetPtr();
1243
1244 // Copy matrices into final structures.
1245 size_t n1, n2;
1246 for (n1 = 0; n1 < nz_loc; ++n1)
1247 {
1248 for (i = 0; i < psize - 1; ++i)
1249 {
1250 for (n2 = 0; n2 < nz_loc; ++n2)
1251 {
1252 for (j = 0; j < psize - 1; ++j)
1253 {
1254 //(*Dh)(i+n1*(psize-1),j+n2*(psize-1)) =
1255 //-DintCinvDTint(i+1+n1*psize,j+1+n2*psize);
1256 Dh_data[(i + n1 * (psize - 1)) +
1257 (j + n2 * (psize - 1)) * Dh->GetRows()] =
1258 -DintCinvDTint(i + 1 + n1 * psize,
1259 j + 1 + n2 * psize);
1260 }
1261 }
1262 }
1263 }
1264
1265 for (n1 = 0; n1 < nz_loc; ++n1)
1266 {
1267 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1268 {
1269 (*Ah)(i, nsize_bndry_p1[n] - nz_loc + n1) =
1270 BCinvDTint_m_DTbnd(i, n1 * psize);
1271 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1, i) =
1272 DintCinvBTtilde_m_Dbnd(n1 * psize, i);
1273 }
1274 }
1275
1276 for (n1 = 0; n1 < nz_loc; ++n1)
1277 {
1278 for (n2 = 0; n2 < nz_loc; ++n2)
1279 {
1280 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1,
1281 nsize_bndry_p1[n] - nz_loc + n2) =
1282 -DintCinvDTint(n1 * psize, n2 * psize);
1283 }
1284 }
1285
1286 for (n1 = 0; n1 < nz_loc; ++n1)
1287 {
1288 for (j = 0; j < psize - 1; ++j)
1289 {
1290 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1291 {
1292 (*Bh)(i, j + n1 * (psize - 1)) =
1293 BCinvDTint_m_DTbnd(i, j + 1 + n1 * psize);
1294 (*Ch)(j + n1 * (psize - 1), i) =
1295 DintCinvBTtilde_m_Dbnd(j + 1 + n1 * psize, i);
1296 }
1297 }
1298 }
1299
1300 for (n1 = 0; n1 < nz_loc; ++n1)
1301 {
1302 for (n2 = 0; n2 < nz_loc; ++n2)
1303 {
1304 for (j = 0; j < psize - 1; ++j)
1305 {
1306 (*Bh)(nsize_bndry_p1[n] - nz_loc + n1,
1307 j + n2 * (psize - 1)) =
1308 -DintCinvDTint(n1 * psize, j + 1 + n2 * psize);
1309 (*Ch)(j + n2 * (psize - 1),
1310 nsize_bndry_p1[n] - nz_loc + n1) =
1311 -DintCinvDTint(j + 1 + n2 * psize, n1 * psize);
1312 }
1313 }
1314 }
1315
1316 // Do static condensation
1317 Dh->Invert();
1318 (*Bh) = (*Bh) * (*Dh);
1319 //(*Ah) = (*Ah) - (*Bh)*(*Ch);
1320 Blas::Dgemm('N', 'N', Bh->GetRows(), Ch->GetColumns(), Bh->GetColumns(),
1321 -1.0, Bh->GetRawPtr(), Bh->GetRows(), Ch->GetRawPtr(),
1322 Ch->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1323
1324 // Set matrices for later inversion. Probably do not need to be
1325 // attached to class
1326 pAh->SetBlock(
1327 n, n,
1329 pBh->SetBlock(
1330 n, n,
1332 pCh->SetBlock(
1333 n, n,
1335 pDh->SetBlock(
1336 n, n,
1338 }
1339 timer.Stop();
1340 std::cout << "Matrix Setup Costs: " << timer.TimePerTest(1) << std::endl;
1341
1342 timer.Start();
1343 // Set up global coupled boundary solver.
1344 // This is a key to define the solution matrix type
1345 // currently we are giving it a argument of eLInearAdvectionReaction
1346 // since this then makes the matrix storage of type eFull
1348 locToGloMap);
1349 mat.m_CoupledBndSys =
1351 AllocateSharedPtr(key, m_fields[0], pAh, pBh, pCh, pDh,
1352 locToGloMap);
1353 mat.m_CoupledBndSys->Initialise(locToGloMap);
1354}
1355
1357{
1358 SolverUtils::AddSummaryItem(s, "Solver Type", "Coupled Linearised NS");
1359}
1360
1361void CoupledLinearNS::v_DoInitialise(bool dumpInitialConditions)
1362{
1363 switch (m_equationType)
1364 {
1365 case eUnsteadyStokes:
1367 {
1368 // LibUtilities::TimeIntegrationMethod intMethod;
1369 // std::string TimeIntStr =
1370 // m_session->GetSolverInfo("TIMEINTEGRATIONMETHOD");
1371 // int i;
1372 // for(i = 0; i < (int)
1373 // LibUtilities::SIZE_TimeIntegrationMethod; ++i)
1374 // {
1375 // if(boost::iequals(LibUtilities::TimeIntegrationMethodMap[i],TimeIntStr))
1376 // {
1377 // intMethod =
1378 // (LibUtilities::TimeIntegrationMethod)i;
1379 // break;
1380 // }
1381 // }
1382 //
1383 // ASSERTL0(i != (int)
1384 // LibUtilities::SIZE_TimeIntegrationMethod, "Invalid
1385 // time integration type.");
1386 //
1387 // m_integrationScheme =
1388 // LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(LibUtilities::TimeIntegrationMethodMap[intMethod]);
1389
1390 // Could defind this from IncNavierStokes class?
1392
1395
1396 // Set initial condition using time t=0
1397
1398 SetInitialConditions(0.0, dumpInitialConditions);
1399 break;
1400 }
1401 case eSteadyStokes:
1402 SetUpCoupledMatrix(0.0);
1403 break;
1404 case eSteadyOseen:
1405 {
1407 for (size_t i = 0; i < m_velocity.size(); ++i)
1408 {
1409 AdvField[i] = Array<OneD, NekDouble>(
1410 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1411 }
1412
1413 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1414 "Advection Velocity section must be defined in "
1415 "session file.");
1416
1417 std::vector<std::string> fieldStr;
1418 for (size_t i = 0; i < m_velocity.size(); ++i)
1419 {
1420 fieldStr.push_back(
1421 m_boundaryConditions->GetVariable(m_velocity[i]));
1422 }
1423 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1424
1425 SetUpCoupledMatrix(0.0, AdvField, false);
1426 }
1427 break;
1429 {
1430 m_session->LoadParameter("KinvisMin", m_kinvisMin);
1431 m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1432 m_session->LoadParameter("Tolerence", m_tol);
1433 m_session->LoadParameter("MaxIteration", m_maxIt);
1434 m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1435 m_session->LoadParameter("Restart", m_Restart);
1436
1438
1439 if (m_Restart == 1)
1440 {
1441 ASSERTL0(m_session->DefinesFunction("Restart"),
1442 "Restart section must be defined in session file.");
1443
1445 for (size_t i = 0; i < m_velocity.size(); ++i)
1446 {
1447 Restart[i] = Array<OneD, NekDouble>(
1448 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1449 }
1450 std::vector<std::string> fieldStr;
1451 for (size_t i = 0; i < m_velocity.size(); ++i)
1452 {
1453 fieldStr.push_back(
1454 m_boundaryConditions->GetVariable(m_velocity[i]));
1455 }
1456 GetFunction("Restart")->Evaluate(fieldStr, Restart);
1457
1458 for (size_t i = 0; i < m_velocity.size(); ++i)
1459 {
1460 m_fields[m_velocity[i]]->FwdTransLocalElmt(
1461 Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1462 }
1463 std::cout << "Saving the RESTART file for m_kinvis = "
1464 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1465 << std::endl;
1466 }
1467 else // We solve the Stokes Problem
1468 {
1469
1470 SetUpCoupledMatrix(0.0);
1471 m_initialStep = true;
1472 m_counter = 1;
1473 // SolveLinearNS(m_ForcingTerm_Coeffs);
1474 Solve();
1475 m_initialStep = false;
1476 std::cout << "Saving the Stokes Flow for m_kinvis = "
1477 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1478 << std::endl;
1479 }
1480 }
1481 break;
1483 {
1484 SetInitialConditions(0.0, dumpInitialConditions);
1485
1487 for (size_t i = 0; i < m_velocity.size(); ++i)
1488 {
1489 AdvField[i] = Array<OneD, NekDouble>(
1490 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1491 }
1492
1493 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1494 "Advection Velocity section must be defined in "
1495 "session file.");
1496
1497 std::vector<std::string> fieldStr;
1498 for (size_t i = 0; i < m_velocity.size(); ++i)
1499 {
1500 fieldStr.push_back(
1501 m_boundaryConditions->GetVariable(m_velocity[i]));
1502 }
1503 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1504
1505 SetUpCoupledMatrix(m_lambda, AdvField, true);
1506 }
1507 break;
1508 case eNoEquationType:
1509 default:
1510 ASSERTL0(false,
1511 "Unknown or undefined equation type for CoupledLinearNS");
1512 }
1513}
1514
1516 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1517 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
1518{
1519 // evaluate convection terms
1520 EvaluateAdvectionTerms(inarray, outarray, time);
1521
1522 for (auto &x : m_forcing)
1523 {
1524 x->Apply(m_fields, outarray, outarray, time);
1525 }
1526}
1527
1529 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1530 Array<OneD, Array<OneD, NekDouble>> &outarray,
1531 [[maybe_unused]] const NekDouble time, const NekDouble aii_Dt)
1532{
1533 size_t i;
1535 NekDouble lambda = 1.0 / aii_Dt;
1536 static NekDouble lambda_store;
1538 // Matrix solution
1539 if (fabs(lambda_store - lambda) > 1e-10)
1540 {
1541 SetUpCoupledMatrix(lambda);
1542 lambda_store = lambda;
1543 }
1544
1545 // Forcing for advection solve
1546 for (i = 0; i < m_velocity.size(); ++i)
1547 {
1548 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1549 m_fields[m_velocity[i]]->SetWaveSpace(true);
1550 m_fields[m_velocity[i]]->IProductWRTBase(
1551 inarray[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1552 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1554 m_fields[m_velocity[i]]->GetCoeffs(), 1,
1555 m_fields[m_velocity[i]]->UpdateCoeffs(), 1);
1556 forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1557 }
1558
1559 SolveLinearNS(forcing);
1560
1561 for (i = 0; i < m_velocity.size(); ++i)
1562 {
1563 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1564 outarray[i]);
1565 }
1566}
1567
1569{
1570 size_t nfields = m_fields.size();
1571 for (size_t k = 0; k < nfields; ++k)
1572 {
1573 // Backward Transformation in physical space for time evolution
1574 m_fields[k]->BwdTrans(m_fields[k]->GetCoeffs(),
1575 m_fields[k]->UpdatePhys());
1576 }
1577}
1578
1580{
1581 size_t nfields = m_fields.size();
1582 for (size_t k = 0; k < nfields; ++k)
1583 {
1584 // Forward Transformation in physical space for time evolution
1585 m_fields[k]->FwdTransLocalElmt(m_fields[k]->GetPhys(),
1586 m_fields[k]->UpdateCoeffs());
1587 }
1588}
1589
1591{
1592 switch (m_equationType)
1593 {
1594 case eUnsteadyStokes:
1596 // AdvanceInTime(m_steps);
1598 break;
1599 case eSteadyStokes:
1600 case eSteadyOseen:
1602 {
1603 Solve();
1604 break;
1605 }
1607 {
1608 LibUtilities::Timer Generaltimer;
1609 Generaltimer.Start();
1610
1611 int Check(0);
1612
1613 // Saving the init datas
1614 Checkpoint_Output(Check);
1615 Check++;
1616
1617 std::cout
1618 << "We execute INITIALLY SolveSteadyNavierStokes for m_kinvis "
1619 "= "
1620 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1621 << std::endl;
1623
1624 while (m_kinvis > m_kinvisMin)
1625 {
1626 if (Check == 1)
1627 {
1628 std::cout
1629 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1630 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1631 << std::endl;
1633 Checkpoint_Output(Check);
1634 Check++;
1635 }
1636
1637 Continuation();
1638
1639 if (m_kinvis > m_kinvisMin)
1640 {
1641 std::cout
1642 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1643 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1644 << std::endl;
1646 Checkpoint_Output(Check);
1647 Check++;
1648 }
1649 }
1650
1651 Generaltimer.Stop();
1652 std::cout << "\nThe total calculation time is : "
1653 << Generaltimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1654
1655 break;
1656 }
1657 case eNoEquationType:
1658 default:
1659 ASSERTL0(false,
1660 "Unknown or undefined equation type for CoupledLinearNS");
1661 }
1662}
1663
1664/** Virtual function to define if operator in DoSolve is negated
1665 * with regard to the strong form. This is currently only used in
1666 * Arnoldi solves. For Coupledd solver this is true since Stokes
1667 * operator is set up as a LHS rather than RHS operation
1668 */
1670{
1671 return true;
1672}
1673
1675{
1676 const size_t ncmpt = m_velocity.size();
1677 Array<OneD, Array<OneD, NekDouble>> forcing_phys(ncmpt);
1678 Array<OneD, Array<OneD, NekDouble>> forcing(ncmpt);
1679
1680 for (size_t i = 0; i < ncmpt; ++i)
1681 {
1682 forcing_phys[i] =
1684 forcing[i] =
1686 }
1687
1688 for (auto &x : m_forcing)
1689 {
1690 const NekDouble time = 0;
1691 x->Apply(m_fields, forcing_phys, forcing_phys, time);
1692 }
1693 for (size_t i = 0; i < ncmpt; ++i)
1694 {
1695 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1696 m_fields[i]->SetWaveSpace(true);
1697 m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1698 m_fields[i]->SetWaveSpace(waveSpace);
1699 }
1700
1701 SolveLinearNS(forcing);
1702}
1703
1705{
1709
1710 for (size_t i = 0; i < m_velocity.size(); ++i)
1711 {
1713 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1716 }
1717
1718 if (m_session->DefinesFunction("ForcingTerm"))
1719 {
1720 std::vector<std::string> fieldStr;
1721 for (size_t i = 0; i < m_velocity.size(); ++i)
1722 {
1723 fieldStr.push_back(
1724 m_boundaryConditions->GetVariable(m_velocity[i]));
1725 }
1726 GetFunction("ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1727 for (size_t i = 0; i < m_velocity.size(); ++i)
1728 {
1729 m_fields[m_velocity[i]]->FwdTransLocalElmt(m_ForcingTerm[i],
1731 }
1732 }
1733 else
1734 {
1735 std::cout
1736 << "'ForcingTerm' section has not been defined in the input file "
1737 "=> forcing=0"
1738 << std::endl;
1739 }
1740}
1741
1743{
1744 LibUtilities::Timer Newtontimer;
1745 Newtontimer.Start();
1746
1749 Array<OneD, Array<OneD, NekDouble>> delta_velocity_Phys(m_velocity.size());
1750 Array<OneD, Array<OneD, NekDouble>> Velocity_Phys(m_velocity.size());
1751 Array<OneD, NekDouble> L2_norm(m_velocity.size(), 1.0);
1752 Array<OneD, NekDouble> Inf_norm(m_velocity.size(), 1.0);
1753
1754 for (size_t i = 0; i < m_velocity.size(); ++i)
1755 {
1756 delta_velocity_Phys[i] = Array<OneD, NekDouble>(
1757 m_fields[m_velocity[i]]->GetTotPoints(), 1.0);
1758 Velocity_Phys[i] = Array<OneD, NekDouble>(
1759 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1760 }
1761
1762 m_counter = 1;
1763
1764 L2Norm(delta_velocity_Phys, L2_norm);
1765
1766 // while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1767 while (std::max(L2_norm[0], L2_norm[1]) > m_tol)
1768 {
1769 if (m_counter == 1)
1770 // At the first Newton step, we use the solution of the
1771 // Stokes problem (if Restart=0 in input file) Or the
1772 // solution of the .rst file (if Restart=1 in input
1773 // file)
1774 {
1775 for (size_t i = 0; i < m_velocity.size(); ++i)
1776 {
1777 RHS_Coeffs[i] = Array<OneD, NekDouble>(
1778 m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1779 RHS_Phys[i] = Array<OneD, NekDouble>(
1780 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1781 }
1782
1783 for (size_t i = 0; i < m_velocity.size(); ++i)
1784 {
1785 m_fields[m_velocity[i]]->BwdTrans(
1786 m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1787 }
1788
1789 m_initialStep = true;
1790 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1791 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1792 SolveLinearNS(RHS_Coeffs);
1793 m_initialStep = false;
1794 }
1795 if (m_counter > 1)
1796 {
1797 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1798
1800 0) // Setting Up the matrix is expensive. We do it at each
1801 // "m_MatrixSetUpStep" step.
1802 {
1803 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1804 }
1805 SolveLinearNS(RHS_Coeffs);
1806 }
1807
1808 for (size_t i = 0; i < m_velocity.size(); ++i)
1809 {
1810 m_fields[m_velocity[i]]->BwdTrans(RHS_Coeffs[i], RHS_Phys[i]);
1811 m_fields[m_velocity[i]]->BwdTrans(
1812 m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1813 }
1814
1815 for (size_t i = 0; i < m_velocity.size(); ++i)
1816 {
1817 Vmath::Vadd(Velocity_Phys[i].size(), Velocity_Phys[i], 1,
1818 delta_velocity_Phys[i], 1, Velocity_Phys[i], 1);
1819 }
1820
1821 // InfNorm(delta_velocity_Phys, Inf_norm);
1822 L2Norm(delta_velocity_Phys, L2_norm);
1823
1824 if (std::max(Inf_norm[0], Inf_norm[1]) > 100)
1825 {
1826 std::cout << "\nThe Newton method has failed at m_kinvis = "
1827 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1828 << std::endl;
1829 ASSERTL0(0, "The Newton method has failed... \n");
1830 }
1831
1832 std::cout << "\n";
1833 m_counter++;
1834 }
1835
1836 if (m_counter > 1) // We save u:=u+\delta u in u->Coeffs
1837 {
1838 for (size_t i = 0; i < m_velocity.size(); ++i)
1839 {
1840 m_fields[m_velocity[i]]->FwdTrans(
1841 Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1842 }
1843 }
1844
1845 Newtontimer.Stop();
1846 std::cout << "We have done " << m_counter - 1 << " iteration(s) in "
1847 << Newtontimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1848}
1849
1851{
1856
1857 std::cout << "We apply the continuation method: " << std::endl;
1858
1859 for (size_t i = 0; i < m_velocity.size(); ++i)
1860 {
1862 0.0);
1863 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1864 u_N[i]);
1865
1867 0.0);
1868 tmp_RHS[i] = Array<OneD, NekDouble>(
1869 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1870
1871 m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1872 Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1873
1874 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1875 m_fields[m_velocity[i]]->SetWaveSpace(true);
1876 m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1877 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1878 }
1879
1880 SetUpCoupledMatrix(0.0, u_N, true);
1881 SolveLinearNS(RHS);
1882
1883 for (size_t i = 0; i < m_velocity.size(); ++i)
1884 {
1885 u_star[i] = Array<OneD, NekDouble>(
1886 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1887 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1888 u_star[i]);
1889
1890 // u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1891 Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1892 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1893
1894 m_fields[m_velocity[i]]->FwdTrans(
1895 u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1896 }
1897
1899}
1900
1902 Array<OneD, NekDouble> &outarray)
1903{
1904 for (size_t i = 0; i < m_velocity.size(); ++i)
1905 {
1906 outarray[i] = 0.0;
1907 for (size_t j = 0; j < inarray[i].size(); ++j)
1908 {
1909 if (inarray[i][j] > outarray[i])
1910 {
1911 outarray[i] = inarray[i][j];
1912 }
1913 }
1914 std::cout << "InfNorm[" << i << "] = " << outarray[i] << std::endl;
1915 }
1916}
1917
1919 Array<OneD, NekDouble> &outarray)
1920{
1921 for (size_t i = 0; i < m_velocity.size(); ++i)
1922 {
1923 outarray[i] = 0.0;
1924 for (size_t j = 0; j < inarray[i].size(); ++j)
1925 {
1926 outarray[i] += inarray[i][j] * inarray[i][j];
1927 }
1928 outarray[i] = std::sqrt(outarray[i]);
1929 std::cout << "L2Norm[" << i << "] = " << outarray[i] << std::endl;
1930 }
1931}
1932
1934 Array<OneD, Array<OneD, NekDouble>> &Velocity,
1935 Array<OneD, Array<OneD, NekDouble>> &outarray)
1936{
1942
1943 for (size_t i = 0; i < m_velocity.size(); ++i)
1944 {
1945 Eval_Adv[i] = Array<OneD, NekDouble>(
1946 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1947 tmp_DerVel[i] = Array<OneD, NekDouble>(
1948 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1949
1950 AdvTerm[i] = Array<OneD, NekDouble>(
1951 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1952 ViscTerm[i] = Array<OneD, NekDouble>(
1953 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1954 Forc[i] = Array<OneD, NekDouble>(
1955 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1956 outarray[i] = Array<OneD, NekDouble>(
1957 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1958
1959 m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1960
1961 Vmath::Smul(tmp_DerVel[i].size(), m_kinvis, tmp_DerVel[i], 1,
1962 tmp_DerVel[i], 1);
1963 }
1964
1965 EvaluateAdvectionTerms(Velocity, Eval_Adv, m_time);
1966
1967 for (size_t i = 0; i < m_velocity.size(); ++i)
1968 {
1969 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1970 m_fields[m_velocity[i]]->SetWaveSpace(true);
1971 m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i],
1972 AdvTerm[i]); //(w, (u.grad)u)
1973 m_fields[m_velocity[i]]->IProductWRTDerivBase(
1974 i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1975 m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i],
1976 Forc[i]); //(w, f)
1977 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1978
1979 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1980 outarray[i], 1);
1981 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1982 outarray[i], 1);
1983
1984 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1985 1);
1986 }
1987}
1988
1991{
1993 returnval =
1995
1996 int nummodes;
1997
1998 for (auto &expMapIter : VelExp)
1999 {
2001
2002 for (size_t i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
2003 {
2004 LibUtilities::BasisKey B = expMapIter.second->m_basisKeyVector[i];
2005 nummodes = B.GetNumModes();
2006 ASSERTL0(nummodes > 3,
2007 "Velocity polynomial space not sufficiently high (>= 4)");
2008 // Should probably set to be an orthogonal basis.
2009 LibUtilities::BasisKey newB(B.GetBasisType(), nummodes - 2,
2010 B.GetPointsKey());
2011 BasisVec.push_back(newB);
2012 }
2013
2014 // Put new expansion into list.
2015 SpatialDomains::ExpansionInfoShPtr expansionElementShPtr =
2017 expMapIter.second->m_geomShPtr, BasisVec);
2018 (*returnval)[expMapIter.first] = expansionElementShPtr;
2019 }
2020
2021 // Save expansion into graph.
2022 m_graph->SetExpansionInfo("p", returnval);
2023
2024 return *returnval;
2025}
2026
2027/**
2028 * @param forcing A list of forcing functions for each velocity
2029 * component
2030 *
2031 * The routine involves two levels of static
2032 * condensations. Initially we require a statically condensed
2033 * forcing function which requires the following manipulation
2034 *
2035 * \f[ {F\_bnd} = {\bf f}_{bnd} -m\_B \,m\_Cinv\, {\bf f}_{int},
2036 * \hspace{1cm} F\_p = m\_D\_{int}\, m\_Cinv\, {\bf f}_{int} \f]
2037 *
2038 * Where \f${\bf f}_{bnd}\f$ denote the forcing degrees of
2039 * freedom of the elemental velocities on the boundary of the
2040 * element, \f${\bf f}_{int}\f$ denote the forcing degrees of
2041 * freedom of the elemental velocities on the interior of the
2042 * element. (see detailed description for more details).
2043 *
2044 * This vector is further manipulated into
2045 *
2046 * \f[ Fh\_{bnd} = \left [ \begin{array}{c} f\_{bnd} -m\_B \,
2047 * m\_Cinv\, {\bf f}_{int}\\ \left [m\_D\_{int} \, m\_Cinv \,{\bf
2048 * f}_{int} \right]_0 \end{array}\right ] \hspace{1cm} [Fh\_p]_{i} =
2049 * \begin{array}{c} [m\_D\_{int} \, m\_Cinv \, {\bf
2050 * f}_{int}]_{i+1} \end{array} \f]
2051 *
2052 * where \f$-{[}m\_D\_{int}^T\, m\_Cinv \,{\bf f}_{int}]_0\f$
2053 * which is corresponds to the mean mode of the pressure degrees
2054 * of freedom and is now added to the boundary system and the
2055 * remainder of the block becomes the interior forcing for the
2056 * inner static condensation (see detailed description for more
2057 * details) which is set up in a GlobalLinSysDirectStaticCond
2058 * class.
2059 *
2060 * Finally we perform the final maniplation of the forcing to
2061 * using hte
2062 * \f[ Fh\_{bnd} = Fh\_{bnd} - m\_Bh \,m\_Chinv \, Fh\_p \f]
2063 *
2064 * We can now call the solver to the global coupled boundary
2065 * system (through the call to #m_CoupledBndSys->Solve) to obtain
2066 * the velocity boundary solution as the mean pressure solution,
2067 * i.e.
2068 *
2069 * \f[ {\cal A}^T(\hat{A} - \hat{C}^T \hat{D}^{-1} \hat{B} ){\cal
2070 * A} \, Bnd = Fh\_{bnd} \f]
2071 *
2072 * Once we know the solution to the above the rest of the pressure
2073 * modes are recoverable thorugh
2074 *
2075 * \f[ Ph = m\_Dhinv\, (Bnd - m\_Ch^T \, Fh_{bnd}) \f]
2076 *
2077 * We can now unpack \f$ Fh\_{bnd} \f$ (last elemental mode) and
2078 * \f$ Ph \f$ into #m_pressure and \f$ F_p\f$ and \f$ Fh\_{bnd}\f$
2079 * into a closed pack list of boundary velocoity degrees of
2080 * freedom stored in \f$ F\_bnd\f$.
2081 *
2082 * Finally the interior velocity degrees of freedom are then
2083 * obtained through the relationship
2084 *
2085 * \f[ F\_{int} = m\_Cinv\ ( F\_{int} + m\_D\_{int}^T\,
2086 * F\_p - m\_Btilde^T\, Bnd) \f]
2087 *
2088 * We then unpack the solution back to the MultiRegion structures
2089 * #m_velocity and #m_pressure
2090 */
2092 const Array<OneD, Array<OneD, NekDouble>> &forcing)
2093{
2094 size_t i;
2097
2098 // Impose Dirichlet conditions on velocity fields
2099 for (i = 0; i < m_velocity.size(); ++i)
2100 {
2101 Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
2102 m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
2103 }
2104
2106 {
2107 int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
2108 for (int n = 0; n < m_npointsZ / 2; ++n)
2109 {
2110 // Get the a Fourier mode of velocity and forcing components.
2111 for (i = 0; i < m_velocity.size(); ++i)
2112 {
2113 vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2 * n);
2114 // Note this needs to correlate with how we pass forcing
2115 force[i] = forcing[i] + 2 * n * ncoeffsplane;
2116 }
2117
2118 SolveLinearNS(force, vel_fields, m_pressure->GetPlane(2 * n), n);
2119 }
2120 for (i = 0; i < m_velocity.size(); ++i)
2121 {
2122 m_fields[m_velocity[i]]->SetPhysState(false);
2123 }
2124 m_pressure->SetPhysState(false);
2125 }
2126 else
2127 {
2128 for (i = 0; i < m_velocity.size(); ++i)
2129 {
2130 vel_fields[i] = m_fields[m_velocity[i]];
2131 // Note this needs to correlate with how we pass forcing
2132 force[i] = forcing[i];
2133 }
2134 SolveLinearNS(force, vel_fields, m_pressure);
2135 }
2136}
2137
2142{
2143 size_t i, j, k, n, cnt, cnt1;
2144 size_t nbnd, nint, offset;
2145 size_t nvel = m_velocity.size();
2146 size_t nel = fields[0]->GetNumElmts();
2147 Array<OneD, unsigned int> bmap, imap;
2148
2149 Array<OneD, NekDouble> f_bnd(m_mat[mode].m_BCinv->GetRows());
2150 NekVector<NekDouble> F_bnd(f_bnd.size(), f_bnd, eWrapper);
2151 Array<OneD, NekDouble> f_int(m_mat[mode].m_BCinv->GetColumns());
2152 NekVector<NekDouble> F_int(f_int.size(), f_int, eWrapper);
2153
2154 size_t nz_loc;
2155 int nplanecoeffs =
2156 fields[m_velocity[0]]
2157 ->GetNcoeffs(); // this is fine since we pass the nplane coeff data.
2158
2159 if (mode) // Homogeneous mode flag
2160 {
2161 nz_loc = 2;
2162 }
2163 else
2164 {
2165 if (m_singleMode)
2166 {
2167 nz_loc = 2;
2168 }
2169 else
2170 {
2171 nz_loc = 1;
2173 {
2175 // Zero fields to set complex mode to zero;
2176 for (i = 0; i < fields.size(); ++i)
2177 {
2178 Vmath::Zero(nplanecoeffs,
2179 tmp = fields[i]->UpdateCoeffs() + nplanecoeffs,
2180 1);
2181 }
2182 Vmath::Zero(2 * pressure->GetNcoeffs(),
2183 pressure->UpdateCoeffs(), 1);
2184 }
2185 }
2186 }
2187
2188 for (k = 0; k < nvel; ++k)
2189 {
2191 std::dynamic_pointer_cast<MultiRegions::ContField>(fields[k]);
2192
2194 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsSign();
2195 const Array<OneD, const int> map =
2196 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsMap();
2197
2198 // Add weak boundary conditions to forcing
2200 fields[k]->GetBndConditions();
2202
2204 {
2205 bndCondExp =
2206 m_fields[k]->GetPlane(2 * mode)->GetBndCondExpansions();
2207 }
2208 else
2209 {
2210 bndCondExp = m_fields[k]->GetBndCondExpansions();
2211 }
2212
2213 for (n = 0; n < nz_loc; ++n)
2214 {
2215 int bndcnt = 0;
2216 for (i = 0; i < bndCondExp.size(); ++i)
2217 {
2218 const Array<OneD, const NekDouble> bndcoeffs =
2219 bndCondExp[i]->GetCoeffs();
2220 size_t nCoeffs = (bndCondExp[i])->GetNcoeffs();
2221
2222 cnt = 0;
2223 if (bndConds[i]->GetBoundaryConditionType() ==
2225 bndConds[i]->GetBoundaryConditionType() ==
2227 {
2228 if (m_locToGloMap[mode]->GetSignChange())
2229 {
2230 for (j = 0; j < nCoeffs; j++)
2231 {
2232 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2233 sign[bndcnt + j] * bndcoeffs[j];
2234 }
2235 }
2236 else
2237 {
2238 for (j = 0; j < nCoeffs; j++)
2239 {
2240 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2241 bndcoeffs[j];
2242 }
2243 }
2244 }
2245
2246 bndcnt += bndCondExp[i]->GetNcoeffs();
2247 }
2248 }
2249 }
2250
2251 Array<OneD, NekDouble> bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(), 0.0);
2252
2253 // Construct f_bnd and f_int and fill in bnd from inarray
2254 // (with Dirichlet BCs imposed)
2255 int bndoffset = 0;
2256 cnt = cnt1 = 0;
2257 for (i = 0; i < nel; ++i) // loop over elements
2258 {
2259 fields[m_velocity[0]]->GetExp(i)->GetBoundaryMap(bmap);
2260 fields[m_velocity[0]]->GetExp(i)->GetInteriorMap(imap);
2261 nbnd = bmap.size();
2262 nint = imap.size();
2263 offset = fields[m_velocity[0]]->GetCoeff_Offset(i);
2264
2265 for (j = 0; j < nvel; ++j) // loop over velocity fields
2266 {
2267 Array<OneD, NekDouble> incoeffs = fields[j]->UpdateCoeffs();
2268
2269 for (n = 0; n < nz_loc; ++n)
2270 {
2271 for (k = 0; k < nbnd; ++k)
2272 {
2273 f_bnd[cnt + k] =
2274 forcing[j][n * nplanecoeffs + offset + bmap[k]];
2275
2276 bnd[bndoffset + (n + j * nz_loc) * nbnd + k] =
2277 incoeffs[n * nplanecoeffs + offset + bmap[k]];
2278 }
2279 for (k = 0; k < nint; ++k)
2280 {
2281 f_int[cnt1 + k] =
2282 forcing[j][n * nplanecoeffs + offset + imap[k]];
2283 }
2284
2285 cnt += nbnd;
2286 cnt1 += nint;
2287 }
2288 }
2289 bndoffset +=
2290 nvel * nz_loc * nbnd + nz_loc * (pressure->GetExp(i)->GetNcoeffs());
2291 }
2292
2293 Array<OneD, NekDouble> f_p(m_mat[mode].m_D_int->GetRows());
2294 NekVector<NekDouble> F_p(f_p.size(), f_p, eWrapper);
2295 NekVector<NekDouble> F_p_tmp(m_mat[mode].m_Cinv->GetRows());
2296
2297 // fbnd does not currently hold the pressure mean
2298 F_bnd = F_bnd - (*m_mat[mode].m_BCinv) * F_int;
2299 F_p_tmp = (*m_mat[mode].m_Cinv) * F_int;
2300 F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
2301
2302 // construct inner forcing
2303 Array<OneD, NekDouble> fh_bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(),
2304 0.0);
2305
2306 offset = cnt = 0;
2307 for (i = 0; i < nel; ++i)
2308 {
2309 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2310
2311 for (j = 0; j < nvel; ++j)
2312 {
2313 for (k = 0; k < nbnd; ++k)
2314 {
2315 fh_bnd[offset + j * nbnd + k] = f_bnd[cnt + k];
2316 }
2317 cnt += nbnd;
2318 }
2319
2320 nint = pressure->GetExp(i)->GetNcoeffs();
2321 offset += nvel * nbnd + nint * nz_loc;
2322 }
2323
2324 offset = cnt1 = 0;
2325 for (i = 0; i < nel; ++i)
2326 {
2327 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2328 nint = pressure->GetExp(i)->GetNcoeffs();
2329
2330 for (n = 0; n < nz_loc; ++n)
2331 {
2332 for (j = 0; j < nint; ++j)
2333 {
2334 fh_bnd[offset + nvel * nbnd + n * nint + j] = f_p[cnt1 + j];
2335 }
2336 cnt1 += nint;
2337 }
2338 offset += nvel * nbnd + nz_loc * nint;
2339 }
2340 m_mat[mode].m_CoupledBndSys->Solve(fh_bnd, bnd, m_locToGloMap[mode]);
2341
2342 // unpack pressure and velocity boundary systems.
2343 offset = cnt = 0;
2344 int totpcoeffs = pressure->GetNcoeffs();
2345 Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2346 for (i = 0; i < nel; ++i)
2347 {
2348 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2349 nint = pressure->GetExp(i)->GetNcoeffs();
2350 for (j = 0; j < nvel; ++j)
2351 {
2352 for (k = 0; k < nbnd; ++k)
2353 {
2354 f_bnd[cnt + k] = bnd[offset + j * nbnd + k];
2355 }
2356 cnt += nbnd;
2357 }
2358 offset += nvel * nbnd + nint * nz_loc;
2359 }
2360
2361 pressure->SetPhysState(false);
2362
2363 offset = cnt = cnt1 = 0;
2364 for (i = 0; i < nel; ++i)
2365 {
2366 nint = pressure->GetExp(i)->GetNcoeffs();
2367 nbnd = fields[0]->GetExp(i)->NumBndryCoeffs();
2368 cnt1 = pressure->GetCoeff_Offset(i);
2369
2370 for (n = 0; n < nz_loc; ++n)
2371 {
2372 for (j = 0; j < nint; ++j)
2373 {
2374 p_coeffs[n * totpcoeffs + cnt1 + j] = f_p[cnt + j] =
2375 bnd[offset + nvel * nz_loc * nbnd + n * nint + j];
2376 }
2377 cnt += nint;
2378 }
2379 offset += (nvel * nbnd + nint) * nz_loc;
2380 }
2381
2382 // Back solve first level of static condensation for interior
2383 // velocity space and store in F_int
2384 F_int = F_int + Transpose(*m_mat[mode].m_D_int) * F_p -
2385 Transpose(*m_mat[mode].m_Btilde) * F_bnd;
2386 F_int = (*m_mat[mode].m_Cinv) * F_int;
2387
2388 // Unpack solution from Bnd and F_int to v_coeffs
2389 cnt = cnt1 = 0;
2390 for (i = 0; i < nel; ++i) // loop over elements
2391 {
2392 fields[0]->GetExp(i)->GetBoundaryMap(bmap);
2393 fields[0]->GetExp(i)->GetInteriorMap(imap);
2394 nbnd = bmap.size();
2395 nint = imap.size();
2396 offset = fields[0]->GetCoeff_Offset(i);
2397
2398 for (j = 0; j < nvel; ++j) // loop over velocity fields
2399 {
2400 for (n = 0; n < nz_loc; ++n)
2401 {
2402 for (k = 0; k < nbnd; ++k)
2403 {
2404 fields[j]->SetCoeff(n * nplanecoeffs + offset + bmap[k],
2405 f_bnd[cnt + k]);
2406 }
2407
2408 for (k = 0; k < nint; ++k)
2409 {
2410 fields[j]->SetCoeff(n * nplanecoeffs + offset + imap[k],
2411 f_int[cnt1 + k]);
2412 }
2413 cnt += nbnd;
2414 cnt1 += nint;
2415 }
2416 }
2417 }
2418
2419 for (j = 0; j < nvel; ++j)
2420 {
2421 fields[j]->SetPhysState(false);
2422 }
2423}
2424
2426{
2427 std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size() + 1);
2428 std::vector<std::string> variables(m_fields.size() + 1);
2429 size_t i;
2430
2431 for (i = 0; i < m_fields.size(); ++i)
2432 {
2433 fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2434 variables[i] = m_boundaryConditions->GetVariable(i);
2435 }
2436
2437 fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2438 // project pressure field to velocity space
2439 if (m_singleMode == true)
2440 {
2441 Array<OneD, NekDouble> tmpfieldcoeffs(m_fields[0]->GetNcoeffs() / 2);
2442 m_pressure->GetPlane(0)->BwdTrans(
2443 m_pressure->GetPlane(0)->GetCoeffs(),
2444 m_pressure->GetPlane(0)->UpdatePhys());
2445 m_pressure->GetPlane(1)->BwdTrans(
2446 m_pressure->GetPlane(1)->GetCoeffs(),
2447 m_pressure->GetPlane(1)->UpdatePhys());
2448 m_fields[0]->GetPlane(0)->FwdTransLocalElmt(
2449 m_pressure->GetPlane(0)->GetPhys(), fieldcoeffs[i]);
2450 m_fields[0]->GetPlane(1)->FwdTransLocalElmt(
2451 m_pressure->GetPlane(1)->GetPhys(), tmpfieldcoeffs);
2452 for (int e = 0; e < m_fields[0]->GetNcoeffs() / 2; e++)
2453 {
2454 fieldcoeffs[i][e + m_fields[0]->GetNcoeffs() / 2] =
2455 tmpfieldcoeffs[e];
2456 }
2457 }
2458 else
2459 {
2460 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
2461 m_fields[0]->FwdTransLocalElmt(m_pressure->GetPhys(), fieldcoeffs[i]);
2462 }
2463 variables[i] = "p";
2464
2465 if (!m_comm->IsParallelInTime())
2466 {
2467 WriteFld(m_sessionName + ".fld", m_fields[0], fieldcoeffs, variables);
2468 }
2469 else
2470 {
2471 std::string newdir = m_sessionName + ".pit";
2472 if (!fs::is_directory(newdir))
2473 {
2474 fs::create_directory(newdir);
2475 }
2476 WriteFld(
2477 newdir + "/" + m_sessionName + "_" +
2478 std::to_string(m_windowPIT * m_comm->GetTimeComm()->GetSize() +
2479 m_comm->GetTimeComm()->GetRank() + 1) +
2480 ".fld",
2481 m_fields[0], fieldcoeffs, variables);
2482 }
2483}
2484
2486{
2487 return m_session->GetVariables().size();
2488}
2489} // 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
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
void L2Norm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
CoupledLinearNS(const LibUtilities::SessionReaderSharedPtr &pSesssion, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
void v_DoSolve(void) override
Solves an unsteady problem.
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
void v_Output(void) override
void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
static std::string solverTypeLookupId
const SpatialDomains::ExpansionInfoMap & GenPressureExp(const SpatialDomains::ExpansionInfoMap &VelExp)
bool v_NegatedOp(void) override
void EvaluateAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
void v_TransPhysToCoeff(void) override
Virtual function for transformation to coefficient space.
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
static std::string className
Name of class.
int v_GetForceDimension(void) override
void InfNorm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
void v_TransCoeffToPhys(void) override
Virtual function for transformation to physical space.
Array< OneD, CoupledSolverMatrices > m_mat
This class is the base class for Navier Stokes problems.
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
NekDouble m_kinvis
Kinematic viscosity.
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);
EquationType m_equationType
equation type;
int m_nConvectiveFields
Number of fields to be convected;.
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)
Describes the specification for a Basis.
Definition: Basis.h:45
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:131
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:137
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:74
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
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.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:65
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
int m_windowPIT
Index of windows for parallel-in-time time iteration.
SOLVER_UTILS_EXPORT int GetNpoints()
NekDouble m_lambda
Lambda constant in real system if one required.
SOLVER_UTILS_EXPORT int GetNcoeffs()
int m_npointsZ
number of points in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
bool m_singleMode
Flag to determine if single homogeneous mode is used.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
enum HomogeneousType m_HomogeneousType
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
SOLVER_UTILS_EXPORT int GetTotPoints()
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
SOLVER_UTILS_EXPORT void v_DoSolve() override
Solves an unsteady problem.
static void Dcopy(const int &n, const double *x, const int &incx, double *y, const int &incy)
BLAS level 1: Copy x to y.
Definition: Blas.hpp:128
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
static const NekDouble kNekUnsetDouble
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
Definition: MeshGraph.h:143
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition: MeshGraph.h:140
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:141
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
StdRegions::ConstFactorMap factors
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
@ eSteadyNavierStokes
@ eUnsteadyStokes
@ eUnsteadyNavierStokes
@ eSteadyLinearisedNS
@ eNoEquationType
@ eCoupledLinearisedNS
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:47
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< CoupledLocalToGlobalC0ContMap > CoupledLocalToGlobalC0ContMapSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
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 Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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
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
DNekScalBlkMatSharedPtr m_D_int
Inner product of pressure system with divergence of the interior velocity space .
DNekScalBlkMatSharedPtr m_D_bnd
Inner product of pressure system with divergence of the boundary velocity space .
DNekScalBlkMatSharedPtr m_Cinv
Interior-Interior Laplaican plus linearised convective terms inverted, i.e. the inverse of .
DNekScalBlkMatSharedPtr m_Btilde
Interior-boundary Laplacian plus linearised convective terms .
DNekScalBlkMatSharedPtr m_BCinv
Boundary-interior Laplacian plus linearised convective terms pre-multiplying Cinv: .
MultiRegions::GlobalLinSysSharedPtr m_CoupledBndSys