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