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