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