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