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