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 // 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 namespace Nektar
45 {
46 
47  string CoupledLinearNS::className = SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction("CoupledLinearisedNS", CoupledLinearNS::create);
48 
49 
50  /**
51  * @class CoupledLinearNS
52  *
53  * Set up expansion field for velocity and pressure, the local to
54  * global mapping arrays and the basic memory definitions for
55  * coupled matrix system
56  */
57  CoupledLinearNS::CoupledLinearNS(const LibUtilities::SessionReaderSharedPtr &pSession):
58  UnsteadySystem(pSession),
59  IncNavierStokes(pSession),
60  m_singleMode(false),
61  m_zeroMode(false)
62  {
63  }
64 
65  void CoupledLinearNS::v_InitObject()
66  {
67  IncNavierStokes::v_InitObject();
68 
69  int i;
70  int expdim = m_graph->GetMeshDimension();
71 
72  // Get Expansion list for orthogonal expansion at p-2
73  const SpatialDomains::ExpansionMap &pressure_exp = GenPressureExp(m_graph->GetExpansions("u"));
74 
75  m_nConvectiveFields = m_fields.num_elements();
76  if(boost::iequals(m_boundaryConditions->GetVariable(m_nConvectiveFields-1), "p"))
77  {
78  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");
79  }
80  // Decide how to declare explist for pressure.
81  if(expdim == 2)
82  {
83  int nz;
84 
85  if(m_HomogeneousType == eHomogeneous1D)
86  {
87  ASSERTL0(m_fields.num_elements() > 2,"Expect to have three at least three components of velocity variables");
88  LibUtilities::BasisKey Homo1DKey = m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
89 
90  m_pressure = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(m_session, Homo1DKey, m_LhomZ, m_useFFT,m_homogen_dealiasing, pressure_exp);
91 
92  ASSERTL1(m_npointsZ%2==0,"Non binary number of planes have been specified");
93  nz = m_npointsZ/2;
94 
95  }
96  else
97  {
98  //m_pressure2 = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(m_session, pressure_exp);
99  m_pressure = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(m_session, pressure_exp);
100  nz = 1;
101  }
102 
103  Array<OneD, MultiRegions::ExpListSharedPtr> velocity(m_velocity.num_elements());
104  for(i =0 ; i < m_velocity.num_elements(); ++i)
105  {
106  velocity[i] = m_fields[m_velocity[i]];
107  }
108 
109  // Set up Array of mappings
110  m_locToGloMap = Array<OneD, CoupledLocalToGlobalC0ContMapSharedPtr> (nz);
111 
112  if(m_session->DefinesSolverInfo("SingleMode"))
113  {
114 
115  ASSERTL0(nz <=2 ,"For single mode calculation can only have nz <= 2");
116  m_singleMode = true;
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 componont 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 
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 = m_fields[m_velocity[0]]->GetOffset_Elmt_Id(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  Timer timer;
460  timer.Start();
461  for(n = 0; n < nel; ++n)
462  {
463  eid = m_fields[m_velocity[0]]->GetOffset_Elmt_Id(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  EvaluateFunction(fieldStr,AdvField,"AdvectionVelocity");
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  EvaluateFunction(fieldStr, Restart, "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  EvaluateFunction(fieldStr,AdvField,"AdvectionVelocity");
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 convectioln terms
1358  EvaluateAdvectionTerms(inarray,outarray);
1359 
1360  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
1361  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
1362  {
1363  (*x)->Apply(m_fields, outarray, outarray, time);
1364  }
1365  }
1366 
1367  void CoupledLinearNS::SolveUnsteadyStokesSystem(const Array<OneD, const Array<OneD, NekDouble> > &inarray,
1368  Array<OneD, Array<OneD, NekDouble> > &outarray,
1369  const NekDouble time,
1370  const NekDouble aii_Dt)
1371  {
1372  int i;
1373  Array<OneD, Array< OneD, NekDouble> > F(m_nConvectiveFields);
1374  NekDouble lambda = 1.0/aii_Dt;
1375  static NekDouble lambda_store;
1376  Array <OneD, Array<OneD, NekDouble> > forcing(m_velocity.num_elements());
1377  // Matrix solution
1378  if(fabs(lambda_store - lambda) > 1e-10)
1379  {
1380  cout << "Setting up Stokes matrix problem [.";
1381  fflush(stdout);
1382  SetUpCoupledMatrix(lambda);
1383  cout << "]" << endl;
1384  lambda_store = lambda;
1385  }
1386 
1387  SetBoundaryConditions(time);
1388 
1389  // Forcing for advection solve
1390  for(i = 0; i < m_velocity.num_elements(); ++i)
1391  {
1392  m_fields[m_velocity[i]]->IProductWRTBase(inarray[i],m_fields[m_velocity[i]]->UpdateCoeffs());
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 
1495  void CoupledLinearNS::Solve(void)
1496  {
1497  const unsigned int ncmpt = m_velocity.num_elements();
1498  Array<OneD, Array<OneD, NekDouble> > forcing_phys(ncmpt);
1499  Array<OneD, Array<OneD, NekDouble> > forcing (ncmpt);
1500 
1501  for(int i = 0; i < ncmpt; ++i)
1502  {
1503  forcing_phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1504  forcing[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNcoeffs(),0.0);
1505  }
1506 
1507  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
1508  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
1509  {
1510  const NekDouble time = 0;
1511  (*x)->Apply(m_fields, forcing_phys, forcing_phys, time);
1512  }
1513  for (unsigned int i = 0; i < ncmpt; ++i)
1514  {
1515  m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1516  }
1517 
1518  SolveLinearNS(forcing);
1519  }
1520 
1521  void CoupledLinearNS::DefineForcingTerm(void)
1522  {
1523  m_ForcingTerm = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1524  m_ForcingTerm_Coeffs = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1525 
1526  for(int i = 0; i < m_velocity.num_elements(); ++i)
1527  {
1528  m_ForcingTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1529  m_ForcingTerm_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1530  }
1531 
1532  if(m_session->DefinesFunction("ForcingTerm"))
1533  {
1534  std::vector<std::string> fieldStr;
1535  for(int i = 0; i < m_velocity.num_elements(); ++i)
1536  {
1537  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1538  }
1539  EvaluateFunction(fieldStr, m_ForcingTerm, "ForcingTerm");
1540  for(int i = 0; i < m_velocity.num_elements(); ++i)
1541  {
1542  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(m_ForcingTerm[i], m_ForcingTerm_Coeffs[i]);
1543  }
1544  }
1545  else
1546  {
1547  cout << "'ForcingTerm' section has not been defined in the input file => forcing=0" << endl;
1548  }
1549  }
1550 
1551  void CoupledLinearNS::SolveSteadyNavierStokes(void)
1552  {
1553  Timer Newtontimer;
1554  Newtontimer.Start();
1555 
1556  Array<OneD, Array<OneD, NekDouble> > RHS_Coeffs(m_velocity.num_elements());
1557  Array<OneD, Array<OneD, NekDouble> > RHS_Phys(m_velocity.num_elements());
1558  Array<OneD, Array<OneD, NekDouble> > delta_velocity_Phys(m_velocity.num_elements());
1559  Array<OneD, Array<OneD, NekDouble> >Velocity_Phys(m_velocity.num_elements());
1560  Array<OneD, NekDouble > L2_norm(m_velocity.num_elements(), 1.0);
1561  Array<OneD, NekDouble > Inf_norm(m_velocity.num_elements(), 1.0);
1562 
1563  for(int i = 0; i < m_velocity.num_elements(); ++i)
1564  {
1565  delta_velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),1.0);
1566  Velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1567  }
1568 
1569  m_counter=1;
1570 
1571  L2Norm(delta_velocity_Phys, L2_norm);
1572 
1573  //while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1574  while(max(L2_norm[0], L2_norm[1]) > m_tol)
1575  {
1576  if(m_counter == 1) //At the first Newton step, we use the solution of the Stokes problem (if Restart=0 in input file)
1577  //Or the solution of the .rst file (if Restart=1 in input file)
1578  {
1579  for(int i = 0; i < m_velocity.num_elements(); ++i)
1580  {
1581  RHS_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1582  RHS_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1583  }
1584 
1585  for(int i = 0; i < m_velocity.num_elements(); ++i)
1586  {
1587  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1588  }
1589 
1590  m_initialStep = true;
1591  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1592  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1593  SolveLinearNS(RHS_Coeffs);
1594  m_initialStep = false;
1595  }
1596  if(m_counter > 1)
1597  {
1598  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1599 
1600  if(m_counter%m_MatrixSetUpStep == 0) //Setting Up the matrix is expensive. We do it at each "m_MatrixSetUpStep" step.
1601  {
1602  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1603  }
1604  SolveLinearNS(RHS_Coeffs);
1605  }
1606 
1607  for(int i = 0; i < m_velocity.num_elements(); ++i)
1608  {
1609  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(RHS_Coeffs[i], RHS_Phys[i]);
1610  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1611  }
1612 
1613  for(int i = 0; i < m_velocity.num_elements(); ++i)
1614  {
1615  Vmath::Vadd(Velocity_Phys[i].num_elements(),Velocity_Phys[i], 1, delta_velocity_Phys[i], 1,
1616  Velocity_Phys[i], 1);
1617  }
1618 
1619  //InfNorm(delta_velocity_Phys, Inf_norm);
1620  L2Norm(delta_velocity_Phys, L2_norm);
1621 
1622  if(max(Inf_norm[0], Inf_norm[1]) > 100)
1623  {
1624  cout<<"\nThe Newton method has failed at m_kinvis = "<<m_kinvis<<" (<=> Re = " << 1/m_kinvis << ")"<< endl;
1625  ASSERTL0(0, "The Newton method has failed... \n");
1626  }
1627 
1628 
1629  cout << "\n";
1630  m_counter++;
1631  }
1632 
1633  if (m_counter > 1) //We save u:=u+\delta u in u->Coeffs
1634  {
1635  for(int i = 0; i < m_velocity.num_elements(); ++i)
1636  {
1637  m_fields[m_velocity[i]]->FwdTrans(Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1638  }
1639  }
1640 
1641  Newtontimer.Stop();
1642  cout<<"We have done "<< m_counter-1 << " iteration(s) in " << Newtontimer.TimePerTest(1)/60 << " minute(s). \n\n";
1643  }
1644 
1645 
1646  void CoupledLinearNS::Continuation(void)
1647  {
1648  Array<OneD, Array<OneD, NekDouble> > u_N(m_velocity.num_elements());
1649  Array<OneD, Array<OneD, NekDouble> > tmp_RHS(m_velocity.num_elements());
1650  Array<OneD, Array<OneD, NekDouble> > RHS(m_velocity.num_elements());
1651  Array<OneD, Array<OneD, NekDouble> > u_star(m_velocity.num_elements());
1652 
1653  cout << "We apply the continuation method: " <<endl;
1654 
1655  for(int i = 0; i < m_velocity.num_elements(); ++i)
1656  {
1657  u_N[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1658  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_N[i]);
1659 
1660  RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1661  tmp_RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1662 
1663  m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1664  Vmath::Smul(tmp_RHS[i].num_elements(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1665 
1666  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1667  }
1668 
1669  SetUpCoupledMatrix(0.0, u_N, true);
1670  SolveLinearNS(RHS);
1671 
1672  for(int i = 0; i < m_velocity.num_elements(); ++i)
1673  {
1674  u_star[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1675  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_star[i]);
1676 
1677  //u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1678  Vmath::Smul(u_star[i].num_elements(), m_kinvis, u_star[i], 1, u_star[i], 1);
1679  Vmath::Vadd(u_star[i].num_elements(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1680 
1681  m_fields[m_velocity[i]]->FwdTrans(u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1682  }
1683 
1684  m_kinvis -= m_kinvis*m_KinvisPercentage/100;
1685  }
1686 
1687 
1688  void CoupledLinearNS::InfNorm(Array<OneD, Array<OneD, NekDouble> > &inarray,
1689  Array<OneD, NekDouble> &outarray)
1690  {
1691  for(int i = 0; i < m_velocity.num_elements(); ++i)
1692  {
1693  outarray[i] = 0.0;
1694  for(int j = 0; j < inarray[i].num_elements(); ++j)
1695  {
1696  if(inarray[i][j] > outarray[i])
1697  {
1698  outarray[i] = inarray[i][j];
1699  }
1700  }
1701  cout << "InfNorm["<<i<<"] = "<< outarray[i] <<endl;
1702  }
1703  }
1704 
1705  void CoupledLinearNS::L2Norm(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  outarray[i] += inarray[i][j]*inarray[i][j];
1714  }
1715  outarray[i]=sqrt(outarray[i]);
1716  cout << "L2Norm["<<i<<"] = "<< outarray[i] <<endl;
1717  }
1718  }
1719 
1720 
1721  void CoupledLinearNS::EvaluateNewtonRHS(Array<OneD, Array<OneD, NekDouble> > &Velocity,
1722  Array<OneD, Array<OneD, NekDouble> > &outarray)
1723  {
1724  Array<OneD, Array<OneD, NekDouble> > Eval_Adv(m_velocity.num_elements());
1725  Array<OneD, Array<OneD, NekDouble> > tmp_DerVel(m_velocity.num_elements());
1726  Array<OneD, Array<OneD, NekDouble> > AdvTerm(m_velocity.num_elements());
1727  Array<OneD, Array<OneD, NekDouble> > ViscTerm(m_velocity.num_elements());
1728  Array<OneD, Array<OneD, NekDouble> > Forc(m_velocity.num_elements());
1729 
1730  for(int i = 0; i < m_velocity.num_elements(); ++i)
1731  {
1732  Eval_Adv[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1733  tmp_DerVel[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1734 
1735  AdvTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1736  ViscTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1737  Forc[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1738  outarray[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1739 
1740  m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1741 
1742  Vmath::Smul(tmp_DerVel[i].num_elements(), m_kinvis, tmp_DerVel[i], 1, tmp_DerVel[i], 1);
1743  }
1744 
1745  EvaluateAdvectionTerms(Velocity, Eval_Adv);
1746 
1747  for(int i = 0; i < m_velocity.num_elements(); ++i)
1748  {
1749  m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i], AdvTerm[i]); //(w, (u.grad)u)
1750  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1751  m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i], Forc[i]); //(w, f)
1752 
1753  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, AdvTerm[i], 1, outarray[i], 1);
1754  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, ViscTerm[i], 1, outarray[i], 1);
1755 
1756  Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1, Forc[i], 1, outarray[i], 1);
1757  }
1758  }
1759 
1760 
1761 
1762  const SpatialDomains::ExpansionMap &CoupledLinearNS::GenPressureExp(const SpatialDomains::ExpansionMap &VelExp)
1763  {
1764  int i;
1765  SpatialDomains::ExpansionMapShPtr returnval;
1766 
1767  returnval = MemoryManager<SpatialDomains::ExpansionMap>::AllocateSharedPtr();
1768 
1769  SpatialDomains::ExpansionMap::const_iterator expMapIter;
1770  int nummodes;
1771 
1772  for (expMapIter = VelExp.begin(); expMapIter != VelExp.end(); ++expMapIter)
1773  {
1774  LibUtilities::BasisKeyVector BasisVec;
1775 
1776  for(i = 0; i < expMapIter->second->m_basisKeyVector.size(); ++i)
1777  {
1778  LibUtilities::BasisKey B = expMapIter->second->m_basisKeyVector[i];
1779  nummodes = B.GetNumModes();
1780  ASSERTL0(nummodes > 3,"Velocity polynomial space not sufficiently high (>= 4)");
1781  // Should probably set to be an orthogonal basis.
1782  LibUtilities::BasisKey newB(B.GetBasisType(),nummodes-2,B.GetPointsKey());
1783  BasisVec.push_back(newB);
1784  }
1785 
1786  // Put new expansion into list.
1787  SpatialDomains::ExpansionShPtr expansionElementShPtr =
1788  MemoryManager<SpatialDomains::Expansion>::AllocateSharedPtr(expMapIter->second->m_geomShPtr, BasisVec);
1789  (*returnval)[expMapIter->first] = expansionElementShPtr;
1790  }
1791 
1792  // Save expansion into graph.
1793  m_graph->SetExpansions("p",returnval);
1794 
1795  return *returnval;
1796  }
1797 
1798  /**
1799  * @param forcing A list of forcing functions for each velocity
1800  * component
1801  *
1802  * The routine involves two levels of static
1803  * condensations. Initially we require a statically condensed
1804  * forcing function which requires the following manipulation
1805  *
1806  * \f[ {F\_bnd} = {\bf f}_{bnd} -m\_B \,m\_Cinv\, {\bf f}_{int},
1807  * \hspace{1cm} F\_p = m\_D\_{int}\, m\_Cinv\, {\bf f}_{int} \f]
1808  *
1809  * Where \f${\bf f}_{bnd}\f$ denote the forcing degrees of
1810  * freedom of the elemental velocities on the boundary of the
1811  * element, \f${\bf f}_{int}\f$ denote the forcing degrees of
1812  * freedom of the elemental velocities on the interior of the
1813  * element. (see detailed description for more details).
1814  *
1815  * This vector is further manipulated into
1816  *
1817  * \f[ Fh\_{bnd} = \left [ \begin{array}{c} f\_{bnd} -m\_B \,
1818  * m\_Cinv\, {\bf f}_{int}\\ \left [m\_D\_{int} \, m\_Cinv \,{\bf
1819  * f}_{int} \right]_0 \end{array}\right ] \hspace{1cm} [Fh\_p]_{i} =
1820  * \begin{array}{c} [m\_D\_{int} \, m\_Cinv \, {\bf
1821  * f}_{int}]_{i+1} \end{array} \f]
1822  *
1823  * where \f$-{[}m\_D\_{int}^T\, m\_Cinv \,{\bf f}_{int}]_0\f$
1824  * which is corresponds to the mean mode of the pressure degrees
1825  * of freedom and is now added to the boundary system and the
1826  * remainder of the block becomes the interior forcing for the
1827  * inner static condensation (see detailed description for more
1828  * details) which is set up in a GlobalLinSysDirectStaticCond
1829  * class.
1830  *
1831  * Finally we perform the final maniplation of the forcing to
1832  * using hte
1833  * \f[ Fh\_{bnd} = Fh\_{bnd} - m\_Bh \,m\_Chinv \, Fh\_p \f]
1834  *
1835  * We can now call the solver to the global coupled boundary
1836  * system (through the call to #m_CoupledBndSys->Solve) to obtain
1837  * the velocity boundary solution as the mean pressure solution,
1838  * i.e.
1839  *
1840  * \f[ {\cal A}^T(\hat{A} - \hat{C}^T \hat{D}^{-1} \hat{B} ){\cal
1841  * A} \, Bnd = Fh\_{bnd} \f]
1842  *
1843  * Once we know the solution to the above the rest of the pressure
1844  * modes are recoverable thorugh
1845  *
1846  * \f[ Ph = m\_Dhinv\, (Bnd - m\_Ch^T \, Fh_{bnd}) \f]
1847  *
1848  * We can now unpack \f$ Fh\_{bnd} \f$ (last elemental mode) and
1849  * \f$ Ph \f$ into #m_pressure and \f$ F_p\f$ and \f$ Fh\_{bnd}\f$
1850  * into a closed pack list of boundary velocoity degrees of
1851  * freedom stored in \f$ F\_bnd\f$.
1852  *
1853  * Finally the interior velocity degrees of freedom are then
1854  * obtained through the relationship
1855  *
1856  * \f[ F\_{int} = m\_Cinv\ ( F\_{int} + m\_D\_{int}^T\,
1857  * F\_p - m\_Btilde^T\, Bnd) \f]
1858  *
1859  * We then unpack the solution back to the MultiRegion structures
1860  * #m_velocity and #m_pressure
1861  */
1862  void CoupledLinearNS::SolveLinearNS(const Array<OneD, Array<OneD, NekDouble> > &forcing)
1863  {
1864  int i,n;
1865  Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.num_elements());
1866  Array<OneD, Array<OneD, NekDouble> > force(m_velocity.num_elements());
1867 
1868  if(m_HomogeneousType == eHomogeneous1D)
1869  {
1870  int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
1871  for(n = 0; n < m_npointsZ/2; ++n)
1872  {
1873  // Get the a Fourier mode of velocity and forcing components.
1874  for(i = 0; i < m_velocity.num_elements(); ++i)
1875  {
1876  vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2*n);
1877  // Note this needs to correlate with how we pass forcing
1878  force[i] = forcing[i] + 2*n*ncoeffsplane;
1879  }
1880 
1881  SolveLinearNS(force,vel_fields,m_pressure->GetPlane(2*n),n);
1882  }
1883  for(i = 0; i < m_velocity.num_elements(); ++i)
1884  {
1885  m_fields[m_velocity[i]]->SetPhysState(false);
1886  }
1887  m_pressure->SetPhysState(false);
1888  }
1889  else
1890  {
1891  for(i = 0; i < m_velocity.num_elements(); ++i)
1892  {
1893  vel_fields[i] = m_fields[m_velocity[i]];
1894  // Note this needs to correlate with how we pass forcing
1895  force[i] = forcing[i];
1896  }
1897  SolveLinearNS(force,vel_fields,m_pressure);
1898  }
1899  }
1900 
1901  void CoupledLinearNS::SolveLinearNS(const Array<OneD, Array<OneD, NekDouble> > &forcing, Array<OneD, MultiRegions::ExpListSharedPtr> &fields, MultiRegions::ExpListSharedPtr &pressure, const int mode)
1902  {
1903  int i,j,k,n,eid,cnt,cnt1;
1904  int nbnd,nint,offset;
1905  int nvel = m_velocity.num_elements();
1906  int nel = fields[0]->GetNumElmts();
1907  Array<OneD, unsigned int> bmap, imap;
1908 
1909  Array<OneD, NekDouble > f_bnd(m_mat[mode].m_BCinv->GetRows());
1910  NekVector< NekDouble > F_bnd(f_bnd.num_elements(), f_bnd, eWrapper);
1911  Array<OneD, NekDouble > f_int(m_mat[mode].m_BCinv->GetColumns());
1912  NekVector< NekDouble > F_int(f_int.num_elements(),f_int, eWrapper);
1913 
1914  int nz_loc;
1915  int nplanecoeffs = fields[m_velocity[0]]->GetNcoeffs();// this is fine since we pass the nplane coeff data.
1916 
1917  if(mode) // Homogeneous mode flag
1918  {
1919  nz_loc = 2;
1920  }
1921  else
1922  {
1923  if(m_singleMode)
1924  {
1925  nz_loc = 2;
1926  }
1927  else
1928  {
1929  nz_loc = 1;
1930  if(m_HomogeneousType == eHomogeneous1D)
1931  {
1932  // Zero fields to set complex mode to zero;
1933  for(i = 0; i < fields.num_elements(); ++i)
1934  {
1935  Vmath::Zero(2*fields[i]->GetNcoeffs(),fields[i]->UpdateCoeffs(),1);
1936  }
1937  Vmath::Zero(2*pressure->GetNcoeffs(),pressure->UpdateCoeffs(),1);
1938  }
1939  }
1940  }
1941 
1942  // Assemble f_bnd and f_int
1943  cnt = cnt1 = 0;
1944  for(i = 0; i < nel; ++i) // loop over elements
1945  {
1946  eid = fields[m_velocity[0]]->GetOffset_Elmt_Id(i);
1947  fields[m_velocity[0]]->GetExp(eid)->GetBoundaryMap(bmap);
1948  fields[m_velocity[0]]->GetExp(eid)->GetInteriorMap(imap);
1949  nbnd = bmap.num_elements();
1950  nint = imap.num_elements();
1951  offset = fields[m_velocity[0]]->GetCoeff_Offset(eid);
1952 
1953  for(j = 0; j < nvel; ++j) // loop over velocity fields
1954  {
1955  for(n = 0; n < nz_loc; ++n)
1956  {
1957  for(k = 0; k < nbnd; ++k)
1958  {
1959  f_bnd[cnt+k] = forcing[j][n*nplanecoeffs +
1960  offset+bmap[k]];
1961  }
1962  for(k = 0; k < nint; ++k)
1963  {
1964  f_int[cnt1+k] = forcing[j][n*nplanecoeffs +
1965  offset+imap[k]];
1966  }
1967  cnt += nbnd;
1968  cnt1 += nint;
1969  }
1970  }
1971  }
1972 
1973  Array<OneD, NekDouble > f_p(m_mat[mode].m_D_int->GetRows());
1974  NekVector< NekDouble > F_p(f_p.num_elements(),f_p,eWrapper);
1975  NekVector< NekDouble > F_p_tmp(m_mat[mode].m_Cinv->GetRows());
1976 
1977  // fbnd does not currently hold the pressure mean
1978  F_bnd = F_bnd - (*m_mat[mode].m_BCinv)*F_int;
1979  F_p_tmp = (*m_mat[mode].m_Cinv)*F_int;
1980  F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
1981 
1982  // construct inner forcing
1983  Array<OneD, NekDouble > bnd (m_locToGloMap[mode]->GetNumGlobalCoeffs(),0.0);
1984  Array<OneD, NekDouble > fh_bnd(m_locToGloMap[mode]->GetNumGlobalCoeffs(),0.0);
1985 
1986  const Array<OneD,const int>& loctoglomap
1987  = m_locToGloMap[mode]->GetLocalToGlobalMap();
1988  const Array<OneD,const NekDouble>& loctoglosign
1989  = m_locToGloMap[mode]->GetLocalToGlobalSign();
1990 
1991  offset = cnt = 0;
1992  for(i = 0; i < nel; ++i)
1993  {
1994  eid = fields[0]->GetOffset_Elmt_Id(i);
1995  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
1996 
1997  for(j = 0; j < nvel; ++j)
1998  {
1999  for(k = 0; k < nbnd; ++k)
2000  {
2001  fh_bnd[loctoglomap[offset+j*nbnd+k]] +=
2002  loctoglosign[offset+j*nbnd+k]*f_bnd[cnt+k];
2003  }
2004  cnt += nbnd;
2005  }
2006 
2007  nint = pressure->GetExp(eid)->GetNcoeffs();
2008  offset += nvel*nbnd + nint*nz_loc;
2009  }
2010 
2011  offset = cnt1 = 0;
2012  for(i = 0; i < nel; ++i)
2013  {
2014  eid = fields[0]->GetOffset_Elmt_Id(i);
2015  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2016  nint = pressure->GetExp(eid)->GetNcoeffs();
2017 
2018  for(n = 0; n < nz_loc; ++n)
2019  {
2020  for(j = 0; j < nint; ++j)
2021  {
2022  fh_bnd[loctoglomap[offset + nvel*nbnd + n*nint+j]] = f_p[cnt1+j];
2023  }
2024  cnt1 += nint;
2025  }
2026  offset += nvel*nbnd + nz_loc*nint;
2027  }
2028 
2029  // Set Weak BC into f_bnd and Dirichlet Dofs in bnd
2030  const Array<OneD,const int>& bndmap
2031  = m_locToGloMap[mode]->GetBndCondCoeffsToGlobalCoeffsMap();
2032 
2033  // Forcing function with weak boundary conditions and
2034  // Dirichlet conditions
2035  int bndcnt=0;
2036 
2037  for(k = 0; k < nvel; ++k)
2038  {
2039  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConds = fields[k]->GetBndConditions();
2040  Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp;
2041  if(m_HomogeneousType == eHomogeneous1D)
2042  {
2043  bndCondExp = m_fields[k]->GetPlane(2*mode)->GetBndCondExpansions();
2044  }
2045  else
2046  {
2047  bndCondExp = m_fields[k]->GetBndCondExpansions();
2048  }
2049 
2050  for(i = 0; i < bndCondExp.num_elements(); ++i)
2051  {
2052  const Array<OneD, const NekDouble > bndCondCoeffs = bndCondExp[i]->GetCoeffs();
2053  cnt = 0;
2054  for(n = 0; n < nz_loc; ++n)
2055  {
2056  if(bndConds[i]->GetBoundaryConditionType()
2057  == SpatialDomains::eDirichlet)
2058  {
2059  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
2060  {
2061  if (m_equationType == eSteadyNavierStokes && m_initialStep == false)
2062  {
2063  //This condition set all the Dirichlet BC at 0 after
2064  //the initial step of the Newton method
2065  bnd[bndmap[bndcnt++]] = 0;
2066  }
2067  else
2068  {
2069  bnd[bndmap[bndcnt++]] = bndCondCoeffs[cnt++];
2070  }
2071  }
2072  }
2073  else
2074  {
2075  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
2076  {
2077  fh_bnd[bndmap[bndcnt++]]
2078  += bndCondCoeffs[cnt++];
2079  }
2080  }
2081  }
2082  }
2083  }
2084 
2085  m_mat[mode].m_CoupledBndSys->Solve(fh_bnd,bnd,m_locToGloMap[mode]);
2086 
2087  // unpack pressure and velocity boundary systems.
2088  offset = cnt = 0;
2089  int totpcoeffs = pressure->GetNcoeffs();
2090  Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2091  for(i = 0; i < nel; ++i)
2092  {
2093  eid = fields[0]->GetOffset_Elmt_Id(i);
2094  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2095  nint = pressure->GetExp(eid)->GetNcoeffs();
2096 
2097  for(j = 0; j < nvel; ++j)
2098  {
2099  for(k = 0; k < nbnd; ++k)
2100  {
2101  f_bnd[cnt+k] = loctoglosign[offset+j*nbnd+k]*bnd[loctoglomap[offset + j*nbnd + k]];
2102  }
2103  cnt += nbnd;
2104  }
2105  offset += nvel*nbnd + nint*nz_loc;
2106  }
2107 
2108  pressure->SetPhysState(false);
2109 
2110  offset = cnt = cnt1 = 0;
2111  for(i = 0; i < nel; ++i)
2112  {
2113  eid = fields[0]->GetOffset_Elmt_Id(i);
2114  nint = pressure->GetExp(eid)->GetNcoeffs();
2115  nbnd = fields[0]->GetExp(eid)->NumBndryCoeffs();
2116  cnt1 = pressure->GetCoeff_Offset(eid);
2117 
2118  for(n = 0; n < nz_loc; ++n)
2119  {
2120  for(j = 0; j < nint; ++j)
2121  {
2122  p_coeffs[n*totpcoeffs + cnt1+j] =
2123  f_p[cnt+j] = bnd[loctoglomap[offset +
2124  (nvel*nz_loc)*nbnd +
2125  n*nint + j]];
2126  }
2127  cnt += nint;
2128  }
2129  offset += (nvel*nbnd + nint)*nz_loc;
2130  }
2131 
2132  // Back solve first level of static condensation for interior
2133  // velocity space and store in F_int
2134  F_int = F_int + Transpose(*m_mat[mode].m_D_int)*F_p
2135  - Transpose(*m_mat[mode].m_Btilde)*F_bnd;
2136  F_int = (*m_mat[mode].m_Cinv)*F_int;
2137 
2138  // Unpack solution from Bnd and F_int to v_coeffs
2139  cnt = cnt1 = 0;
2140  for(i = 0; i < nel; ++i) // loop over elements
2141  {
2142  eid = fields[m_velocity[0]]->GetOffset_Elmt_Id(i);
2143  fields[0]->GetExp(eid)->GetBoundaryMap(bmap);
2144  fields[0]->GetExp(eid)->GetInteriorMap(imap);
2145  nbnd = bmap.num_elements();
2146  nint = imap.num_elements();
2147  offset = fields[0]->GetCoeff_Offset(eid);
2148 
2149  for(j = 0; j < nvel; ++j) // loop over velocity fields
2150  {
2151  for(n = 0; n < nz_loc; ++n)
2152  {
2153  for(k = 0; k < nbnd; ++k)
2154  {
2155  fields[j]->SetCoeff(n*nplanecoeffs +
2156  offset+bmap[k],
2157  f_bnd[cnt+k]);
2158  }
2159 
2160  for(k = 0; k < nint; ++k)
2161  {
2162  fields[j]->SetCoeff(n*nplanecoeffs +
2163  offset+imap[k],
2164  f_int[cnt1+k]);
2165  }
2166  cnt += nbnd;
2167  cnt1 += nint;
2168  }
2169  }
2170  }
2171 
2172  for(j = 0; j < nvel; ++j)
2173  {
2174  fields[j]->SetPhysState(false);
2175  }
2176  }
2177 
2178  void CoupledLinearNS::v_Output(void)
2179  {
2180  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_fields.num_elements()+1);
2181  std::vector<std::string> variables(m_fields.num_elements()+1);
2182  int i;
2183 
2184  for(i = 0; i < m_fields.num_elements(); ++i)
2185  {
2186  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2187  variables[i] = m_boundaryConditions->GetVariable(i);
2188  }
2189 
2190  fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2191  // project pressure field to velocity space
2192  if(m_singleMode==true)
2193  {
2194  Array<OneD, NekDouble > tmpfieldcoeffs (m_fields[0]->GetNcoeffs()/2);
2195  m_pressure->GetPlane(0)->BwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetCoeffs(), m_pressure->GetPlane(0)->UpdatePhys());
2196  m_pressure->GetPlane(1)->BwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetCoeffs(), m_pressure->GetPlane(1)->UpdatePhys());
2197  m_fields[0]->GetPlane(0)->FwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetPhys(),fieldcoeffs[i]);
2198  m_fields[0]->GetPlane(1)->FwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetPhys(),tmpfieldcoeffs);
2199  for(int e=0; e<m_fields[0]->GetNcoeffs()/2; e++)
2200  {
2201  fieldcoeffs[i][e+m_fields[0]->GetNcoeffs()/2] = tmpfieldcoeffs[e];
2202  }
2203  }
2204  else
2205  {
2206  m_pressure->BwdTrans_IterPerExp(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
2207  m_fields[0]->FwdTrans_IterPerExp(m_pressure->GetPhys(),fieldcoeffs[i]);
2208  }
2209  variables[i] = "p";
2210 
2211  std::string outname = m_sessionName + ".fld";
2212 
2213  WriteFld(outname,m_fields[0],fieldcoeffs,variables);
2214  }
2215 
2216  int CoupledLinearNS::v_GetForceDimension()
2217  {
2218  return m_session->GetVariables().size();
2219  }
2220 }
2221 
2222 /**
2223  * $Log: CoupledLinearNS.cpp,v $
2224  **/
Nektar::IncNavierStokes::m_equationType
EquationType m_equationType
equation type;
Definition: IncNavierStokes.h:179
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:135
TimeIntegrationWrapper.h
Nektar::IncNavierStokes::SetBoundaryConditions
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
Definition: IncNavierStokes.cpp:360
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:168
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineImplicitSolve
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationScheme.h:202
Nektar::GetExtrapolateFactory
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:48
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:67
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:438
Nektar::SolverUtils::EquationSystem::m_LhomZ
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition: EquationSystem.h:498
Nektar::CoupledLinearNS::v_Output
virtual void v_Output(void)
Definition: CoupledLinearNS.cpp:2178
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:1762
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:163
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1953
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:143
Nektar::CoupledLinearNS::m_mat
Array< OneD, CoupledSolverMatrices > m_mat
Definition: CoupledLinearNS.h:187
Nektar::SolverUtils::EquationSystem::m_useFFT
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
Definition: EquationSystem.h:456
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
Nektar::CoupledLinearNS::L2Norm
void L2Norm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: CoupledLinearNS.cpp:1705
Nektar::SolverUtils::EquationSystem::m_npointsZ
int m_npointsZ
number of points in Z direction (if homogeneous)
Definition: EquationSystem.h:502
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:430
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:159
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
Nektar::CoupledLinearNS::v_TransCoeffToPhys
virtual void v_TransCoeffToPhys(void)
Virtual function for transformation to physical space.
Definition: CoupledLinearNS.cpp:1406
Nektar::CoupledLinearNS::m_singleMode
bool m_singleMode
Identify if a single mode is required for stability analysis.
Definition: CoupledLinearNS.h:170
Nektar::DNekScalMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:73
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:824
Nektar::CoupledLinearNS::EvaluateNewtonRHS
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: CoupledLinearNS.cpp:1721
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:65
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:199
Timer.h
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationScheme.h:178
Nektar::SpatialDomains::ExpansionMapShPtr
boost::shared_ptr< ExpansionMap > ExpansionMapShPtr
Definition: MeshGraph.h:175
Nektar::CoupledLinearNS::SolveLinearNS
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Definition: CoupledLinearNS.cpp:1862
Nektar::IncNavierStokes::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Definition: IncNavierStokes.h:156
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
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:1208
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:426
CoupledLinearNS.h
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:221
Nektar::CoupledLinearNS::v_GenerateSummary
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Definition: CoupledLinearNS.cpp:1203
Nektar::CoupledLinearNS::InfNorm
void InfNorm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: CoupledLinearNS.cpp:1688
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::LibUtilities::BasisKeyVector
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
Definition: FoundationsFwd.hpp:69
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:2216
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:684
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:461
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:81
MatrixKey.h
Nektar::CoupledLinearNS::className
static std::string className
Name of class.
Definition: CoupledLinearNS.h:103
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:329
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:1993
Nektar::CoupledLinearNS::m_locToGloMap
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
Definition: CoupledLinearNS.h:161
Nektar::IncNavierStokes::EvaluateAdvectionTerms
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
Definition: IncNavierStokes.cpp:318
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:1353
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:420
Nektar::SpatialDomains::ExpansionShPtr
boost::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:170
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:412
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:768
Nektar::CoupledLinearNS::create
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
Definition: CoupledLinearNS.h:95
Nektar::CoupledLinearNS::m_zeroMode
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
Definition: CoupledLinearNS.h:172
Nektar::IncNavierStokes
This class is the base class for Navier Stokes problems.
Definition: IncNavierStokes.h:104
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition: EquationSystem.h:428
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:742
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:166
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:785
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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:1367
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1791
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:165
Nektar::IncNavierStokes::v_InitObject
virtual void v_InitObject()
Init object for UnsteadySystem class.
Definition: IncNavierStokes.cpp:65
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:494
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:285
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:169
Nektar::SpatialDomains::ExpansionMap
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:171
Nektar::coupledSolverMatrices::m_BCinv
DNekScalBlkMatSharedPtr m_BCinv
Boundary-interior Laplacian plus linearised convective terms pre-multiplying Cinv: ...
Definition: CoupledLinearNS.h:59
Nektar::CoupledLinearNS::CoupledLinearNS
CoupledLinearNS(const LibUtilities::SessionReaderSharedPtr &pSesssion)
Definition: CoupledLinearNS.cpp:57
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 Tue Dec 8 2015 23:41:27 for Nektar++ by   doxygen 1.8.9.1