Nektar++
CoupledLinearNS.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: CoupledLinearNS.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Coupled Solver for the Linearised Incompressible
32// Navier Stokes equations
33///////////////////////////////////////////////////////////////////////////////
34
35#include <boost/algorithm/string.hpp>
36
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 // LibUtilities::TimeIntegrationMethod intMethod;
1369 // std::string TimeIntStr =
1370 // m_session->GetSolverInfo("TIMEINTEGRATIONMETHOD");
1371 // int i;
1372 // for(i = 0; i < (int)
1373 // LibUtilities::SIZE_TimeIntegrationMethod; ++i)
1374 // {
1375 // if(boost::iequals(LibUtilities::TimeIntegrationMethodMap[i],TimeIntStr))
1376 // {
1377 // intMethod =
1378 // (LibUtilities::TimeIntegrationMethod)i;
1379 // break;
1380 // }
1381 // }
1382 //
1383 // ASSERTL0(i != (int)
1384 // LibUtilities::SIZE_TimeIntegrationMethod, "Invalid
1385 // time integration type.");
1386 //
1387 // m_integrationScheme =
1388 // LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(LibUtilities::TimeIntegrationMethodMap[intMethod]);
1389
1390 // Could defind this from IncNavierStokes class?
1391 m_ode.DefineOdeRhs(&CoupledLinearNS::EvaluateAdvection, this);
1392
1393 m_ode.DefineImplicitSolve(
1394 &CoupledLinearNS::SolveUnsteadyStokesSystem, this);
1395
1396 // Set initial condition using time t=0
1397
1398 SetInitialConditions(0.0, dumpInitialConditions);
1399 break;
1400 }
1401 case eSteadyStokes:
1402 SetUpCoupledMatrix(0.0);
1403 break;
1404 case eSteadyOseen:
1405 {
1406 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1407 for (size_t i = 0; i < m_velocity.size(); ++i)
1408 {
1409 AdvField[i] = Array<OneD, NekDouble>(
1410 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1411 }
1412
1413 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1414 "Advection Velocity section must be defined in "
1415 "session file.");
1416
1417 std::vector<std::string> fieldStr;
1418 for (size_t i = 0; i < m_velocity.size(); ++i)
1419 {
1420 fieldStr.push_back(
1421 m_boundaryConditions->GetVariable(m_velocity[i]));
1422 }
1423 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1424
1425 SetUpCoupledMatrix(0.0, AdvField, false);
1426 }
1427 break;
1428 case eSteadyNavierStokes:
1429 {
1430 m_session->LoadParameter("KinvisMin", m_kinvisMin);
1431 m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1432 m_session->LoadParameter("Tolerence", m_tol);
1433 m_session->LoadParameter("MaxIteration", m_maxIt);
1434 m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1435 m_session->LoadParameter("Restart", m_Restart);
1436
1437 DefineForcingTerm();
1438
1439 if (m_Restart == 1)
1440 {
1441 ASSERTL0(m_session->DefinesFunction("Restart"),
1442 "Restart section must be defined in session file.");
1443
1444 Array<OneD, Array<OneD, NekDouble>> Restart(m_velocity.size());
1445 for (size_t i = 0; i < m_velocity.size(); ++i)
1446 {
1447 Restart[i] = Array<OneD, NekDouble>(
1448 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1449 }
1450 std::vector<std::string> fieldStr;
1451 for (size_t i = 0; i < m_velocity.size(); ++i)
1452 {
1453 fieldStr.push_back(
1454 m_boundaryConditions->GetVariable(m_velocity[i]));
1455 }
1456 GetFunction("Restart")->Evaluate(fieldStr, Restart);
1457
1458 for (size_t i = 0; i < m_velocity.size(); ++i)
1459 {
1460 m_fields[m_velocity[i]]->FwdTransLocalElmt(
1461 Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1462 }
1463 std::cout << "Saving the RESTART file for m_kinvis = "
1464 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1465 << std::endl;
1466 }
1467 else // We solve the Stokes Problem
1468 {
1469
1470 SetUpCoupledMatrix(0.0);
1471 m_initialStep = true;
1472 m_counter = 1;
1473 // SolveLinearNS(m_ForcingTerm_Coeffs);
1474 Solve();
1475 m_initialStep = false;
1476 std::cout << "Saving the Stokes Flow for m_kinvis = "
1477 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1478 << std::endl;
1479 }
1480 }
1481 break;
1482 case eSteadyLinearisedNS:
1483 {
1484 SetInitialConditions(0.0, dumpInitialConditions);
1485
1486 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1487 for (size_t i = 0; i < m_velocity.size(); ++i)
1488 {
1489 AdvField[i] = Array<OneD, NekDouble>(
1490 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1491 }
1492
1493 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1494 "Advection Velocity section must be defined in "
1495 "session file.");
1496
1497 std::vector<std::string> fieldStr;
1498 for (size_t i = 0; i < m_velocity.size(); ++i)
1499 {
1500 fieldStr.push_back(
1501 m_boundaryConditions->GetVariable(m_velocity[i]));
1502 }
1503 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1504
1505 SetUpCoupledMatrix(m_lambda, AdvField, true);
1506 }
1507 break;
1508 case eNoEquationType:
1509 default:
1510 ASSERTL0(false,
1511 "Unknown or undefined equation type for CoupledLinearNS");
1512 }
1513}
1514
1515void CoupledLinearNS::EvaluateAdvection(
1516 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1517 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
1518{
1519 // evaluate convection terms
1520 EvaluateAdvectionTerms(inarray, outarray, time);
1521
1522 for (auto &x : m_forcing)
1523 {
1524 x->Apply(m_fields, outarray, outarray, time);
1525 }
1526}
1527
1528void CoupledLinearNS::SolveUnsteadyStokesSystem(
1529 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1530 Array<OneD, Array<OneD, NekDouble>> &outarray,
1531 [[maybe_unused]] const NekDouble time, const NekDouble aii_Dt)
1532{
1533 size_t i;
1534 Array<OneD, Array<OneD, NekDouble>> F(m_nConvectiveFields);
1535 NekDouble lambda = 1.0 / aii_Dt;
1536 static NekDouble lambda_store;
1537 Array<OneD, Array<OneD, NekDouble>> forcing(m_velocity.size());
1538 // Matrix solution
1539 if (fabs(lambda_store - lambda) > 1e-10)
1540 {
1541 SetUpCoupledMatrix(lambda);
1542 lambda_store = lambda;
1543 }
1544
1545 // Forcing for advection solve
1546 for (i = 0; i < m_velocity.size(); ++i)
1547 {
1548 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1549 m_fields[m_velocity[i]]->SetWaveSpace(true);
1550 m_fields[m_velocity[i]]->IProductWRTBase(
1551 inarray[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1552 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1553 Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(), lambda,
1554 m_fields[m_velocity[i]]->GetCoeffs(), 1,
1555 m_fields[m_velocity[i]]->UpdateCoeffs(), 1);
1556 forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1557 }
1558
1559 SolveLinearNS(forcing);
1560
1561 for (i = 0; i < m_velocity.size(); ++i)
1562 {
1563 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1564 outarray[i]);
1565 }
1566}
1567
1568void CoupledLinearNS::v_TransCoeffToPhys(void)
1569{
1570 size_t nfields = m_fields.size();
1571 for (size_t k = 0; k < nfields; ++k)
1572 {
1573 // Backward Transformation in physical space for time evolution
1574 m_fields[k]->BwdTrans(m_fields[k]->GetCoeffs(),
1575 m_fields[k]->UpdatePhys());
1576 }
1577}
1578
1579void CoupledLinearNS::v_TransPhysToCoeff(void)
1580{
1581 size_t nfields = m_fields.size();
1582 for (size_t k = 0; k < nfields; ++k)
1583 {
1584 // Forward Transformation in physical space for time evolution
1585 m_fields[k]->FwdTransLocalElmt(m_fields[k]->GetPhys(),
1586 m_fields[k]->UpdateCoeffs());
1587 }
1588}
1589
1590void CoupledLinearNS::v_DoSolve(void)
1591{
1592 switch (m_equationType)
1593 {
1594 case eUnsteadyStokes:
1595 case eUnsteadyNavierStokes:
1596 // AdvanceInTime(m_steps);
1597 UnsteadySystem::v_DoSolve();
1598 break;
1599 case eSteadyStokes:
1600 case eSteadyOseen:
1601 case eSteadyLinearisedNS:
1602 {
1603 Solve();
1604 break;
1605 }
1606 case eSteadyNavierStokes:
1607 {
1608 LibUtilities::Timer Generaltimer;
1609 Generaltimer.Start();
1610
1611 int Check(0);
1612
1613 // Saving the init datas
1614 Checkpoint_Output(Check);
1615 Check++;
1616
1617 std::cout
1618 << "We execute INITIALLY SolveSteadyNavierStokes for m_kinvis "
1619 "= "
1620 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1621 << std::endl;
1622 SolveSteadyNavierStokes();
1623
1624 while (m_kinvis > m_kinvisMin)
1625 {
1626 if (Check == 1)
1627 {
1628 std::cout
1629 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1630 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1631 << std::endl;
1632 SolveSteadyNavierStokes();
1633 Checkpoint_Output(Check);
1634 Check++;
1635 }
1636
1637 Continuation();
1638
1639 if (m_kinvis > m_kinvisMin)
1640 {
1641 std::cout
1642 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1643 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1644 << std::endl;
1645 SolveSteadyNavierStokes();
1646 Checkpoint_Output(Check);
1647 Check++;
1648 }
1649 }
1650
1651 Generaltimer.Stop();
1652 std::cout << "\nThe total calculation time is : "
1653 << Generaltimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1654
1655 break;
1656 }
1657 case eNoEquationType:
1658 default:
1659 ASSERTL0(false,
1660 "Unknown or undefined equation type for CoupledLinearNS");
1661 }
1662}
1663
1664/** Virtual function to define if operator in DoSolve is negated
1665 * with regard to the strong form. This is currently only used in
1666 * Arnoldi solves. For Coupledd solver this is true since Stokes
1667 * operator is set up as a LHS rather than RHS operation
1668 */
1669bool CoupledLinearNS::v_NegatedOp(void)
1670{
1671 return true;
1672}
1673
1674void CoupledLinearNS::Solve(void)
1675{
1676 const size_t ncmpt = m_velocity.size();
1677 Array<OneD, Array<OneD, NekDouble>> forcing_phys(ncmpt);
1678 Array<OneD, Array<OneD, NekDouble>> forcing(ncmpt);
1679
1680 for (size_t i = 0; i < ncmpt; ++i)
1681 {
1682 forcing_phys[i] =
1683 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1684 forcing[i] =
1685 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNcoeffs(), 0.0);
1686 }
1687
1688 for (auto &x : m_forcing)
1689 {
1690 const NekDouble time = 0;
1691 x->Apply(m_fields, forcing_phys, forcing_phys, time);
1692 }
1693 for (size_t i = 0; i < ncmpt; ++i)
1694 {
1695 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1696 m_fields[i]->SetWaveSpace(true);
1697 m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1698 m_fields[i]->SetWaveSpace(waveSpace);
1699 }
1700
1701 SolveLinearNS(forcing);
1702}
1703
1704void CoupledLinearNS::DefineForcingTerm(void)
1705{
1706 m_ForcingTerm = Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1707 m_ForcingTerm_Coeffs =
1708 Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1709
1710 for (size_t i = 0; i < m_velocity.size(); ++i)
1711 {
1712 m_ForcingTerm[i] = Array<OneD, NekDouble>(
1713 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1714 m_ForcingTerm_Coeffs[i] =
1715 Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1716 }
1717
1718 if (m_session->DefinesFunction("ForcingTerm"))
1719 {
1720 std::vector<std::string> fieldStr;
1721 for (size_t i = 0; i < m_velocity.size(); ++i)
1722 {
1723 fieldStr.push_back(
1724 m_boundaryConditions->GetVariable(m_velocity[i]));
1725 }
1726 GetFunction("ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1727 for (size_t i = 0; i < m_velocity.size(); ++i)
1728 {
1729 m_fields[m_velocity[i]]->FwdTransLocalElmt(m_ForcingTerm[i],
1730 m_ForcingTerm_Coeffs[i]);
1731 }
1732 }
1733 else
1734 {
1735 std::cout
1736 << "'ForcingTerm' section has not been defined in the input file "
1737 "=> forcing=0"
1738 << std::endl;
1739 }
1740}
1741
1742void CoupledLinearNS::SolveSteadyNavierStokes(void)
1743{
1744 LibUtilities::Timer Newtontimer;
1745 Newtontimer.Start();
1746
1747 Array<OneD, Array<OneD, NekDouble>> RHS_Coeffs(m_velocity.size());
1748 Array<OneD, Array<OneD, NekDouble>> RHS_Phys(m_velocity.size());
1749 Array<OneD, Array<OneD, NekDouble>> delta_velocity_Phys(m_velocity.size());
1750 Array<OneD, Array<OneD, NekDouble>> Velocity_Phys(m_velocity.size());
1751 Array<OneD, NekDouble> L2_norm(m_velocity.size(), 1.0);
1752 Array<OneD, NekDouble> Inf_norm(m_velocity.size(), 1.0);
1753
1754 for (size_t i = 0; i < m_velocity.size(); ++i)
1755 {
1756 delta_velocity_Phys[i] = Array<OneD, NekDouble>(
1757 m_fields[m_velocity[i]]->GetTotPoints(), 1.0);
1758 Velocity_Phys[i] = Array<OneD, NekDouble>(
1759 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1760 }
1761
1762 m_counter = 1;
1763
1764 L2Norm(delta_velocity_Phys, L2_norm);
1765
1766 // while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1767 while (std::max(L2_norm[0], L2_norm[1]) > m_tol)
1768 {
1769 if (m_counter == 1)
1770 // At the first Newton step, we use the solution of the
1771 // Stokes problem (if Restart=0 in input file) Or the
1772 // solution of the .rst file (if Restart=1 in input
1773 // file)
1774 {
1775 for (size_t i = 0; i < m_velocity.size(); ++i)
1776 {
1777 RHS_Coeffs[i] = Array<OneD, NekDouble>(
1778 m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1779 RHS_Phys[i] = Array<OneD, NekDouble>(
1780 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1781 }
1782
1783 for (size_t i = 0; i < m_velocity.size(); ++i)
1784 {
1785 m_fields[m_velocity[i]]->BwdTrans(
1786 m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1787 }
1788
1789 m_initialStep = true;
1790 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1791 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1792 SolveLinearNS(RHS_Coeffs);
1793 m_initialStep = false;
1794 }
1795 if (m_counter > 1)
1796 {
1797 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1798
1799 if (m_counter % m_MatrixSetUpStep ==
1800 0) // Setting Up the matrix is expensive. We do it at each
1801 // "m_MatrixSetUpStep" step.
1802 {
1803 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1804 }
1805 SolveLinearNS(RHS_Coeffs);
1806 }
1807
1808 for (size_t i = 0; i < m_velocity.size(); ++i)
1809 {
1810 m_fields[m_velocity[i]]->BwdTrans(RHS_Coeffs[i], RHS_Phys[i]);
1811 m_fields[m_velocity[i]]->BwdTrans(
1812 m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1813 }
1814
1815 for (size_t i = 0; i < m_velocity.size(); ++i)
1816 {
1817 Vmath::Vadd(Velocity_Phys[i].size(), Velocity_Phys[i], 1,
1818 delta_velocity_Phys[i], 1, Velocity_Phys[i], 1);
1819 }
1820
1821 // InfNorm(delta_velocity_Phys, Inf_norm);
1822 L2Norm(delta_velocity_Phys, L2_norm);
1823
1824 if (std::max(Inf_norm[0], Inf_norm[1]) > 100)
1825 {
1826 std::cout << "\nThe Newton method has failed at m_kinvis = "
1827 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1828 << std::endl;
1829 ASSERTL0(0, "The Newton method has failed... \n");
1830 }
1831
1832 std::cout << "\n";
1833 m_counter++;
1834 }
1835
1836 if (m_counter > 1) // We save u:=u+\delta u in u->Coeffs
1837 {
1838 for (size_t i = 0; i < m_velocity.size(); ++i)
1839 {
1840 m_fields[m_velocity[i]]->FwdTrans(
1841 Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1842 }
1843 }
1844
1845 Newtontimer.Stop();
1846 std::cout << "We have done " << m_counter - 1 << " iteration(s) in "
1847 << Newtontimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1848}
1849
1850void CoupledLinearNS::Continuation(void)
1851{
1852 Array<OneD, Array<OneD, NekDouble>> u_N(m_velocity.size());
1853 Array<OneD, Array<OneD, NekDouble>> tmp_RHS(m_velocity.size());
1854 Array<OneD, Array<OneD, NekDouble>> RHS(m_velocity.size());
1855 Array<OneD, Array<OneD, NekDouble>> u_star(m_velocity.size());
1856
1857 std::cout << "We apply the continuation method: " << std::endl;
1858
1859 for (size_t i = 0; i < m_velocity.size(); ++i)
1860 {
1861 u_N[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1862 0.0);
1863 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1864 u_N[i]);
1865
1866 RHS[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1867 0.0);
1868 tmp_RHS[i] = Array<OneD, NekDouble>(
1869 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1870
1871 m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1872 Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1873
1874 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1875 m_fields[m_velocity[i]]->SetWaveSpace(true);
1876 m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1877 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1878 }
1879
1880 SetUpCoupledMatrix(0.0, u_N, true);
1881 SolveLinearNS(RHS);
1882
1883 for (size_t i = 0; i < m_velocity.size(); ++i)
1884 {
1885 u_star[i] = Array<OneD, NekDouble>(
1886 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1887 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1888 u_star[i]);
1889
1890 // u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1891 Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1892 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1893
1894 m_fields[m_velocity[i]]->FwdTrans(
1895 u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1896 }
1897
1898 m_kinvis -= m_kinvis * m_KinvisPercentage / 100;
1899}
1900
1901void CoupledLinearNS::InfNorm(Array<OneD, Array<OneD, NekDouble>> &inarray,
1902 Array<OneD, NekDouble> &outarray)
1903{
1904 for (size_t i = 0; i < m_velocity.size(); ++i)
1905 {
1906 outarray[i] = 0.0;
1907 for (size_t j = 0; j < inarray[i].size(); ++j)
1908 {
1909 if (inarray[i][j] > outarray[i])
1910 {
1911 outarray[i] = inarray[i][j];
1912 }
1913 }
1914 std::cout << "InfNorm[" << i << "] = " << outarray[i] << std::endl;
1915 }
1916}
1917
1918void CoupledLinearNS::L2Norm(Array<OneD, Array<OneD, NekDouble>> &inarray,
1919 Array<OneD, NekDouble> &outarray)
1920{
1921 for (size_t i = 0; i < m_velocity.size(); ++i)
1922 {
1923 outarray[i] = 0.0;
1924 for (size_t j = 0; j < inarray[i].size(); ++j)
1925 {
1926 outarray[i] += inarray[i][j] * inarray[i][j];
1927 }
1928 outarray[i] = std::sqrt(outarray[i]);
1929 std::cout << "L2Norm[" << i << "] = " << outarray[i] << std::endl;
1930 }
1931}
1932
1933void CoupledLinearNS::EvaluateNewtonRHS(
1934 Array<OneD, Array<OneD, NekDouble>> &Velocity,
1935 Array<OneD, Array<OneD, NekDouble>> &outarray)
1936{
1937 Array<OneD, Array<OneD, NekDouble>> Eval_Adv(m_velocity.size());
1938 Array<OneD, Array<OneD, NekDouble>> tmp_DerVel(m_velocity.size());
1939 Array<OneD, Array<OneD, NekDouble>> AdvTerm(m_velocity.size());
1940 Array<OneD, Array<OneD, NekDouble>> ViscTerm(m_velocity.size());
1941 Array<OneD, Array<OneD, NekDouble>> Forc(m_velocity.size());
1942
1943 for (size_t i = 0; i < m_velocity.size(); ++i)
1944 {
1945 Eval_Adv[i] = Array<OneD, NekDouble>(
1946 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1947 tmp_DerVel[i] = Array<OneD, NekDouble>(
1948 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1949
1950 AdvTerm[i] = Array<OneD, NekDouble>(
1951 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1952 ViscTerm[i] = Array<OneD, NekDouble>(
1953 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1954 Forc[i] = Array<OneD, NekDouble>(
1955 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1956 outarray[i] = Array<OneD, NekDouble>(
1957 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1958
1959 m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1960
1961 Vmath::Smul(tmp_DerVel[i].size(), m_kinvis, tmp_DerVel[i], 1,
1962 tmp_DerVel[i], 1);
1963 }
1964
1965 EvaluateAdvectionTerms(Velocity, Eval_Adv, m_time);
1966
1967 for (size_t i = 0; i < m_velocity.size(); ++i)
1968 {
1969 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1970 m_fields[m_velocity[i]]->SetWaveSpace(true);
1971 m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i],
1972 AdvTerm[i]); //(w, (u.grad)u)
1973 m_fields[m_velocity[i]]->IProductWRTDerivBase(
1974 i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1975 m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i],
1976 Forc[i]); //(w, f)
1977 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1978
1979 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1980 outarray[i], 1);
1981 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1982 outarray[i], 1);
1983
1984 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1985 1);
1986 }
1987}
1988
1989const SpatialDomains::ExpansionInfoMap &CoupledLinearNS::GenPressureExp(
1990 const SpatialDomains::ExpansionInfoMap &VelExp)
1991{
1992 SpatialDomains::ExpansionInfoMapShPtr returnval;
1993 returnval =
1994 MemoryManager<SpatialDomains::ExpansionInfoMap>::AllocateSharedPtr();
1995
1996 int nummodes;
1997
1998 for (auto &expMapIter : VelExp)
1999 {
2000 LibUtilities::BasisKeyVector BasisVec;
2001
2002 for (size_t i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
2003 {
2004 LibUtilities::BasisKey B = expMapIter.second->m_basisKeyVector[i];
2005 nummodes = B.GetNumModes();
2006 ASSERTL0(nummodes > 3,
2007 "Velocity polynomial space not sufficiently high (>= 4)");
2008 // Should probably set to be an orthogonal basis.
2009 LibUtilities::BasisKey newB(B.GetBasisType(), nummodes - 2,
2010 B.GetPointsKey());
2011 BasisVec.push_back(newB);
2012 }
2013
2014 // Put new expansion into list.
2015 SpatialDomains::ExpansionInfoShPtr expansionElementShPtr =
2016 MemoryManager<SpatialDomains::ExpansionInfo>::AllocateSharedPtr(
2017 expMapIter.second->m_geomShPtr, BasisVec);
2018 (*returnval)[expMapIter.first] = expansionElementShPtr;
2019 }
2020
2021 // Save expansion into graph.
2022 m_graph->SetExpansionInfo("p", returnval);
2023
2024 return *returnval;
2025}
2026
2027/**
2028 * @param forcing A list of forcing functions for each velocity
2029 * component
2030 *
2031 * The routine involves two levels of static
2032 * condensations. Initially we require a statically condensed
2033 * forcing function which requires the following manipulation
2034 *
2035 * \f[ {F\_bnd} = {\bf f}_{bnd} -m\_B \,m\_Cinv\, {\bf f}_{int},
2036 * \hspace{1cm} F\_p = m\_D\_{int}\, m\_Cinv\, {\bf f}_{int} \f]
2037 *
2038 * Where \f${\bf f}_{bnd}\f$ denote the forcing degrees of
2039 * freedom of the elemental velocities on the boundary of the
2040 * element, \f${\bf f}_{int}\f$ denote the forcing degrees of
2041 * freedom of the elemental velocities on the interior of the
2042 * element. (see detailed description for more details).
2043 *
2044 * This vector is further manipulated into
2045 *
2046 * \f[ Fh\_{bnd} = \left [ \begin{array}{c} f\_{bnd} -m\_B \,
2047 * m\_Cinv\, {\bf f}_{int}\\ \left [m\_D\_{int} \, m\_Cinv \,{\bf
2048 * f}_{int} \right]_0 \end{array}\right ] \hspace{1cm} [Fh\_p]_{i} =
2049 * \begin{array}{c} [m\_D\_{int} \, m\_Cinv \, {\bf
2050 * f}_{int}]_{i+1} \end{array} \f]
2051 *
2052 * where \f$-{[}m\_D\_{int}^T\, m\_Cinv \,{\bf f}_{int}]_0\f$
2053 * which is corresponds to the mean mode of the pressure degrees
2054 * of freedom and is now added to the boundary system and the
2055 * remainder of the block becomes the interior forcing for the
2056 * inner static condensation (see detailed description for more
2057 * details) which is set up in a GlobalLinSysDirectStaticCond
2058 * class.
2059 *
2060 * Finally we perform the final maniplation of the forcing to
2061 * using hte
2062 * \f[ Fh\_{bnd} = Fh\_{bnd} - m\_Bh \,m\_Chinv \, Fh\_p \f]
2063 *
2064 * We can now call the solver to the global coupled boundary
2065 * system (through the call to #m_CoupledBndSys->Solve) to obtain
2066 * the velocity boundary solution as the mean pressure solution,
2067 * i.e.
2068 *
2069 * \f[ {\cal A}^T(\hat{A} - \hat{C}^T \hat{D}^{-1} \hat{B} ){\cal
2070 * A} \, Bnd = Fh\_{bnd} \f]
2071 *
2072 * Once we know the solution to the above the rest of the pressure
2073 * modes are recoverable thorugh
2074 *
2075 * \f[ Ph = m\_Dhinv\, (Bnd - m\_Ch^T \, Fh_{bnd}) \f]
2076 *
2077 * We can now unpack \f$ Fh\_{bnd} \f$ (last elemental mode) and
2078 * \f$ Ph \f$ into #m_pressure and \f$ F_p\f$ and \f$ Fh\_{bnd}\f$
2079 * into a closed pack list of boundary velocoity degrees of
2080 * freedom stored in \f$ F\_bnd\f$.
2081 *
2082 * Finally the interior velocity degrees of freedom are then
2083 * obtained through the relationship
2084 *
2085 * \f[ F\_{int} = m\_Cinv\ ( F\_{int} + m\_D\_{int}^T\,
2086 * F\_p - m\_Btilde^T\, Bnd) \f]
2087 *
2088 * We then unpack the solution back to the MultiRegion structures
2089 * #m_velocity and #m_pressure
2090 */
2091void CoupledLinearNS::SolveLinearNS(
2092 const Array<OneD, Array<OneD, NekDouble>> &forcing)
2093{
2094 size_t i;
2095 Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.size());
2096 Array<OneD, Array<OneD, NekDouble>> force(m_velocity.size());
2097
2098 // Impose Dirichlet conditions on velocity fields
2099 for (i = 0; i < m_velocity.size(); ++i)
2100 {
2101 Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
2102 m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
2103 }
2104
2105 if (m_HomogeneousType == eHomogeneous1D)
2106 {
2107 int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
2108 for (int n = 0; n < m_npointsZ / 2; ++n)
2109 {
2110 // Get the a Fourier mode of velocity and forcing components.
2111 for (i = 0; i < m_velocity.size(); ++i)
2112 {
2113 vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2 * n);
2114 // Note this needs to correlate with how we pass forcing
2115 force[i] = forcing[i] + 2 * n * ncoeffsplane;
2116 }
2117
2118 SolveLinearNS(force, vel_fields, m_pressure->GetPlane(2 * n), n);
2119 }
2120 for (i = 0; i < m_velocity.size(); ++i)
2121 {
2122 m_fields[m_velocity[i]]->SetPhysState(false);
2123 }
2124 m_pressure->SetPhysState(false);
2125 }
2126 else
2127 {
2128 for (i = 0; i < m_velocity.size(); ++i)
2129 {
2130 vel_fields[i] = m_fields[m_velocity[i]];
2131 // Note this needs to correlate with how we pass forcing
2132 force[i] = forcing[i];
2133 }
2134 SolveLinearNS(force, vel_fields, m_pressure);
2135 }
2136}
2137
2138void CoupledLinearNS::SolveLinearNS(
2139 Array<OneD, Array<OneD, NekDouble>> &forcing,
2140 Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
2141 MultiRegions::ExpListSharedPtr &pressure, const int mode)
2142{
2143 size_t i, j, k, n, cnt, cnt1;
2144 size_t nbnd, nint, offset;
2145 size_t nvel = m_velocity.size();
2146 size_t nel = fields[0]->GetNumElmts();
2147 Array<OneD, unsigned int> bmap, imap;
2148
2149 Array<OneD, NekDouble> f_bnd(m_mat[mode].m_BCinv->GetRows());
2150 NekVector<NekDouble> F_bnd(f_bnd.size(), f_bnd, eWrapper);
2151 Array<OneD, NekDouble> f_int(m_mat[mode].m_BCinv->GetColumns());
2152 NekVector<NekDouble> F_int(f_int.size(), f_int, eWrapper);
2153
2154 size_t nz_loc;
2155 int nplanecoeffs =
2156 fields[m_velocity[0]]
2157 ->GetNcoeffs(); // this is fine since we pass the nplane coeff data.
2158
2159 if (mode) // Homogeneous mode flag
2160 {
2161 nz_loc = 2;
2162 }
2163 else
2164 {
2165 if (m_singleMode)
2166 {
2167 nz_loc = 2;
2168 }
2169 else
2170 {
2171 nz_loc = 1;
2172 if (m_HomogeneousType == eHomogeneous1D)
2173 {
2174 Array<OneD, NekDouble> tmp;
2175 // Zero fields to set complex mode to zero;
2176 for (i = 0; i < fields.size(); ++i)
2177 {
2178 Vmath::Zero(nplanecoeffs,
2179 tmp = fields[i]->UpdateCoeffs() + nplanecoeffs,
2180 1);
2181 }
2182 Vmath::Zero(2 * pressure->GetNcoeffs(),
2183 pressure->UpdateCoeffs(), 1);
2184 }
2185 }
2186 }
2187
2188 for (k = 0; k < nvel; ++k)
2189 {
2190 MultiRegions::ContFieldSharedPtr cfield =
2191 std::dynamic_pointer_cast<MultiRegions::ContField>(fields[k]);
2192
2193 Array<OneD, NekDouble> sign =
2194 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsSign();
2195 const Array<OneD, const int> map =
2196 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsMap();
2197
2198 // Add weak boundary conditions to forcing
2199 const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConds =
2200 fields[k]->GetBndConditions();
2201 Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp;
2202
2203 if (m_HomogeneousType == eHomogeneous1D)
2204 {
2205 bndCondExp =
2206 m_fields[k]->GetPlane(2 * mode)->GetBndCondExpansions();
2207 }
2208 else
2209 {
2210 bndCondExp = m_fields[k]->GetBndCondExpansions();
2211 }
2212
2213 for (n = 0; n < nz_loc; ++n)
2214 {
2215 int bndcnt = 0;
2216 for (i = 0; i < bndCondExp.size(); ++i)
2217 {
2218 const Array<OneD, const NekDouble> bndcoeffs =
2219 bndCondExp[i]->GetCoeffs();
2220 size_t nCoeffs = (bndCondExp[i])->GetNcoeffs();
2221
2222 cnt = 0;
2223 if (bndConds[i]->GetBoundaryConditionType() ==
2224 SpatialDomains::eNeumann ||
2225 bndConds[i]->GetBoundaryConditionType() ==
2226 SpatialDomains::eRobin)
2227 {
2228 if (m_locToGloMap[mode]->GetSignChange())
2229 {
2230 for (j = 0; j < nCoeffs; j++)
2231 {
2232 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2233 sign[bndcnt + j] * bndcoeffs[j];
2234 }
2235 }
2236 else
2237 {
2238 for (j = 0; j < nCoeffs; j++)
2239 {
2240 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2241 bndcoeffs[j];
2242 }
2243 }
2244 }
2245
2246 bndcnt += bndCondExp[i]->GetNcoeffs();
2247 }
2248 }
2249 }
2250
2251 Array<OneD, NekDouble> bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(), 0.0);
2252
2253 // Construct f_bnd and f_int and fill in bnd from inarray
2254 // (with Dirichlet BCs imposed)
2255 int bndoffset = 0;
2256 cnt = cnt1 = 0;
2257 for (i = 0; i < nel; ++i) // loop over elements
2258 {
2259 fields[m_velocity[0]]->GetExp(i)->GetBoundaryMap(bmap);
2260 fields[m_velocity[0]]->GetExp(i)->GetInteriorMap(imap);
2261 nbnd = bmap.size();
2262 nint = imap.size();
2263 offset = fields[m_velocity[0]]->GetCoeff_Offset(i);
2264
2265 for (j = 0; j < nvel; ++j) // loop over velocity fields
2266 {
2267 Array<OneD, NekDouble> incoeffs = fields[j]->UpdateCoeffs();
2268
2269 for (n = 0; n < nz_loc; ++n)
2270 {
2271 for (k = 0; k < nbnd; ++k)
2272 {
2273 f_bnd[cnt + k] =
2274 forcing[j][n * nplanecoeffs + offset + bmap[k]];
2275
2276 bnd[bndoffset + (n + j * nz_loc) * nbnd + k] =
2277 incoeffs[n * nplanecoeffs + offset + bmap[k]];
2278 }
2279 for (k = 0; k < nint; ++k)
2280 {
2281 f_int[cnt1 + k] =
2282 forcing[j][n * nplanecoeffs + offset + imap[k]];
2283 }
2284
2285 cnt += nbnd;
2286 cnt1 += nint;
2287 }
2288 }
2289 bndoffset +=
2290 nvel * nz_loc * nbnd + nz_loc * (pressure->GetExp(i)->GetNcoeffs());
2291 }
2292
2293 Array<OneD, NekDouble> f_p(m_mat[mode].m_D_int->GetRows());
2294 NekVector<NekDouble> F_p(f_p.size(), f_p, eWrapper);
2295 NekVector<NekDouble> F_p_tmp(m_mat[mode].m_Cinv->GetRows());
2296
2297 // fbnd does not currently hold the pressure mean
2298 F_bnd = F_bnd - (*m_mat[mode].m_BCinv) * F_int;
2299 F_p_tmp = (*m_mat[mode].m_Cinv) * F_int;
2300 F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
2301
2302 // construct inner forcing
2303 Array<OneD, NekDouble> fh_bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(),
2304 0.0);
2305
2306 offset = cnt = 0;
2307 for (i = 0; i < nel; ++i)
2308 {
2309 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2310
2311 for (j = 0; j < nvel; ++j)
2312 {
2313 for (k = 0; k < nbnd; ++k)
2314 {
2315 fh_bnd[offset + j * nbnd + k] = f_bnd[cnt + k];
2316 }
2317 cnt += nbnd;
2318 }
2319
2320 nint = pressure->GetExp(i)->GetNcoeffs();
2321 offset += nvel * nbnd + nint * nz_loc;
2322 }
2323
2324 offset = cnt1 = 0;
2325 for (i = 0; i < nel; ++i)
2326 {
2327 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2328 nint = pressure->GetExp(i)->GetNcoeffs();
2329
2330 for (n = 0; n < nz_loc; ++n)
2331 {
2332 for (j = 0; j < nint; ++j)
2333 {
2334 fh_bnd[offset + nvel * nbnd + n * nint + j] = f_p[cnt1 + j];
2335 }
2336 cnt1 += nint;
2337 }
2338 offset += nvel * nbnd + nz_loc * nint;
2339 }
2340 m_mat[mode].m_CoupledBndSys->Solve(fh_bnd, bnd, m_locToGloMap[mode]);
2341
2342 // unpack pressure and velocity boundary systems.
2343 offset = cnt = 0;
2344 int totpcoeffs = pressure->GetNcoeffs();
2345 Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2346 for (i = 0; i < nel; ++i)
2347 {
2348 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2349 nint = pressure->GetExp(i)->GetNcoeffs();
2350 for (j = 0; j < nvel; ++j)
2351 {
2352 for (k = 0; k < nbnd; ++k)
2353 {
2354 f_bnd[cnt + k] = bnd[offset + j * nbnd + k];
2355 }
2356 cnt += nbnd;
2357 }
2358 offset += nvel * nbnd + nint * nz_loc;
2359 }
2360
2361 pressure->SetPhysState(false);
2362
2363 offset = cnt = cnt1 = 0;
2364 for (i = 0; i < nel; ++i)
2365 {
2366 nint = pressure->GetExp(i)->GetNcoeffs();
2367 nbnd = fields[0]->GetExp(i)->NumBndryCoeffs();
2368 cnt1 = pressure->GetCoeff_Offset(i);
2369
2370 for (n = 0; n < nz_loc; ++n)
2371 {
2372 for (j = 0; j < nint; ++j)
2373 {
2374 p_coeffs[n * totpcoeffs + cnt1 + j] = f_p[cnt + j] =
2375 bnd[offset + nvel * nz_loc * nbnd + n * nint + j];
2376 }
2377 cnt += nint;
2378 }
2379 offset += (nvel * nbnd + nint) * nz_loc;
2380 }
2381
2382 // Back solve first level of static condensation for interior
2383 // velocity space and store in F_int
2384 F_int = F_int + Transpose(*m_mat[mode].m_D_int) * F_p -
2385 Transpose(*m_mat[mode].m_Btilde) * F_bnd;
2386 F_int = (*m_mat[mode].m_Cinv) * F_int;
2387
2388 // Unpack solution from Bnd and F_int to v_coeffs
2389 cnt = cnt1 = 0;
2390 for (i = 0; i < nel; ++i) // loop over elements
2391 {
2392 fields[0]->GetExp(i)->GetBoundaryMap(bmap);
2393 fields[0]->GetExp(i)->GetInteriorMap(imap);
2394 nbnd = bmap.size();
2395 nint = imap.size();
2396 offset = fields[0]->GetCoeff_Offset(i);
2397
2398 for (j = 0; j < nvel; ++j) // loop over velocity fields
2399 {
2400 for (n = 0; n < nz_loc; ++n)
2401 {
2402 for (k = 0; k < nbnd; ++k)
2403 {
2404 fields[j]->SetCoeff(n * nplanecoeffs + offset + bmap[k],
2405 f_bnd[cnt + k]);
2406 }
2407
2408 for (k = 0; k < nint; ++k)
2409 {
2410 fields[j]->SetCoeff(n * nplanecoeffs + offset + imap[k],
2411 f_int[cnt1 + k]);
2412 }
2413 cnt += nbnd;
2414 cnt1 += nint;
2415 }
2416 }
2417 }
2418
2419 for (j = 0; j < nvel; ++j)
2420 {
2421 fields[j]->SetPhysState(false);
2422 }
2423}
2424
2425void CoupledLinearNS::v_Output(void)
2426{
2427 std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size() + 1);
2428 std::vector<std::string> variables(m_fields.size() + 1);
2429 size_t i;
2430
2431 for (i = 0; i < m_fields.size(); ++i)
2432 {
2433 fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2434 variables[i] = m_boundaryConditions->GetVariable(i);
2435 }
2436
2437 fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2438 // project pressure field to velocity space
2439 if (m_singleMode == true)
2440 {
2441 Array<OneD, NekDouble> tmpfieldcoeffs(m_fields[0]->GetNcoeffs() / 2);
2442 m_pressure->GetPlane(0)->BwdTrans(
2443 m_pressure->GetPlane(0)->GetCoeffs(),
2444 m_pressure->GetPlane(0)->UpdatePhys());
2445 m_pressure->GetPlane(1)->BwdTrans(
2446 m_pressure->GetPlane(1)->GetCoeffs(),
2447 m_pressure->GetPlane(1)->UpdatePhys());
2448 m_fields[0]->GetPlane(0)->FwdTransLocalElmt(
2449 m_pressure->GetPlane(0)->GetPhys(), fieldcoeffs[i]);
2450 m_fields[0]->GetPlane(1)->FwdTransLocalElmt(
2451 m_pressure->GetPlane(1)->GetPhys(), tmpfieldcoeffs);
2452 for (int e = 0; e < m_fields[0]->GetNcoeffs() / 2; e++)
2453 {
2454 fieldcoeffs[i][e + m_fields[0]->GetNcoeffs() / 2] =
2455 tmpfieldcoeffs[e];
2456 }
2457 }
2458 else
2459 {
2460 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
2461 m_fields[0]->FwdTransLocalElmt(m_pressure->GetPhys(), fieldcoeffs[i]);
2462 }
2463 variables[i] = "p";
2464
2465 if (!m_comm->IsParallelInTime())
2466 {
2467 WriteFld(m_sessionName + ".fld", m_fields[0], fieldcoeffs, variables);
2468 }
2469 else
2470 {
2471 std::string newdir = m_sessionName + ".pit";
2472 if (!fs::is_directory(newdir))
2473 {
2474 fs::create_directory(newdir);
2475 }
2476 WriteFld(
2477 newdir + "/" + m_sessionName + "_" +
2478 std::to_string(m_windowPIT * m_comm->GetTimeComm()->GetSize() +
2479 m_comm->GetTimeComm()->GetRank() + 1) +
2480 ".fld",
2481 m_fields[0], fieldcoeffs, variables);
2482 }
2483}
2484
2485int CoupledLinearNS::v_GetForceDimension()
2486{
2487 return m_session->GetVariables().size();
2488}
2489} // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
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:1918
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
Print a summary of time stepping parameters.
Definition: CoupledLinearNS.cpp:1356
Nektar::CoupledLinearNS::EvaluateNewtonRHS
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: CoupledLinearNS.cpp:1933
Nektar::CoupledLinearNS::v_DoSolve
void v_DoSolve(void) override
Solves an unsteady problem.
Definition: CoupledLinearNS.cpp:1590
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:1528
Nektar::CoupledLinearNS::SolveLinearNS
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Definition: CoupledLinearNS.cpp:2091
Nektar::CoupledLinearNS::v_Output
void v_Output(void) override
Definition: CoupledLinearNS.cpp:2425
Nektar::CoupledLinearNS::v_DoInitialise
void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
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:1989
Nektar::CoupledLinearNS::v_NegatedOp
bool v_NegatedOp(void) override
Definition: CoupledLinearNS.cpp:1669
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:1515
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:1579
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:2485
Nektar::CoupledLinearNS::InfNorm
void InfNorm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: CoupledLinearNS.cpp:1901
Nektar::CoupledLinearNS::v_TransCoeffToPhys
void v_TransCoeffToPhys(void) override
Virtual function for transformation to physical space.
Definition: CoupledLinearNS.cpp:1568
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:197
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: BasicUtils/NekFactory.hpp:143
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::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:84
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:1253
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1315
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_HomogeneousType
enum HomogeneousType m_HomogeneousType
Definition: EquationSystem.h:556
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:741
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
CG_Iterations.pressure
pressure
Definition: CG_Iterations.py:229
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:268
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:143
Nektar::SpatialDomains::ExpansionInfoShPtr
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition: MeshGraph.h:140
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar::SpatialDomains::ExpansionInfoMap
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:141
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1794
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
Nektar::VarcoeffHashingTest::factors
StdRegions::ConstFactorMap factors
Definition: TestVarcoeffHashing.cpp:52
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:66
Nektar::eSteadyNavierStokes
@ eSteadyNavierStokes
Definition: IncNavierStokes.h:76
Nektar::eUnsteadyStokes
@ eUnsteadyStokes
Definition: IncNavierStokes.h:73
Nektar::eUnsteadyNavierStokes
@ eUnsteadyNavierStokes
Definition: IncNavierStokes.h:75
Nektar::eSteadyLinearisedNS
@ eSteadyLinearisedNS
Definition: IncNavierStokes.h:72
Nektar::eNoEquationType
@ eNoEquationType
Definition: IncNavierStokes.h:69
Nektar::eCoupledLinearisedNS
@ eCoupledLinearisedNS
Definition: IncNavierStokes.h:57
Nektar::NullNekDoubleArrayOfArray
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
Definition: BasicUtils/SharedArray.hpp:920
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
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
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 Jun 15 2025 07:52:34 for Nektar++ by doxygen 1.9.4