Nektar++
Loading...
Searching...
No Matches
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
37#include <IncNavierStokesSolver/EquationSystems/CoupledLinearNS.h>
38#include <LibUtilities/BasicUtils/Timer.h>
39#include <LocalRegions/MatrixKey.h>
40#include <MultiRegions/ContField.h>
41#include <MultiRegions/GlobalLinSysDirectStaticCond.h>
42
43namespace Nektar
44{
45
46std::string CoupledLinearNS::className =
47 SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
48 "CoupledLinearisedNS", CoupledLinearNS::create);
49
50std::string CoupledLinearNS::solverTypeLookupId =
51 LibUtilities::SessionReader::RegisterEnumValue(
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 */
61CoupledLinearNS::CoupledLinearNS(
62 const LibUtilities::SessionReaderSharedPtr &pSession,
63 const SpatialDomains::MeshGraphSharedPtr &pGraph)
64 : UnsteadySystem(pSession, pGraph), IncNavierStokes(pSession, pGraph),
65 m_zeroMode(false)
66{
67}
68
69void CoupledLinearNS::v_InitObject(bool DeclareField)
70{
71 IncNavierStokes::v_InitObject(DeclareField);
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
80 m_nConvectiveFields = m_fields.size();
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
93 if (m_HomogeneousType == eHomogeneous1D)
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
100 m_pressure = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::
101 AllocateSharedPtr(m_session, Homo1DKey, m_LhomZ, m_useFFT,
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 =
111 MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
112 m_session, pressure_exp);
113 nz = 1;
114 }
115
116 Array<OneD, MultiRegions::ExpListSharedPtr> velocity(m_velocity.size());
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
123 m_locToGloMap = Array<OneD, CoupledLocalToGlobalC0ContMapSharedPtr>(nz);
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] =
136 MemoryManager<CoupledLocalToGlobalC0ContMap>::AllocateSharedPtr(
137 m_session, m_graph, m_boundaryConditions, velocity,
138 m_pressure, nz_loc, m_zeroMode);
139 }
140 else
141 {
142 // base mode
143 int nz_loc = 1;
144 m_locToGloMap[0] =
145 MemoryManager<CoupledLocalToGlobalC0ContMap>::AllocateSharedPtr(
146 m_session, m_graph, m_boundaryConditions, velocity,
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] =
155 MemoryManager<CoupledLocalToGlobalC0ContMap>::
156 AllocateSharedPtr(
157 m_session, m_graph, m_boundaryConditions, velocity,
158 m_pressure->GetPlane(2), nz_loc, false);
159 // Note high order modes cannot be singular
160 for (i = 2; i < nz; ++i)
161 {
162 m_locToGloMap[i] = m_locToGloMap[1];
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
179 if (m_equationType == eUnsteadyNavierStokes)
180 {
181 std::string vExtrapolation = "Standard";
182
183 if (m_session->DefinesSolverInfo("Extrapolation"))
184 {
185 vExtrapolation = m_session->GetSolverInfo("Extrapolation");
186 }
187
188 m_extrapolation = GetExtrapolateFactory().CreateInstance(
189 vExtrapolation, m_session, m_fields, m_pressure, m_velocity,
190 m_advObject);
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 */
200void CoupledLinearNS::SetUpCoupledMatrix(
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,
215 NekConstants::kNekUnsetDouble);
216 nz = 1;
217 m_mat = Array<OneD, CoupledSolverMatrices>(nz);
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
248 m_mat = Array<OneD, CoupledSolverMatrices>(nz);
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
402void CoupledLinearNS::SetUpCoupledMatrix(
403 const NekDouble lambda, const Array<OneD, Array<OneD, NekDouble>> &Advfield,
404 bool IsLinearNSEquation, const int HomogeneousMode,
405 CoupledSolverMatrices &mat,
406 CoupledLocalToGlobalC0ContMapSharedPtr &locToGloMap,
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;
423 StdRegions::StdExpansionSharedPtr locExp;
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;
469 DNekScalBlkMatSharedPtr pAh =
470 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
471 nsize_bndry_p1, nsize_bndry_p1, blkmatStorage);
472 mat.m_BCinv = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
473 nsize_bndry, nsize_int, blkmatStorage);
474 mat.m_Btilde = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
475 nsize_bndry, nsize_int, blkmatStorage);
476 mat.m_Cinv = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
477 nsize_int, nsize_int, blkmatStorage);
478
479 mat.m_D_bnd = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
480 nsize_p, nsize_bndry, blkmatStorage);
481
482 mat.m_D_int = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
483 nsize_p, nsize_int, blkmatStorage);
484
485 // Final level static condensation matrices.
486 DNekScalBlkMatSharedPtr pBh =
487 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
488 nsize_bndry_p1, nsize_p_m1, blkmatStorage);
489 DNekScalBlkMatSharedPtr pCh =
490 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
491 nsize_p_m1, nsize_bndry_p1, blkmatStorage);
492 DNekScalBlkMatSharedPtr pDh =
493 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nsize_p_m1, nsize_p_m1,
494 blkmatStorage);
495
496 LibUtilities::Timer timer;
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];
504 DNekMatSharedPtr Ah =
505 MemoryManager<DNekMat>::AllocateSharedPtr(k, k, zero);
506 Array<OneD, NekDouble> Ah_data = Ah->GetPtr();
507 int AhRows = k;
508 DNekMatSharedPtr B =
509 MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nint, zero);
510 Array<OneD, NekDouble> B_data = B->GetPtr();
511 DNekMatSharedPtr C =
512 MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nint, zero);
513 Array<OneD, NekDouble> C_data = C->GetPtr();
514 DNekMatSharedPtr D =
515 MemoryManager<DNekMat>::AllocateSharedPtr(nint, nint, zero);
516 Array<OneD, NekDouble> D_data = D->GetPtr();
517
518 DNekMatSharedPtr Dbnd = MemoryManager<DNekMat>::AllocateSharedPtr(
519 nsize_p[n], nsize_bndry[n], zero);
520 DNekMatSharedPtr Dint = MemoryManager<DNekMat>::AllocateSharedPtr(
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);
526 StdRegions::ConstFactorMap factors;
527 factors[StdRegions::eFactorLambda] = lambda / m_kinvis;
528 LocalRegions::MatrixKey helmkey(
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
549 DNekScalBlkMatSharedPtr CondMat;
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
615 Vmath::Smul(
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 {
647 Blas::Dcopy(
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
674 Vmath::Smul(
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 {
710 Blas::Dcopy(
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 {
740 LocalRegions::MatrixKey masskey(
741 StdRegions::eMass, locExp->DetShapeType(), *locExp);
742 MassMat = locExp->as<LocalRegions::Expansion>()->GetLocMatrix(
743 masskey);
744 }
745
746 Array<OneD, NekDouble> Advtmp;
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(
778 MultiRegions::DirCartesianMap[nv1],
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 {
950 Blas::Dcopy(
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 {
1149 Blas::Dcopy(
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,
1204 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, B));
1205 mat.m_Btilde->SetBlock(
1206 n, n,
1207 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, C));
1208 mat.m_Cinv->SetBlock(
1209 n, n,
1210 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, D));
1211 mat.m_D_bnd->SetBlock(
1212 n, n,
1213 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Dbnd));
1214 mat.m_D_int->SetBlock(
1215 n, n,
1216 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Dint));
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.
1236 DNekMatSharedPtr Bh = MemoryManager<DNekMat>::AllocateSharedPtr(
1237 nsize_bndry_p1[n], nsize_p_m1[n], zero);
1238 DNekMatSharedPtr Ch = MemoryManager<DNekMat>::AllocateSharedPtr(
1239 nsize_p_m1[n], nsize_bndry_p1[n], zero);
1240 DNekMatSharedPtr Dh = MemoryManager<DNekMat>::AllocateSharedPtr(
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,
1328 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Ah));
1329 pBh->SetBlock(
1330 n, n,
1331 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Bh));
1332 pCh->SetBlock(
1333 n, n,
1334 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Ch));
1335 pDh->SetBlock(
1336 n, n,
1337 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Dh));
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
1347 MultiRegions::GlobalLinSysKey key(StdRegions::eLinearAdvectionReaction,
1348 locToGloMap);
1349 mat.m_CoupledBndSys =
1350 MemoryManager<MultiRegions::GlobalLinSysDirectStaticCond>::
1351 AllocateSharedPtr(key, m_fields[0], pAh, pBh, pCh, pDh,
1352 locToGloMap);
1353 mat.m_CoupledBndSys->Initialise(locToGloMap);
1354}
1355
1356void CoupledLinearNS::v_GenerateSummary(SolverUtils::SummaryList &s)
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:
1366 case eUnsteadyNavierStokes:
1367 {
1368 // Could defind this from IncNavierStokes class?
1369 m_ode.DefineOdeRhs(&CoupledLinearNS::EvaluateAdvection, this);
1370
1371 m_ode.DefineImplicitSolve(
1372 &CoupledLinearNS::SolveUnsteadyStokesSystem, this);
1373
1374 // Set initial condition using time t=0
1375
1376 SetInitialConditions(0.0, dumpInitialConditions);
1377 break;
1378 }
1379 case eSteadyStokes:
1380 SetUpCoupledMatrix(0.0);
1381 break;
1382 case eSteadyOseen:
1383 {
1384 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1385 for (size_t i = 0; i < m_velocity.size(); ++i)
1386 {
1387 AdvField[i] = Array<OneD, NekDouble>(
1388 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1389 }
1390
1391 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1392 "Advection Velocity section must be defined in "
1393 "session file.");
1394
1395 std::vector<std::string> fieldStr;
1396 for (size_t i = 0; i < m_velocity.size(); ++i)
1397 {
1398 fieldStr.push_back(
1399 m_boundaryConditions->GetVariable(m_velocity[i]));
1400 }
1401 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1402
1403 SetUpCoupledMatrix(0.0, AdvField, false);
1404 }
1405 break;
1406 case eSteadyNavierStokes:
1407 {
1408 m_session->LoadParameter("KinvisMin", m_kinvisMin);
1409 m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1410 m_session->LoadParameter("Tolerence", m_tol);
1411 m_session->LoadParameter("MaxIteration", m_maxIt);
1412 m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1413 m_session->LoadParameter("Restart", m_Restart);
1414
1415 DefineForcingTerm();
1416
1417 if (m_Restart == 1)
1418 {
1419 ASSERTL0(m_session->DefinesFunction("Restart"),
1420 "Restart section must be defined in session file.");
1421
1422 Array<OneD, Array<OneD, NekDouble>> Restart(m_velocity.size());
1423 for (size_t i = 0; i < m_velocity.size(); ++i)
1424 {
1425 Restart[i] = Array<OneD, NekDouble>(
1426 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1427 }
1428 std::vector<std::string> fieldStr;
1429 for (size_t i = 0; i < m_velocity.size(); ++i)
1430 {
1431 fieldStr.push_back(
1432 m_boundaryConditions->GetVariable(m_velocity[i]));
1433 }
1434 GetFunction("Restart")->Evaluate(fieldStr, Restart);
1435
1436 for (size_t i = 0; i < m_velocity.size(); ++i)
1437 {
1438 m_fields[m_velocity[i]]->FwdTransLocalElmt(
1439 Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1440 }
1441 std::cout << "Saving the RESTART file for m_kinvis = "
1442 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1443 << std::endl;
1444 }
1445 else // We solve the Stokes Problem
1446 {
1447
1448 SetUpCoupledMatrix(0.0);
1449 m_initialStep = true;
1450 m_counter = 1;
1451 // SolveLinearNS(m_ForcingTerm_Coeffs);
1452 Solve();
1453 m_initialStep = false;
1454 std::cout << "Saving the Stokes Flow for m_kinvis = "
1455 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1456 << std::endl;
1457 }
1458 }
1459 break;
1460 case eSteadyLinearisedNS:
1461 {
1462 SetInitialConditions(0.0, dumpInitialConditions);
1463
1464 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1465 for (size_t i = 0; i < m_velocity.size(); ++i)
1466 {
1467 AdvField[i] = Array<OneD, NekDouble>(
1468 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1469 }
1470
1471 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1472 "Advection Velocity section must be defined in "
1473 "session file.");
1474
1475 std::vector<std::string> fieldStr;
1476 for (size_t i = 0; i < m_velocity.size(); ++i)
1477 {
1478 fieldStr.push_back(
1479 m_boundaryConditions->GetVariable(m_velocity[i]));
1480 }
1481 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1482
1483 SetUpCoupledMatrix(m_lambda, AdvField, true);
1484 }
1485 break;
1486 case eNoEquationType:
1487 default:
1488 ASSERTL0(false,
1489 "Unknown or undefined equation type for CoupledLinearNS");
1490 }
1491}
1492
1493void CoupledLinearNS::EvaluateAdvection(
1494 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1495 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
1496{
1497 // evaluate convection terms
1498 EvaluateAdvectionTerms(inarray, outarray, time);
1499
1500 for (auto &x : m_forcing)
1501 {
1502 x->Apply(m_fields, outarray, outarray, time);
1503 }
1504}
1505
1506void CoupledLinearNS::SolveUnsteadyStokesSystem(
1507 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1508 Array<OneD, Array<OneD, NekDouble>> &outarray,
1509 [[maybe_unused]] const NekDouble time, const NekDouble aii_Dt)
1510{
1511 size_t i;
1512 Array<OneD, Array<OneD, NekDouble>> F(m_nConvectiveFields);
1513 NekDouble lambda = 1.0 / aii_Dt;
1514 static NekDouble lambda_store;
1515 Array<OneD, Array<OneD, NekDouble>> forcing(m_velocity.size());
1516 // Matrix solution
1517 if (fabs(lambda_store - lambda) > 1e-10)
1518 {
1519 SetUpCoupledMatrix(lambda);
1520 lambda_store = lambda;
1521 }
1522
1523 // Forcing for advection solve
1524 for (i = 0; i < m_velocity.size(); ++i)
1525 {
1526 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1527 m_fields[m_velocity[i]]->SetWaveSpace(true);
1528 m_fields[m_velocity[i]]->IProductWRTBase(
1529 inarray[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1530 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1531 Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(), lambda,
1532 m_fields[m_velocity[i]]->GetCoeffs(), 1,
1533 m_fields[m_velocity[i]]->UpdateCoeffs(), 1);
1534 forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1535 }
1536
1537 SolveLinearNS(forcing);
1538
1539 for (i = 0; i < m_velocity.size(); ++i)
1540 {
1541 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1542 outarray[i]);
1543 }
1544}
1545
1546void CoupledLinearNS::v_TransCoeffToPhys(void)
1547{
1548 size_t nfields = m_fields.size();
1549 for (size_t k = 0; k < nfields; ++k)
1550 {
1551 // Backward Transformation in physical space for time evolution
1552 m_fields[k]->BwdTrans(m_fields[k]->GetCoeffs(),
1553 m_fields[k]->UpdatePhys());
1554 }
1555}
1556
1557void CoupledLinearNS::v_TransPhysToCoeff(void)
1558{
1559 size_t nfields = m_fields.size();
1560 for (size_t k = 0; k < nfields; ++k)
1561 {
1562 // Forward Transformation in physical space for time evolution
1563 m_fields[k]->FwdTransLocalElmt(m_fields[k]->GetPhys(),
1564 m_fields[k]->UpdateCoeffs());
1565 }
1566}
1567
1568void CoupledLinearNS::v_DoSolve(void)
1569{
1570 switch (m_equationType)
1571 {
1572 case eUnsteadyStokes:
1573 case eUnsteadyNavierStokes:
1574 // AdvanceInTime(m_steps);
1575 UnsteadySystem::v_DoSolve();
1576 break;
1577 case eSteadyStokes:
1578 case eSteadyOseen:
1579 case eSteadyLinearisedNS:
1580 {
1581 Solve();
1582 break;
1583 }
1584 case eSteadyNavierStokes:
1585 {
1586 LibUtilities::Timer Generaltimer;
1587 Generaltimer.Start();
1588
1589 int Check(0);
1590
1591 // Saving the init datas
1592 Checkpoint_Output(Check);
1593 Check++;
1594
1595 std::cout
1596 << "We execute INITIALLY SolveSteadyNavierStokes for m_kinvis "
1597 "= "
1598 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1599 << std::endl;
1600 SolveSteadyNavierStokes();
1601
1602 while (m_kinvis > m_kinvisMin)
1603 {
1604 if (Check == 1)
1605 {
1606 std::cout
1607 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1608 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1609 << std::endl;
1610 SolveSteadyNavierStokes();
1611 Checkpoint_Output(Check);
1612 Check++;
1613 }
1614
1615 Continuation();
1616
1617 if (m_kinvis > m_kinvisMin)
1618 {
1619 std::cout
1620 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1621 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1622 << std::endl;
1623 SolveSteadyNavierStokes();
1624 Checkpoint_Output(Check);
1625 Check++;
1626 }
1627 }
1628
1629 Generaltimer.Stop();
1630 std::cout << "\nThe total calculation time is : "
1631 << Generaltimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1632
1633 break;
1634 }
1635 case eNoEquationType:
1636 default:
1637 ASSERTL0(false,
1638 "Unknown or undefined equation type for CoupledLinearNS");
1639 }
1640}
1641
1642/** Virtual function to define if operator in DoSolve is negated
1643 * with regard to the strong form. This is currently only used in
1644 * Arnoldi solves. For Coupledd solver this is true since Stokes
1645 * operator is set up as a LHS rather than RHS operation
1646 */
1647bool CoupledLinearNS::v_NegatedOp(void)
1648{
1649 return true;
1650}
1651
1652void CoupledLinearNS::Solve(void)
1653{
1654 const size_t ncmpt = m_velocity.size();
1655 Array<OneD, Array<OneD, NekDouble>> forcing_phys(ncmpt);
1656 Array<OneD, Array<OneD, NekDouble>> forcing(ncmpt);
1657
1658 for (size_t i = 0; i < ncmpt; ++i)
1659 {
1660 forcing_phys[i] =
1661 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1662 forcing[i] =
1663 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNcoeffs(), 0.0);
1664 }
1665
1666 for (auto &x : m_forcing)
1667 {
1668 const NekDouble time = 0;
1669 x->Apply(m_fields, forcing_phys, forcing_phys, time);
1670 }
1671 for (size_t i = 0; i < ncmpt; ++i)
1672 {
1673 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1674 m_fields[i]->SetWaveSpace(true);
1675 m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1676 m_fields[i]->SetWaveSpace(waveSpace);
1677 }
1678
1679 SolveLinearNS(forcing);
1680}
1681
1682void CoupledLinearNS::DefineForcingTerm(void)
1683{
1684 m_ForcingTerm = Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1685 m_ForcingTerm_Coeffs =
1686 Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1687
1688 for (size_t i = 0; i < m_velocity.size(); ++i)
1689 {
1690 m_ForcingTerm[i] = Array<OneD, NekDouble>(
1691 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1692 m_ForcingTerm_Coeffs[i] =
1693 Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1694 }
1695
1696 if (m_session->DefinesFunction("ForcingTerm"))
1697 {
1698 std::vector<std::string> fieldStr;
1699 for (size_t i = 0; i < m_velocity.size(); ++i)
1700 {
1701 fieldStr.push_back(
1702 m_boundaryConditions->GetVariable(m_velocity[i]));
1703 }
1704 GetFunction("ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1705 for (size_t i = 0; i < m_velocity.size(); ++i)
1706 {
1707 m_fields[m_velocity[i]]->FwdTransLocalElmt(m_ForcingTerm[i],
1708 m_ForcingTerm_Coeffs[i]);
1709 }
1710 }
1711 else
1712 {
1713 std::cout
1714 << "'ForcingTerm' section has not been defined in the input file "
1715 "=> forcing=0"
1716 << std::endl;
1717 }
1718}
1719
1720void CoupledLinearNS::SolveSteadyNavierStokes(void)
1721{
1722 LibUtilities::Timer Newtontimer;
1723 Newtontimer.Start();
1724
1725 Array<OneD, Array<OneD, NekDouble>> RHS_Coeffs(m_velocity.size());
1726 Array<OneD, Array<OneD, NekDouble>> RHS_Phys(m_velocity.size());
1727 Array<OneD, Array<OneD, NekDouble>> delta_velocity_Phys(m_velocity.size());
1728 Array<OneD, Array<OneD, NekDouble>> Velocity_Phys(m_velocity.size());
1729 Array<OneD, NekDouble> L2_norm(m_velocity.size(), 1.0);
1730 Array<OneD, NekDouble> Inf_norm(m_velocity.size(), 1.0);
1731
1732 for (size_t i = 0; i < m_velocity.size(); ++i)
1733 {
1734 delta_velocity_Phys[i] = Array<OneD, NekDouble>(
1735 m_fields[m_velocity[i]]->GetTotPoints(), 1.0);
1736 Velocity_Phys[i] = Array<OneD, NekDouble>(
1737 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1738 }
1739
1740 m_counter = 1;
1741
1742 L2Norm(delta_velocity_Phys, L2_norm);
1743
1744 // while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1745 while (std::max(L2_norm[0], L2_norm[1]) > m_tol)
1746 {
1747 if (m_counter == 1)
1748 // At the first Newton step, we use the solution of the
1749 // Stokes problem (if Restart=0 in input file) Or the
1750 // solution of the .rst file (if Restart=1 in input
1751 // file)
1752 {
1753 for (size_t i = 0; i < m_velocity.size(); ++i)
1754 {
1755 RHS_Coeffs[i] = Array<OneD, NekDouble>(
1756 m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1757 RHS_Phys[i] = Array<OneD, NekDouble>(
1758 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1759 }
1760
1761 for (size_t i = 0; i < m_velocity.size(); ++i)
1762 {
1763 m_fields[m_velocity[i]]->BwdTrans(
1764 m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1765 }
1766
1767 m_initialStep = true;
1768 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1769 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1770 SolveLinearNS(RHS_Coeffs);
1771 m_initialStep = false;
1772 }
1773 if (m_counter > 1)
1774 {
1775 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1776
1777 if (m_counter % m_MatrixSetUpStep ==
1778 0) // Setting Up the matrix is expensive. We do it at each
1779 // "m_MatrixSetUpStep" step.
1780 {
1781 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1782 }
1783 SolveLinearNS(RHS_Coeffs);
1784 }
1785
1786 for (size_t i = 0; i < m_velocity.size(); ++i)
1787 {
1788 m_fields[m_velocity[i]]->BwdTrans(RHS_Coeffs[i], RHS_Phys[i]);
1789 m_fields[m_velocity[i]]->BwdTrans(
1790 m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1791 }
1792
1793 for (size_t i = 0; i < m_velocity.size(); ++i)
1794 {
1795 Vmath::Vadd(Velocity_Phys[i].size(), Velocity_Phys[i], 1,
1796 delta_velocity_Phys[i], 1, Velocity_Phys[i], 1);
1797 }
1798
1799 // InfNorm(delta_velocity_Phys, Inf_norm);
1800 L2Norm(delta_velocity_Phys, L2_norm);
1801
1802 if (std::max(Inf_norm[0], Inf_norm[1]) > 100)
1803 {
1804 std::cout << "\nThe Newton method has failed at m_kinvis = "
1805 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1806 << std::endl;
1807 ASSERTL0(0, "The Newton method has failed... \n");
1808 }
1809
1810 std::cout << "\n";
1811 m_counter++;
1812 }
1813
1814 if (m_counter > 1) // We save u:=u+\delta u in u->Coeffs
1815 {
1816 for (size_t i = 0; i < m_velocity.size(); ++i)
1817 {
1818 m_fields[m_velocity[i]]->FwdTrans(
1819 Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1820 }
1821 }
1822
1823 Newtontimer.Stop();
1824 std::cout << "We have done " << m_counter - 1 << " iteration(s) in "
1825 << Newtontimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1826}
1827
1828void CoupledLinearNS::Continuation(void)
1829{
1830 Array<OneD, Array<OneD, NekDouble>> u_N(m_velocity.size());
1831 Array<OneD, Array<OneD, NekDouble>> tmp_RHS(m_velocity.size());
1832 Array<OneD, Array<OneD, NekDouble>> RHS(m_velocity.size());
1833 Array<OneD, Array<OneD, NekDouble>> u_star(m_velocity.size());
1834
1835 std::cout << "We apply the continuation method: " << std::endl;
1836
1837 for (size_t i = 0; i < m_velocity.size(); ++i)
1838 {
1839 u_N[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1840 0.0);
1841 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1842 u_N[i]);
1843
1844 RHS[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1845 0.0);
1846 tmp_RHS[i] = Array<OneD, NekDouble>(
1847 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1848
1849 m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1850 Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1851
1852 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1853 m_fields[m_velocity[i]]->SetWaveSpace(true);
1854 m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1855 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1856 }
1857
1858 SetUpCoupledMatrix(0.0, u_N, true);
1859 SolveLinearNS(RHS);
1860
1861 for (size_t i = 0; i < m_velocity.size(); ++i)
1862 {
1863 u_star[i] = Array<OneD, NekDouble>(
1864 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1865 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1866 u_star[i]);
1867
1868 // u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1869 Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1870 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1871
1872 m_fields[m_velocity[i]]->FwdTrans(
1873 u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1874 }
1875
1876 m_kinvis -= m_kinvis * m_KinvisPercentage / 100;
1877}
1878
1879void CoupledLinearNS::InfNorm(Array<OneD, Array<OneD, NekDouble>> &inarray,
1880 Array<OneD, NekDouble> &outarray)
1881{
1882 for (size_t i = 0; i < m_velocity.size(); ++i)
1883 {
1884 outarray[i] = 0.0;
1885 for (size_t j = 0; j < inarray[i].size(); ++j)
1886 {
1887 if (inarray[i][j] > outarray[i])
1888 {
1889 outarray[i] = inarray[i][j];
1890 }
1891 }
1892 std::cout << "InfNorm[" << i << "] = " << outarray[i] << std::endl;
1893 }
1894}
1895
1896void CoupledLinearNS::L2Norm(Array<OneD, Array<OneD, NekDouble>> &inarray,
1897 Array<OneD, NekDouble> &outarray)
1898{
1899 for (size_t i = 0; i < m_velocity.size(); ++i)
1900 {
1901 outarray[i] = 0.0;
1902 for (size_t j = 0; j < inarray[i].size(); ++j)
1903 {
1904 outarray[i] += inarray[i][j] * inarray[i][j];
1905 }
1906 outarray[i] = std::sqrt(outarray[i]);
1907 std::cout << "L2Norm[" << i << "] = " << outarray[i] << std::endl;
1908 }
1909}
1910
1911void CoupledLinearNS::EvaluateNewtonRHS(
1912 Array<OneD, Array<OneD, NekDouble>> &Velocity,
1913 Array<OneD, Array<OneD, NekDouble>> &outarray)
1914{
1915 Array<OneD, Array<OneD, NekDouble>> Eval_Adv(m_velocity.size());
1916 Array<OneD, Array<OneD, NekDouble>> tmp_DerVel(m_velocity.size());
1917 Array<OneD, Array<OneD, NekDouble>> AdvTerm(m_velocity.size());
1918 Array<OneD, Array<OneD, NekDouble>> ViscTerm(m_velocity.size());
1919 Array<OneD, Array<OneD, NekDouble>> Forc(m_velocity.size());
1920
1921 for (size_t i = 0; i < m_velocity.size(); ++i)
1922 {
1923 Eval_Adv[i] = Array<OneD, NekDouble>(
1924 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1925 tmp_DerVel[i] = Array<OneD, NekDouble>(
1926 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1927
1928 AdvTerm[i] = Array<OneD, NekDouble>(
1929 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1930 ViscTerm[i] = Array<OneD, NekDouble>(
1931 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1932 Forc[i] = Array<OneD, NekDouble>(
1933 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1934 outarray[i] = Array<OneD, NekDouble>(
1935 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1936
1937 m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1938
1939 Vmath::Smul(tmp_DerVel[i].size(), m_kinvis, tmp_DerVel[i], 1,
1940 tmp_DerVel[i], 1);
1941 }
1942
1943 EvaluateAdvectionTerms(Velocity, Eval_Adv, m_time);
1944
1945 for (size_t i = 0; i < m_velocity.size(); ++i)
1946 {
1947 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1948 m_fields[m_velocity[i]]->SetWaveSpace(true);
1949 m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i],
1950 AdvTerm[i]); //(w, (u.grad)u)
1951 m_fields[m_velocity[i]]->IProductWRTDerivBase(
1952 i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1953 m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i],
1954 Forc[i]); //(w, f)
1955 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1956
1957 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1958 outarray[i], 1);
1959 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1960 outarray[i], 1);
1961
1962 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1963 1);
1964 }
1965}
1966
1967const SpatialDomains::ExpansionInfoMap &CoupledLinearNS::GenPressureExp(
1968 const SpatialDomains::ExpansionInfoMap &VelExp)
1969{
1970 SpatialDomains::ExpansionInfoMapShPtr returnval;
1971 returnval =
1972 MemoryManager<SpatialDomains::ExpansionInfoMap>::AllocateSharedPtr();
1973
1974 int nummodes;
1975
1976 for (auto &expMapIter : VelExp)
1977 {
1978 LibUtilities::BasisKeyVector BasisVec;
1979
1980 for (size_t i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
1981 {
1982 LibUtilities::BasisKey B = expMapIter.second->m_basisKeyVector[i];
1983 nummodes = B.GetNumModes();
1984 ASSERTL0(nummodes > 3,
1985 "Velocity polynomial space not sufficiently high (>= 4)");
1986 // Should probably set to be an orthogonal basis.
1987 LibUtilities::BasisKey newB(B.GetBasisType(), nummodes - 2,
1988 B.GetPointsKey());
1989 BasisVec.push_back(newB);
1990 }
1991
1992 // Put new expansion into list.
1993 SpatialDomains::ExpansionInfoShPtr expansionElementShPtr =
1994 MemoryManager<SpatialDomains::ExpansionInfo>::AllocateSharedPtr(
1995 expMapIter.second->m_geomPtr, BasisVec);
1996 (*returnval)[expMapIter.first] = expansionElementShPtr;
1997 }
1998
1999 // Save expansion into graph.
2000 m_graph->SetExpansionInfo("p", returnval);
2001
2002 return *returnval;
2003}
2004
2005/**
2006 * @param forcing A list of forcing functions for each velocity
2007 * component
2008 *
2009 * The routine involves two levels of static
2010 * condensations. Initially we require a statically condensed
2011 * forcing function which requires the following manipulation
2012 *
2013 * \f[ {F\_bnd} = {\bf f}_{bnd} -m\_B \,m\_Cinv\, {\bf f}_{int},
2014 * \hspace{1cm} F\_p = m\_D\_{int}\, m\_Cinv\, {\bf f}_{int} \f]
2015 *
2016 * Where \f${\bf f}_{bnd}\f$ denote the forcing degrees of
2017 * freedom of the elemental velocities on the boundary of the
2018 * element, \f${\bf f}_{int}\f$ denote the forcing degrees of
2019 * freedom of the elemental velocities on the interior of the
2020 * element. (see detailed description for more details).
2021 *
2022 * This vector is further manipulated into
2023 *
2024 * \f[ Fh\_{bnd} = \left [ \begin{array}{c} f\_{bnd} -m\_B \,
2025 * m\_Cinv\, {\bf f}_{int}\\ \left [m\_D\_{int} \, m\_Cinv \,{\bf
2026 * f}_{int} \right]_0 \end{array}\right ] \hspace{1cm} [Fh\_p]_{i} =
2027 * \begin{array}{c} [m\_D\_{int} \, m\_Cinv \, {\bf
2028 * f}_{int}]_{i+1} \end{array} \f]
2029 *
2030 * where \f$-{[}m\_D\_{int}^T\, m\_Cinv \,{\bf f}_{int}]_0\f$
2031 * which is corresponds to the mean mode of the pressure degrees
2032 * of freedom and is now added to the boundary system and the
2033 * remainder of the block becomes the interior forcing for the
2034 * inner static condensation (see detailed description for more
2035 * details) which is set up in a GlobalLinSysDirectStaticCond
2036 * class.
2037 *
2038 * Finally we perform the final maniplation of the forcing to
2039 * using hte
2040 * \f[ Fh\_{bnd} = Fh\_{bnd} - m\_Bh \,m\_Chinv \, Fh\_p \f]
2041 *
2042 * We can now call the solver to the global coupled boundary
2043 * system (through the call to #m_CoupledBndSys->Solve) to obtain
2044 * the velocity boundary solution as the mean pressure solution,
2045 * i.e.
2046 *
2047 * \f[ {\cal A}^T(\hat{A} - \hat{C}^T \hat{D}^{-1} \hat{B} ){\cal
2048 * A} \, Bnd = Fh\_{bnd} \f]
2049 *
2050 * Once we know the solution to the above the rest of the pressure
2051 * modes are recoverable thorugh
2052 *
2053 * \f[ Ph = m\_Dhinv\, (Bnd - m\_Ch^T \, Fh_{bnd}) \f]
2054 *
2055 * We can now unpack \f$ Fh\_{bnd} \f$ (last elemental mode) and
2056 * \f$ Ph \f$ into #m_pressure and \f$ F_p\f$ and \f$ Fh\_{bnd}\f$
2057 * into a closed pack list of boundary velocoity degrees of
2058 * freedom stored in \f$ F\_bnd\f$.
2059 *
2060 * Finally the interior velocity degrees of freedom are then
2061 * obtained through the relationship
2062 *
2063 * \f[ F\_{int} = m\_Cinv\ ( F\_{int} + m\_D\_{int}^T\,
2064 * F\_p - m\_Btilde^T\, Bnd) \f]
2065 *
2066 * We then unpack the solution back to the MultiRegion structures
2067 * #m_velocity and #m_pressure
2068 */
2069void CoupledLinearNS::SolveLinearNS(
2070 const Array<OneD, Array<OneD, NekDouble>> &forcing)
2071{
2072 size_t i;
2073 Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.size());
2074 Array<OneD, Array<OneD, NekDouble>> force(m_velocity.size());
2075
2076 // Impose Dirichlet conditions on velocity fields
2077 for (i = 0; i < m_velocity.size(); ++i)
2078 {
2079 Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
2080 m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
2081 }
2082
2083 if (m_HomogeneousType == eHomogeneous1D)
2084 {
2085 int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
2086 for (int n = 0; n < m_npointsZ / 2; ++n)
2087 {
2088 // Get the a Fourier mode of velocity and forcing components.
2089 for (i = 0; i < m_velocity.size(); ++i)
2090 {
2091 vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2 * n);
2092 // Note this needs to correlate with how we pass forcing
2093 force[i] = forcing[i] + 2 * n * ncoeffsplane;
2094 }
2095
2096 SolveLinearNS(force, vel_fields, m_pressure->GetPlane(2 * n), n);
2097 }
2098 for (i = 0; i < m_velocity.size(); ++i)
2099 {
2100 m_fields[m_velocity[i]]->SetPhysState(false);
2101 }
2102 m_pressure->SetPhysState(false);
2103 }
2104 else
2105 {
2106 for (i = 0; i < m_velocity.size(); ++i)
2107 {
2108 vel_fields[i] = m_fields[m_velocity[i]];
2109 // Note this needs to correlate with how we pass forcing
2110 force[i] = forcing[i];
2111 }
2112 SolveLinearNS(force, vel_fields, m_pressure);
2113 }
2114}
2115
2116void CoupledLinearNS::SolveLinearNS(
2117 Array<OneD, Array<OneD, NekDouble>> &forcing,
2118 Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
2119 MultiRegions::ExpListSharedPtr &pressure, const int mode)
2120{
2121 size_t i, j, k, n, cnt, cnt1;
2122 size_t nbnd, nint, offset;
2123 size_t nvel = m_velocity.size();
2124 size_t nel = fields[0]->GetNumElmts();
2125 Array<OneD, unsigned int> bmap, imap;
2126
2127 Array<OneD, NekDouble> f_bnd(m_mat[mode].m_BCinv->GetRows());
2128 NekVector<NekDouble> F_bnd(f_bnd.size(), f_bnd, eWrapper);
2129 Array<OneD, NekDouble> f_int(m_mat[mode].m_BCinv->GetColumns());
2130 NekVector<NekDouble> F_int(f_int.size(), f_int, eWrapper);
2131
2132 size_t nz_loc;
2133 int nplanecoeffs =
2134 fields[m_velocity[0]]
2135 ->GetNcoeffs(); // this is fine since we pass the nplane coeff data.
2136
2137 if (mode) // Homogeneous mode flag
2138 {
2139 nz_loc = 2;
2140 }
2141 else
2142 {
2143 if (m_singleMode)
2144 {
2145 nz_loc = 2;
2146 }
2147 else
2148 {
2149 nz_loc = 1;
2150 if (m_HomogeneousType == eHomogeneous1D)
2151 {
2152 Array<OneD, NekDouble> tmp;
2153 // Zero fields to set complex mode to zero;
2154 for (i = 0; i < fields.size(); ++i)
2155 {
2156 Vmath::Zero(nplanecoeffs,
2157 tmp = fields[i]->UpdateCoeffs() + nplanecoeffs,
2158 1);
2159 }
2160 Vmath::Zero(2 * pressure->GetNcoeffs(),
2161 pressure->UpdateCoeffs(), 1);
2162 }
2163 }
2164 }
2165
2166 for (k = 0; k < nvel; ++k)
2167 {
2168 MultiRegions::ContFieldSharedPtr cfield =
2169 std::dynamic_pointer_cast<MultiRegions::ContField>(fields[k]);
2170
2171 Array<OneD, NekDouble> sign =
2172 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsSign();
2173 const Array<OneD, const int> map =
2174 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsMap();
2175
2176 // Add weak boundary conditions to forcing
2177 const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConds =
2178 fields[k]->GetBndConditions();
2179 Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp;
2180
2181 if (m_HomogeneousType == eHomogeneous1D)
2182 {
2183 bndCondExp =
2184 m_fields[k]->GetPlane(2 * mode)->GetBndCondExpansions();
2185 }
2186 else
2187 {
2188 bndCondExp = m_fields[k]->GetBndCondExpansions();
2189 }
2190
2191 for (n = 0; n < nz_loc; ++n)
2192 {
2193 int bndcnt = 0;
2194 for (i = 0; i < bndCondExp.size(); ++i)
2195 {
2196 const Array<OneD, const NekDouble> bndcoeffs =
2197 bndCondExp[i]->GetCoeffs();
2198 size_t nCoeffs = (bndCondExp[i])->GetNcoeffs();
2199
2200 cnt = 0;
2201 if (bndConds[i]->GetBoundaryConditionType() ==
2202 SpatialDomains::eNeumann ||
2203 bndConds[i]->GetBoundaryConditionType() ==
2204 SpatialDomains::eRobin)
2205 {
2206 if (m_locToGloMap[mode]->GetSignChange())
2207 {
2208 for (j = 0; j < nCoeffs; j++)
2209 {
2210 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2211 sign[bndcnt + j] * bndcoeffs[j];
2212 }
2213 }
2214 else
2215 {
2216 for (j = 0; j < nCoeffs; j++)
2217 {
2218 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2219 bndcoeffs[j];
2220 }
2221 }
2222 }
2223
2224 bndcnt += bndCondExp[i]->GetNcoeffs();
2225 }
2226 }
2227 }
2228
2229 Array<OneD, NekDouble> bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(), 0.0);
2230
2231 // Construct f_bnd and f_int and fill in bnd from inarray
2232 // (with Dirichlet BCs imposed)
2233 int bndoffset = 0;
2234 cnt = cnt1 = 0;
2235 for (i = 0; i < nel; ++i) // loop over elements
2236 {
2237 fields[m_velocity[0]]->GetExp(i)->GetBoundaryMap(bmap);
2238 fields[m_velocity[0]]->GetExp(i)->GetInteriorMap(imap);
2239 nbnd = bmap.size();
2240 nint = imap.size();
2241 offset = fields[m_velocity[0]]->GetCoeff_Offset(i);
2242
2243 for (j = 0; j < nvel; ++j) // loop over velocity fields
2244 {
2245 Array<OneD, NekDouble> incoeffs = fields[j]->UpdateCoeffs();
2246
2247 for (n = 0; n < nz_loc; ++n)
2248 {
2249 for (k = 0; k < nbnd; ++k)
2250 {
2251 f_bnd[cnt + k] =
2252 forcing[j][n * nplanecoeffs + offset + bmap[k]];
2253
2254 bnd[bndoffset + (n + j * nz_loc) * nbnd + k] =
2255 incoeffs[n * nplanecoeffs + offset + bmap[k]];
2256 }
2257 for (k = 0; k < nint; ++k)
2258 {
2259 f_int[cnt1 + k] =
2260 forcing[j][n * nplanecoeffs + offset + imap[k]];
2261 }
2262
2263 cnt += nbnd;
2264 cnt1 += nint;
2265 }
2266 }
2267 bndoffset +=
2268 nvel * nz_loc * nbnd + nz_loc * (pressure->GetExp(i)->GetNcoeffs());
2269 }
2270
2271 Array<OneD, NekDouble> f_p(m_mat[mode].m_D_int->GetRows());
2272 NekVector<NekDouble> F_p(f_p.size(), f_p, eWrapper);
2273 NekVector<NekDouble> F_p_tmp(m_mat[mode].m_Cinv->GetRows());
2274
2275 // fbnd does not currently hold the pressure mean
2276 F_bnd = F_bnd - (*m_mat[mode].m_BCinv) * F_int;
2277 F_p_tmp = (*m_mat[mode].m_Cinv) * F_int;
2278 F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
2279
2280 // construct inner forcing
2281 Array<OneD, NekDouble> fh_bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(),
2282 0.0);
2283
2284 offset = cnt = 0;
2285 for (i = 0; i < nel; ++i)
2286 {
2287 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2288
2289 for (j = 0; j < nvel; ++j)
2290 {
2291 for (k = 0; k < nbnd; ++k)
2292 {
2293 fh_bnd[offset + j * nbnd + k] = f_bnd[cnt + k];
2294 }
2295 cnt += nbnd;
2296 }
2297
2298 nint = pressure->GetExp(i)->GetNcoeffs();
2299 offset += nvel * nbnd + nint * nz_loc;
2300 }
2301
2302 offset = cnt1 = 0;
2303 for (i = 0; i < nel; ++i)
2304 {
2305 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2306 nint = pressure->GetExp(i)->GetNcoeffs();
2307
2308 for (n = 0; n < nz_loc; ++n)
2309 {
2310 for (j = 0; j < nint; ++j)
2311 {
2312 fh_bnd[offset + nvel * nbnd + n * nint + j] = f_p[cnt1 + j];
2313 }
2314 cnt1 += nint;
2315 }
2316 offset += nvel * nbnd + nz_loc * nint;
2317 }
2318 m_mat[mode].m_CoupledBndSys->Solve(fh_bnd, bnd, m_locToGloMap[mode]);
2319
2320 // unpack pressure and velocity boundary systems.
2321 offset = cnt = 0;
2322 int totpcoeffs = pressure->GetNcoeffs();
2323 Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2324 for (i = 0; i < nel; ++i)
2325 {
2326 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2327 nint = pressure->GetExp(i)->GetNcoeffs();
2328 for (j = 0; j < nvel; ++j)
2329 {
2330 for (k = 0; k < nbnd; ++k)
2331 {
2332 f_bnd[cnt + k] = bnd[offset + j * nbnd + k];
2333 }
2334 cnt += nbnd;
2335 }
2336 offset += nvel * nbnd + nint * nz_loc;
2337 }
2338
2339 pressure->SetPhysState(false);
2340
2341 offset = cnt = cnt1 = 0;
2342 for (i = 0; i < nel; ++i)
2343 {
2344 nint = pressure->GetExp(i)->GetNcoeffs();
2345 nbnd = fields[0]->GetExp(i)->NumBndryCoeffs();
2346 cnt1 = pressure->GetCoeff_Offset(i);
2347
2348 for (n = 0; n < nz_loc; ++n)
2349 {
2350 for (j = 0; j < nint; ++j)
2351 {
2352 p_coeffs[n * totpcoeffs + cnt1 + j] = f_p[cnt + j] =
2353 bnd[offset + nvel * nz_loc * nbnd + n * nint + j];
2354 }
2355 cnt += nint;
2356 }
2357 offset += (nvel * nbnd + nint) * nz_loc;
2358 }
2359
2360 // Back solve first level of static condensation for interior
2361 // velocity space and store in F_int
2362 F_int = F_int + Transpose(*m_mat[mode].m_D_int) * F_p -
2363 Transpose(*m_mat[mode].m_Btilde) * F_bnd;
2364 F_int = (*m_mat[mode].m_Cinv) * F_int;
2365
2366 // Unpack solution from Bnd and F_int to v_coeffs
2367 cnt = cnt1 = 0;
2368 for (i = 0; i < nel; ++i) // loop over elements
2369 {
2370 fields[0]->GetExp(i)->GetBoundaryMap(bmap);
2371 fields[0]->GetExp(i)->GetInteriorMap(imap);
2372 nbnd = bmap.size();
2373 nint = imap.size();
2374 offset = fields[0]->GetCoeff_Offset(i);
2375
2376 for (j = 0; j < nvel; ++j) // loop over velocity fields
2377 {
2378 for (n = 0; n < nz_loc; ++n)
2379 {
2380 for (k = 0; k < nbnd; ++k)
2381 {
2382 fields[j]->SetCoeff(n * nplanecoeffs + offset + bmap[k],
2383 f_bnd[cnt + k]);
2384 }
2385
2386 for (k = 0; k < nint; ++k)
2387 {
2388 fields[j]->SetCoeff(n * nplanecoeffs + offset + imap[k],
2389 f_int[cnt1 + k]);
2390 }
2391 cnt += nbnd;
2392 cnt1 += nint;
2393 }
2394 }
2395 }
2396
2397 for (j = 0; j < nvel; ++j)
2398 {
2399 fields[j]->SetPhysState(false);
2400 }
2401}
2402
2403void CoupledLinearNS::v_Output(void)
2404{
2405 std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size() + 1);
2406 std::vector<std::string> variables(m_fields.size() + 1);
2407 size_t i;
2408
2409 for (i = 0; i < m_fields.size(); ++i)
2410 {
2411 fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2412 variables[i] = m_boundaryConditions->GetVariable(i);
2413 }
2414
2415 fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2416 // project pressure field to velocity space
2417 if (m_singleMode == true)
2418 {
2419 Array<OneD, NekDouble> tmpfieldcoeffs(m_fields[0]->GetNcoeffs() / 2);
2420 m_pressure->GetPlane(0)->BwdTrans(
2421 m_pressure->GetPlane(0)->GetCoeffs(),
2422 m_pressure->GetPlane(0)->UpdatePhys());
2423 m_pressure->GetPlane(1)->BwdTrans(
2424 m_pressure->GetPlane(1)->GetCoeffs(),
2425 m_pressure->GetPlane(1)->UpdatePhys());
2426 m_fields[0]->GetPlane(0)->FwdTransLocalElmt(
2427 m_pressure->GetPlane(0)->GetPhys(), fieldcoeffs[i]);
2428 m_fields[0]->GetPlane(1)->FwdTransLocalElmt(
2429 m_pressure->GetPlane(1)->GetPhys(), tmpfieldcoeffs);
2430 for (int e = 0; e < m_fields[0]->GetNcoeffs() / 2; e++)
2431 {
2432 fieldcoeffs[i][e + m_fields[0]->GetNcoeffs() / 2] =
2433 tmpfieldcoeffs[e];
2434 }
2435 }
2436 else
2437 {
2438 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
2439 m_fields[0]->FwdTransLocalElmt(m_pressure->GetPhys(), fieldcoeffs[i]);
2440 }
2441 variables[i] = "p";
2442
2443 if (!m_comm->IsParallelInTime())
2444 {
2445 WriteFld(m_sessionName + ".fld", m_fields[0], fieldcoeffs, variables);
2446 }
2447 else
2448 {
2449 std::string newdir = m_sessionName + ".pit";
2450 if (!fs::is_directory(newdir))
2451 {
2452 fs::create_directory(newdir);
2453 }
2454 WriteFld(
2455 newdir + "/" + m_sessionName + "_" +
2456 std::to_string(m_windowPIT * m_comm->GetTimeComm()->GetSize() +
2457 m_comm->GetTimeComm()->GetRank() + 1) +
2458 ".fld",
2459 m_fields[0], fieldcoeffs, variables);
2460 }
2461}
2462
2463int CoupledLinearNS::v_GetForceDimension()
2464{
2465 return m_session->GetVariables().size();
2466}
2467} // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition ErrorUtil.hpp:250
sign
#define sign(a, b)
return the sign(b)*a
Definition Polylib.cpp:47
Nektar::CoupledLinearNS::L2Norm
void L2Norm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition CoupledLinearNS.cpp:1896
Nektar::CoupledLinearNS::CoupledLinearNS
CoupledLinearNS(const LibUtilities::SessionReaderSharedPtr &pSesssion, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition CoupledLinearNS.cpp:61
Nektar::CoupledLinearNS::v_GenerateSummary
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Virtual function for generating summary information.
Definition CoupledLinearNS.cpp:1356
Nektar::CoupledLinearNS::EvaluateNewtonRHS
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition CoupledLinearNS.cpp:1911
Nektar::CoupledLinearNS::v_DoSolve
void v_DoSolve(void) override
Virtual function for solve implementation.
Definition CoupledLinearNS.cpp:1568
Nektar::CoupledLinearNS::m_locToGloMap
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
Definition CoupledLinearNS.h:162
Nektar::CoupledLinearNS::SolveUnsteadyStokesSystem
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
Definition CoupledLinearNS.cpp:1506
Nektar::CoupledLinearNS::SolveLinearNS
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Definition CoupledLinearNS.cpp:2069
Nektar::CoupledLinearNS::v_Output
void v_Output(void) override
Definition CoupledLinearNS.cpp:2403
Nektar::CoupledLinearNS::v_DoInitialise
void v_DoInitialise(bool dumpInitialConditions=true) override
Virtual function for initialisation implementation.
Definition CoupledLinearNS.cpp:1361
Nektar::CoupledLinearNS::solverTypeLookupId
static std::string solverTypeLookupId
Definition CoupledLinearNS.h:170
Nektar::CoupledLinearNS::GenPressureExp
const SpatialDomains::ExpansionInfoMap & GenPressureExp(const SpatialDomains::ExpansionInfoMap &VelExp)
Definition CoupledLinearNS.cpp:1967
Nektar::CoupledLinearNS::v_NegatedOp
bool v_NegatedOp(void) override
Definition CoupledLinearNS.cpp:1647
Nektar::CoupledLinearNS::EvaluateAdvection
void EvaluateAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition CoupledLinearNS.cpp:1493
Nektar::CoupledLinearNS::create
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition CoupledLinearNS.h:92
Nektar::CoupledLinearNS::v_InitObject
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
Definition CoupledLinearNS.cpp:69
Nektar::CoupledLinearNS::m_zeroMode
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
Definition CoupledLinearNS.h:174
Nektar::CoupledLinearNS::m_ForcingTerm_Coeffs
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
Definition CoupledLinearNS.h:160
Nektar::CoupledLinearNS::m_ForcingTerm
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
Definition CoupledLinearNS.h:159
Nektar::CoupledLinearNS::v_TransPhysToCoeff
void v_TransPhysToCoeff(void) override
Virtual function for transformation to coefficient space.
Definition CoupledLinearNS.cpp:1557
Nektar::CoupledLinearNS::SetUpCoupledMatrix
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
Definition CoupledLinearNS.cpp:200
Nektar::CoupledLinearNS::className
static std::string className
Name of class.
Definition CoupledLinearNS.h:102
Nektar::CoupledLinearNS::v_GetForceDimension
int v_GetForceDimension(void) override
Definition CoupledLinearNS.cpp:2463
Nektar::CoupledLinearNS::InfNorm
void InfNorm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition CoupledLinearNS.cpp:1879
Nektar::CoupledLinearNS::v_TransCoeffToPhys
void v_TransCoeffToPhys(void) override
Virtual function for transformation to physical space.
Definition CoupledLinearNS.cpp:1546
Nektar::CoupledLinearNS::m_mat
Array< OneD, CoupledSolverMatrices > m_mat
Definition CoupledLinearNS.h:187
Nektar::IncNavierStokes
This class is the base class for Navier Stokes problems.
Definition IncNavierStokes.h:148
Nektar::IncNavierStokes::m_pressure
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition IncNavierStokes.h:181
Nektar::IncNavierStokes::m_kinvis
NekDouble m_kinvis
Kinematic viscosity.
Definition IncNavierStokes.h:183
Nektar::IncNavierStokes::v_InitObject
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
Definition IncNavierStokes.cpp:89
Nektar::IncNavierStokes::m_extrapolation
ExtrapolateSharedPtr m_extrapolation
Definition IncNavierStokes.h:161
Nektar::IncNavierStokes::m_velocity
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition IncNavierStokes.h:178
Nektar::IncNavierStokes::m_equationType
EquationType m_equationType
equation type;
Definition IncNavierStokes.h:188
Nektar::IncNavierStokes::m_nConvectiveFields
int m_nConvectiveFields
Number of fields to be convected;.
Definition IncNavierStokes.h:174
Nektar::IncNavierStokes::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Definition IncNavierStokes.h:171
Nektar::IncNavierStokes::EvaluateAdvectionTerms
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition IncNavierStokes.cpp:298
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition Basis.h:45
Nektar::LibUtilities::BasisKey::GetBasisType
BasisType GetBasisType() const
Return type of expansion basis.
Definition Basis.h:131
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition Basis.h:137
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition Basis.h:74
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition BasicUtils/NekFactory.hpp:196
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition BasicUtils/NekFactory.hpp:142
Nektar::LibUtilities::SessionReader::RegisterEnumValue
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
Definition SessionReader.h:576
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineImplicitSolve
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
Definition TimeIntegrationSchemeOperators.h:100
Nektar::LibUtilities::Timer::TimePerTest
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition Timer.cpp:65
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition NekMemoryManager.hpp:81
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition NekMemoryManager.hpp:166
Nektar::SolverUtils::AdvectionSystem::m_advObject
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
Definition AdvectionSystem.h:67
Nektar::SolverUtils::EquationSystem::m_useFFT
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
Definition EquationSystem.h:514
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition EquationSystem.h:471
Nektar::SolverUtils::EquationSystem::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition EquationSystem.h:457
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition EquationSystem.h:475
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition EquationSystem.h:467
Nektar::SolverUtils::EquationSystem::m_windowPIT
int m_windowPIT
Index of windows for parallel-in-time time iteration.
Definition EquationSystem.h:502
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition EquationSystem.h:351
Nektar::SolverUtils::EquationSystem::m_lambda
NekDouble m_lambda
Lambda constant in real system if one required.
Definition EquationSystem.h:484
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition EquationSystem.h:283
Nektar::SolverUtils::EquationSystem::m_npointsZ
int m_npointsZ
number of points in Z direction (if homogeneous)
Definition EquationSystem.h:564
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition EquationSystem.cpp:1260
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition EquationSystem.cpp:1322
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition EquationSystem.h:473
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::m_singleMode
bool m_singleMode
Flag to determine if single homogeneous mode is used.
Definition EquationSystem.h:508
Nektar::SolverUtils::EquationSystem::SetInitialConditions
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
Definition EquationSystem.h:171
Nektar::SolverUtils::EquationSystem::m_homogen_dealiasing
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
Definition EquationSystem.h:519
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition EquationSystem.h:341
Nektar::SolverUtils::EquationSystem::m_boundaryConditions
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Definition EquationSystem.h:469
Nektar::SolverUtils::EquationSystem::m_LhomZ
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition EquationSystem.h:560
Nektar::SolverUtils::EquationSystem::GetFunction
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Definition EquationSystem.cpp:748
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition UnsteadySystem.h:47
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition UnsteadySystem.h:89
Nektar::SolverUtils::UnsteadySystem::v_DoSolve
SOLVER_UTILS_EXPORT void v_DoSolve() override
Solves an unsteady problem.
Definition UnsteadySystem.cpp:205
Blas::Dcopy
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
Blas::Dgemm
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
Nektar::LibUtilities::BasisKeyVector
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
Definition FoundationsFwd.hpp:68
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition SessionReader.h:119
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition PointsType.h:59
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition ExpList.h:87
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition LocTraceToTraceMap.h:56
Nektar::MultiRegions::ContFieldSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition ContField.h:278
Nektar::NekConstants::kNekUnsetDouble
static const NekDouble kNekUnsetDouble
Definition NektarUnivConsts.hpp:44
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition Misc.h:46
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition EquationSystem.cpp:99
Nektar::SolverUtils::AddSummaryItem
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
Nektar::SpatialDomains::ExpansionInfoMapShPtr
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
Definition MeshGraph.h:186
Nektar::SpatialDomains::ExpansionInfoShPtr
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition MeshGraph.h:183
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition MeshGraph.h:217
Nektar::SpatialDomains::ExpansionInfoMap
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition MeshGraph.h:184
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition StdExpansion.h:1783
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition StdRegions.hpp:430
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition NekMatrixFwd.hpp:66
Nektar::eSteadyNavierStokes
@ eSteadyNavierStokes
Definition IncNavierStokes.h:76
Nektar::eUnsteadyNavierStokes
@ eUnsteadyNavierStokes
Definition IncNavierStokes.h:75
Nektar::eSteadyLinearisedNS
@ eSteadyLinearisedNS
Definition IncNavierStokes.h:72
Nektar::eCoupledLinearisedNS
@ eCoupledLinearisedNS
Definition IncNavierStokes.h:57
Nektar::NullNekDoubleArrayOfArray
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
Definition BasicUtils/SharedArray.hpp:929
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition NekTypeDefs.hpp:79
Nektar::GetExtrapolateFactory
ExtrapolateFactory & GetExtrapolateFactory()
Definition Extrapolate.cpp:47
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition BlockMatrix.hpp:235
Nektar::CoupledLocalToGlobalC0ContMapSharedPtr
std::shared_ptr< CoupledLocalToGlobalC0ContMap > CoupledLocalToGlobalC0ContMapSharedPtr
Definition CoupledLocalToGlobalC0ContMap.h:66
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition NekTypeDefs.hpp:75
Vmath::Vmul
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
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition Vmath.hpp:292
Vmath::Vadd
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
Vmath::Smul
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
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition Vmath.hpp:273
Vmath::Vsub
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
Nektar::coupledSolverMatrices::m_D_int
DNekScalBlkMatSharedPtr m_D_int
Inner product of pressure system with divergence of the interior velocity space .
Definition CoupledLinearNS.h:81
Nektar::coupledSolverMatrices::m_D_bnd
DNekScalBlkMatSharedPtr m_D_bnd
Inner product of pressure system with divergence of the boundary velocity space .
Definition CoupledLinearNS.h:75
Nektar::coupledSolverMatrices::m_Cinv
DNekScalBlkMatSharedPtr m_Cinv
Interior-Interior Laplaican plus linearised convective terms inverted, i.e. the inverse of .
Definition CoupledLinearNS.h:69
Nektar::coupledSolverMatrices::m_Btilde
DNekScalBlkMatSharedPtr m_Btilde
Interior-boundary Laplacian plus linearised convective terms .
Definition CoupledLinearNS.h:62
Nektar::coupledSolverMatrices::m_BCinv
DNekScalBlkMatSharedPtr m_BCinv
Boundary-interior Laplacian plus linearised convective terms pre-multiplying Cinv: .
Definition CoupledLinearNS.h:56
Nektar::coupledSolverMatrices::m_CoupledBndSys
MultiRegions::GlobalLinSysSharedPtr m_CoupledBndSys
Definition CoupledLinearNS.h:83
Generated on Sun Nov 9 2025 20:17:01 for Nektar++ by doxygen 1.9.8