Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CoupledLinearNS.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File CoupledLinearNS.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Coupled Solver for the Linearised Incompressible
33 // Navier Stokes equations
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/algorithm/string.hpp>
37 
38 #include <LibUtilities/TimeIntegration/TimeIntegrationWrapper.h>
39 #include <IncNavierStokesSolver/EquationSystems/CoupledLinearNS.h>
40 #include <LibUtilities/BasicUtils/Timer.h>
41 #include <LocalRegions/MatrixKey.h>
42 #include <MultiRegions/GlobalLinSysDirectStaticCond.h>
43 
44 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  if(m_HomogeneousType == eHomogeneous1D)
1517  {
1518  if(!m_singleMode) // Do not want to Fourier transoform in case of single mode stability an
1519  {
1520  m_fields[i]->HomogeneousFwdTrans(forcing[i], forcing[i]);
1521  }
1522  }
1523  }
1524 
1525  SolveLinearNS(forcing);
1526  }
1527 
1528  void CoupledLinearNS::DefineForcingTerm(void)
1529  {
1530  m_ForcingTerm = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1531  m_ForcingTerm_Coeffs = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1532 
1533  for(int i = 0; i < m_velocity.num_elements(); ++i)
1534  {
1535  m_ForcingTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1536  m_ForcingTerm_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1537  }
1538 
1539  if(m_session->DefinesFunction("ForcingTerm"))
1540  {
1541  std::vector<std::string> fieldStr;
1542  for(int i = 0; i < m_velocity.num_elements(); ++i)
1543  {
1544  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1545  }
1546  EvaluateFunction(fieldStr, m_ForcingTerm, "ForcingTerm");
1547  for(int i = 0; i < m_velocity.num_elements(); ++i)
1548  {
1549  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(m_ForcingTerm[i], m_ForcingTerm_Coeffs[i]);
1550  }
1551  }
1552  else
1553  {
1554  cout << "'ForcingTerm' section has not been defined in the input file => forcing=0" << endl;
1555  }
1556  }
1557 
1558  void CoupledLinearNS::SolveSteadyNavierStokes(void)
1559  {
1560  Timer Newtontimer;
1561  Newtontimer.Start();
1562 
1563  Array<OneD, Array<OneD, NekDouble> > RHS_Coeffs(m_velocity.num_elements());
1564  Array<OneD, Array<OneD, NekDouble> > RHS_Phys(m_velocity.num_elements());
1565  Array<OneD, Array<OneD, NekDouble> > delta_velocity_Phys(m_velocity.num_elements());
1566  Array<OneD, Array<OneD, NekDouble> >Velocity_Phys(m_velocity.num_elements());
1567  Array<OneD, NekDouble > L2_norm(m_velocity.num_elements(), 1.0);
1568  Array<OneD, NekDouble > Inf_norm(m_velocity.num_elements(), 1.0);
1569 
1570  for(int i = 0; i < m_velocity.num_elements(); ++i)
1571  {
1572  delta_velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),1.0);
1573  Velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1574  }
1575 
1576  m_counter=1;
1577 
1578  L2Norm(delta_velocity_Phys, L2_norm);
1579 
1580  //while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1581  while(max(L2_norm[0], L2_norm[1]) > m_tol)
1582  {
1583  if(m_counter == 1) //At the first Newton step, we use the solution of the Stokes problem (if Restart=0 in input file)
1584  //Or the solution of the .rst file (if Restart=1 in input file)
1585  {
1586  for(int i = 0; i < m_velocity.num_elements(); ++i)
1587  {
1588  RHS_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1589  RHS_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1590  }
1591 
1592  for(int i = 0; i < m_velocity.num_elements(); ++i)
1593  {
1594  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1595  }
1596 
1597  m_initialStep = true;
1598  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1599  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1600  SolveLinearNS(RHS_Coeffs);
1601  m_initialStep = false;
1602  }
1603  if(m_counter > 1)
1604  {
1605  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1606 
1607  if(m_counter%m_MatrixSetUpStep == 0) //Setting Up the matrix is expensive. We do it at each "m_MatrixSetUpStep" step.
1608  {
1609  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1610  }
1611  SolveLinearNS(RHS_Coeffs);
1612  }
1613 
1614  for(int i = 0; i < m_velocity.num_elements(); ++i)
1615  {
1616  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(RHS_Coeffs[i], RHS_Phys[i]);
1617  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1618  }
1619 
1620  for(int i = 0; i < m_velocity.num_elements(); ++i)
1621  {
1622  Vmath::Vadd(Velocity_Phys[i].num_elements(),Velocity_Phys[i], 1, delta_velocity_Phys[i], 1,
1623  Velocity_Phys[i], 1);
1624  }
1625 
1626  //InfNorm(delta_velocity_Phys, Inf_norm);
1627  L2Norm(delta_velocity_Phys, L2_norm);
1628 
1629  if(max(Inf_norm[0], Inf_norm[1]) > 100)
1630  {
1631  cout<<"\nThe Newton method has failed at m_kinvis = "<<m_kinvis<<" (<=> Re = " << 1/m_kinvis << ")"<< endl;
1632  ASSERTL0(0, "The Newton method has failed... \n");
1633  }
1634 
1635 
1636  cout << "\n";
1637  m_counter++;
1638  }
1639 
1640  if (m_counter > 1) //We save u:=u+\delta u in u->Coeffs
1641  {
1642  for(int i = 0; i < m_velocity.num_elements(); ++i)
1643  {
1644  m_fields[m_velocity[i]]->FwdTrans(Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1645  }
1646  }
1647 
1648  Newtontimer.Stop();
1649  cout<<"We have done "<< m_counter-1 << " iteration(s) in " << Newtontimer.TimePerTest(1)/60 << " minute(s). \n\n";
1650  }
1651 
1652 
1653  void CoupledLinearNS::Continuation(void)
1654  {
1655  Array<OneD, Array<OneD, NekDouble> > u_N(m_velocity.num_elements());
1656  Array<OneD, Array<OneD, NekDouble> > tmp_RHS(m_velocity.num_elements());
1657  Array<OneD, Array<OneD, NekDouble> > RHS(m_velocity.num_elements());
1658  Array<OneD, Array<OneD, NekDouble> > u_star(m_velocity.num_elements());
1659 
1660  cout << "We apply the continuation method: " <<endl;
1661 
1662  for(int i = 0; i < m_velocity.num_elements(); ++i)
1663  {
1664  u_N[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1665  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_N[i]);
1666 
1667  RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1668  tmp_RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1669 
1670  m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1671  Vmath::Smul(tmp_RHS[i].num_elements(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1672 
1673  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1674  }
1675 
1676  SetUpCoupledMatrix(0.0, u_N, true);
1677  SolveLinearNS(RHS);
1678 
1679  for(int i = 0; i < m_velocity.num_elements(); ++i)
1680  {
1681  u_star[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1682  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_star[i]);
1683 
1684  //u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1685  Vmath::Smul(u_star[i].num_elements(), m_kinvis, u_star[i], 1, u_star[i], 1);
1686  Vmath::Vadd(u_star[i].num_elements(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1687 
1688  m_fields[m_velocity[i]]->FwdTrans(u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1689  }
1690 
1691  m_kinvis -= m_kinvis*m_KinvisPercentage/100;
1692  }
1693 
1694 
1695  void CoupledLinearNS::InfNorm(Array<OneD, Array<OneD, NekDouble> > &inarray,
1696  Array<OneD, NekDouble> &outarray)
1697  {
1698  for(int i = 0; i < m_velocity.num_elements(); ++i)
1699  {
1700  outarray[i] = 0.0;
1701  for(int j = 0; j < inarray[i].num_elements(); ++j)
1702  {
1703  if(inarray[i][j] > outarray[i])
1704  {
1705  outarray[i] = inarray[i][j];
1706  }
1707  }
1708  cout << "InfNorm["<<i<<"] = "<< outarray[i] <<endl;
1709  }
1710  }
1711 
1712  void CoupledLinearNS::L2Norm(Array<OneD, Array<OneD, NekDouble> > &inarray,
1713  Array<OneD, NekDouble> &outarray)
1714  {
1715  for(int i = 0; i < m_velocity.num_elements(); ++i)
1716  {
1717  outarray[i] = 0.0;
1718  for(int j = 0; j < inarray[i].num_elements(); ++j)
1719  {
1720  outarray[i] += inarray[i][j]*inarray[i][j];
1721  }
1722  outarray[i]=sqrt(outarray[i]);
1723  cout << "L2Norm["<<i<<"] = "<< outarray[i] <<endl;
1724  }
1725  }
1726 
1727 
1728  void CoupledLinearNS::EvaluateNewtonRHS(Array<OneD, Array<OneD, NekDouble> > &Velocity,
1729  Array<OneD, Array<OneD, NekDouble> > &outarray)
1730  {
1731  Array<OneD, Array<OneD, NekDouble> > Eval_Adv(m_velocity.num_elements());
1732  Array<OneD, Array<OneD, NekDouble> > tmp_DerVel(m_velocity.num_elements());
1733  Array<OneD, Array<OneD, NekDouble> > AdvTerm(m_velocity.num_elements());
1734  Array<OneD, Array<OneD, NekDouble> > ViscTerm(m_velocity.num_elements());
1735  Array<OneD, Array<OneD, NekDouble> > Forc(m_velocity.num_elements());
1736 
1737  for(int i = 0; i < m_velocity.num_elements(); ++i)
1738  {
1739  Eval_Adv[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1740  tmp_DerVel[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1741 
1742  AdvTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1743  ViscTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1744  Forc[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1745  outarray[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1746 
1747  m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1748 
1749  Vmath::Smul(tmp_DerVel[i].num_elements(), m_kinvis, tmp_DerVel[i], 1, tmp_DerVel[i], 1);
1750  }
1751 
1752  EvaluateAdvectionTerms(Velocity, Eval_Adv);
1753 
1754  for(int i = 0; i < m_velocity.num_elements(); ++i)
1755  {
1756  m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i], AdvTerm[i]); //(w, (u.grad)u)
1757  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1758  m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i], Forc[i]); //(w, f)
1759 
1760  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, AdvTerm[i], 1, outarray[i], 1);
1761  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, ViscTerm[i], 1, outarray[i], 1);
1762 
1763  Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1, Forc[i], 1, outarray[i], 1);
1764  }
1765  }
1766 
1767 
1768 
1769  const SpatialDomains::ExpansionMap &CoupledLinearNS::GenPressureExp(const SpatialDomains::ExpansionMap &VelExp)
1770  {
1771  int i;
1772  SpatialDomains::ExpansionMapShPtr returnval;
1773 
1774  returnval = MemoryManager<SpatialDomains::ExpansionMap>::AllocateSharedPtr();
1775 
1776  SpatialDomains::ExpansionMap::const_iterator expMapIter;
1777  int nummodes;
1778 
1779  for (expMapIter = VelExp.begin(); expMapIter != VelExp.end(); ++expMapIter)
1780  {
1781  LibUtilities::BasisKeyVector BasisVec;
1782 
1783  for(i = 0; i < expMapIter->second->m_basisKeyVector.size(); ++i)
1784  {
1785  LibUtilities::BasisKey B = expMapIter->second->m_basisKeyVector[i];
1786  nummodes = B.GetNumModes();
1787  ASSERTL0(nummodes > 3,"Velocity polynomial space not sufficiently high (>= 4)");
1788  // Should probably set to be an orthogonal basis.
1789  LibUtilities::BasisKey newB(B.GetBasisType(),nummodes-2,B.GetPointsKey());
1790  BasisVec.push_back(newB);
1791  }
1792 
1793  // Put new expansion into list.
1794  SpatialDomains::ExpansionShPtr expansionElementShPtr =
1795  MemoryManager<SpatialDomains::Expansion>::AllocateSharedPtr(expMapIter->second->m_geomShPtr, BasisVec);
1796  (*returnval)[expMapIter->first] = expansionElementShPtr;
1797  }
1798 
1799  // Save expansion into graph.
1800  m_graph->SetExpansions("p",returnval);
1801 
1802  return *returnval;
1803  }
1804 
1805  /**
1806  * @param forcing A list of forcing functions for each velocity
1807  * component
1808  *
1809  * The routine involves two levels of static
1810  * condensations. Initially we require a statically condensed
1811  * forcing function which requires the following manipulation
1812  *
1813  * \f[ {F\_bnd} = {\bf f}_{bnd} -m\_B \,m\_Cinv\, {\bf f}_{int},
1814  * \hspace{1cm} F\_p = m\_D\_{int}\, m\_Cinv\, {\bf f}_{int} \f]
1815  *
1816  * Where \f${\bf f}_{bnd}\f$ denote the forcing degrees of
1817  * freedom of the elemental velocities on the boundary of the
1818  * element, \f${\bf f}_{int}\f$ denote the forcing degrees of
1819  * freedom of the elemental velocities on the interior of the
1820  * element. (see detailed description for more details).
1821  *
1822  * This vector is further manipulated into
1823  *
1824  * \f[ Fh\_{bnd} = \left [ \begin{array}{c} f\_{bnd} -m\_B \,
1825  * m\_Cinv\, {\bf f}_{int}\\ \left [m\_D\_{int} \, m\_Cinv \,{\bf
1826  * f}_{int} \right]_0 \end{array}\right ] \hspace{1cm} [Fh\_p]_{i} =
1827  * \begin{array}{c} [m\_D\_{int} \, m\_Cinv \, {\bf
1828  * f}_{int}]_{i+1} \end{array} \f]
1829  *
1830  * where \f$-{[}m\_D\_{int}^T\, m\_Cinv \,{\bf f}_{int}]_0\f$
1831  * which is corresponds to the mean mode of the pressure degrees
1832  * of freedom and is now added to the boundary system and the
1833  * remainder of the block becomes the interior forcing for the
1834  * inner static condensation (see detailed description for more
1835  * details) which is set up in a GlobalLinSysDirectStaticCond
1836  * class.
1837  *
1838  * Finally we perform the final maniplation of the forcing to
1839  * using hte
1840  * \f[ Fh\_{bnd} = Fh\_{bnd} - m\_Bh \,m\_Chinv \, Fh\_p \f]
1841  *
1842  * We can now call the solver to the global coupled boundary
1843  * system (through the call to #m_CoupledBndSys->Solve) to obtain
1844  * the velocity boundary solution as the mean pressure solution,
1845  * i.e.
1846  *
1847  * \f[ {\cal A}^T(\hat{A} - \hat{C}^T \hat{D}^{-1} \hat{B} ){\cal
1848  * A} \, Bnd = Fh\_{bnd} \f]
1849  *
1850  * Once we know the solution to the above the rest of the pressure
1851  * modes are recoverable thorugh
1852  *
1853  * \f[ Ph = m\_Dhinv\, (Bnd - m\_Ch^T \, Fh_{bnd}) \f]
1854  *
1855  * We can now unpack \f$ Fh\_{bnd} \f$ (last elemental mode) and
1856  * \f$ Ph \f$ into #m_pressure and \f$ F_p\f$ and \f$ Fh\_{bnd}\f$
1857  * into a closed pack list of boundary velocoity degrees of
1858  * freedom stored in \f$ F\_bnd\f$.
1859  *
1860  * Finally the interior velocity degrees of freedom are then
1861  * obtained through the relationship
1862  *
1863  * \f[ F\_{int} = m\_Cinv\ ( F\_{int} + m\_D\_{int}^T\,
1864  * F\_p - m\_Btilde^T\, Bnd) \f]
1865  *
1866  * We then unpack the solution back to the MultiRegion structures
1867  * #m_velocity and #m_pressure
1868  */
1869  void CoupledLinearNS::SolveLinearNS(const Array<OneD, Array<OneD, NekDouble> > &forcing)
1870  {
1871  int i,n;
1872  Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.num_elements());
1873  Array<OneD, Array<OneD, NekDouble> > force(m_velocity.num_elements());
1874 
1875  if(m_HomogeneousType == eHomogeneous1D)
1876  {
1877  int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
1878  for(n = 0; n < m_npointsZ/2; ++n)
1879  {
1880  // Get the a Fourier mode of velocity and forcing components.
1881  for(i = 0; i < m_velocity.num_elements(); ++i)
1882  {
1883  vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2*n);
1884  // Note this needs to correlate with how we pass forcing
1885  force[i] = forcing[i] + 2*n*ncoeffsplane;
1886  }
1887 
1888  SolveLinearNS(force,vel_fields,m_pressure->GetPlane(2*n),n);
1889  }
1890  for(i = 0; i < m_velocity.num_elements(); ++i)
1891  {
1892  m_fields[m_velocity[i]]->SetPhysState(false);
1893  }
1894  m_pressure->SetPhysState(false);
1895  }
1896  else
1897  {
1898  for(i = 0; i < m_velocity.num_elements(); ++i)
1899  {
1900  vel_fields[i] = m_fields[m_velocity[i]];
1901  // Note this needs to correlate with how we pass forcing
1902  force[i] = forcing[i];
1903  }
1904  SolveLinearNS(force,vel_fields,m_pressure);
1905  }
1906  }
1907 
1908  void CoupledLinearNS::SolveLinearNS(const Array<OneD, Array<OneD, NekDouble> > &forcing, Array<OneD, MultiRegions::ExpListSharedPtr> &fields, MultiRegions::ExpListSharedPtr &pressure, const int mode)
1909  {
1910  int i,j,k,n,eid,cnt,cnt1;
1911  int nbnd,nint,offset;
1912  int nvel = m_velocity.num_elements();
1913  int nel = fields[0]->GetNumElmts();
1914  Array<OneD, unsigned int> bmap, imap;
1915 
1916  Array<OneD, NekDouble > f_bnd(m_mat[mode].m_BCinv->GetRows());
1917  NekVector< NekDouble > F_bnd(f_bnd.num_elements(), f_bnd, eWrapper);
1918  Array<OneD, NekDouble > f_int(m_mat[mode].m_BCinv->GetColumns());
1919  NekVector< NekDouble > F_int(f_int.num_elements(),f_int, eWrapper);
1920 
1921  int nz_loc;
1922  int nplanecoeffs = fields[m_velocity[0]]->GetNcoeffs();// this is fine since we pass the nplane coeff data.
1923 
1924  if(mode) // Homogeneous mode flag
1925  {
1926  nz_loc = 2;
1927  }
1928  else
1929  {
1930  if(m_singleMode)
1931  {
1932  nz_loc = 2;
1933  }
1934  else
1935  {
1936  nz_loc = 1;
1937  if(m_HomogeneousType == eHomogeneous1D)
1938  {
1939  // Zero fields to set complex mode to zero;
1940  for(i = 0; i < fields.num_elements(); ++i)
1941  {
1942  Vmath::Zero(2*fields[i]->GetNcoeffs(),fields[i]->UpdateCoeffs(),1);
1943  }
1944  Vmath::Zero(2*pressure->GetNcoeffs(),pressure->UpdateCoeffs(),1);
1945  }
1946  }
1947  }
1948 
1949  // Assemble f_bnd and f_int
1950  cnt = cnt1 = 0;
1951  for(i = 0; i < nel; ++i) // loop over elements
1952  {
1953  eid = fields[m_velocity[0]]->GetOffset_Elmt_Id(i);
1954  fields[m_velocity[0]]->GetExp(eid)->GetBoundaryMap(bmap);
1955  fields[m_velocity[0]]->GetExp(eid)->GetInteriorMap(imap);
1956  nbnd = bmap.num_elements();
1957  nint = imap.num_elements();
1958  offset = fields[m_velocity[0]]->GetCoeff_Offset(eid);
1959 
1960  for(j = 0; j < nvel; ++j) // loop over velocity fields
1961  {
1962  for(n = 0; n < nz_loc; ++n)
1963  {
1964  for(k = 0; k < nbnd; ++k)
1965  {
1966  f_bnd[cnt+k] = forcing[j][n*nplanecoeffs +
1967  offset+bmap[k]];
1968  }
1969  for(k = 0; k < nint; ++k)
1970  {
1971  f_int[cnt1+k] = forcing[j][n*nplanecoeffs +
1972  offset+imap[k]];
1973  }
1974  cnt += nbnd;
1975  cnt1 += nint;
1976  }
1977  }
1978  }
1979 
1980  Array<OneD, NekDouble > f_p(m_mat[mode].m_D_int->GetRows());
1981  NekVector< NekDouble > F_p(f_p.num_elements(),f_p,eWrapper);
1982  NekVector< NekDouble > F_p_tmp(m_mat[mode].m_Cinv->GetRows());
1983 
1984  // fbnd does not currently hold the pressure mean
1985  F_bnd = F_bnd - (*m_mat[mode].m_BCinv)*F_int;
1986  F_p_tmp = (*m_mat[mode].m_Cinv)*F_int;
1987  F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
1988 
1989  // construct inner forcing
1990  Array<OneD, NekDouble > bnd (m_locToGloMap[mode]->GetNumGlobalCoeffs(),0.0);
1991  Array<OneD, NekDouble > fh_bnd(m_locToGloMap[mode]->GetNumGlobalCoeffs(),0.0);
1992 
1993  const Array<OneD,const int>& loctoglomap
1994  = m_locToGloMap[mode]->GetLocalToGlobalMap();
1995  const Array<OneD,const NekDouble>& loctoglosign
1996  = m_locToGloMap[mode]->GetLocalToGlobalSign();
1997 
1998  offset = cnt = 0;
1999  for(i = 0; i < nel; ++i)
2000  {
2001  eid = fields[0]->GetOffset_Elmt_Id(i);
2002  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2003 
2004  for(j = 0; j < nvel; ++j)
2005  {
2006  for(k = 0; k < nbnd; ++k)
2007  {
2008  fh_bnd[loctoglomap[offset+j*nbnd+k]] +=
2009  loctoglosign[offset+j*nbnd+k]*f_bnd[cnt+k];
2010  }
2011  cnt += nbnd;
2012  }
2013 
2014  nint = pressure->GetExp(eid)->GetNcoeffs();
2015  offset += nvel*nbnd + nint*nz_loc;
2016  }
2017 
2018  offset = cnt1 = 0;
2019  for(i = 0; i < nel; ++i)
2020  {
2021  eid = fields[0]->GetOffset_Elmt_Id(i);
2022  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2023  nint = pressure->GetExp(eid)->GetNcoeffs();
2024 
2025  for(n = 0; n < nz_loc; ++n)
2026  {
2027  for(j = 0; j < nint; ++j)
2028  {
2029  fh_bnd[loctoglomap[offset + nvel*nbnd + n*nint+j]] = f_p[cnt1+j];
2030  }
2031  cnt1 += nint;
2032  }
2033  offset += nvel*nbnd + nz_loc*nint;
2034  }
2035 
2036  // Set Weak BC into f_bnd and Dirichlet Dofs in bnd
2037  const Array<OneD,const int>& bndmap
2038  = m_locToGloMap[mode]->GetBndCondCoeffsToGlobalCoeffsMap();
2039 
2040  // Forcing function with weak boundary conditions and
2041  // Dirichlet conditions
2042  int bndcnt=0;
2043 
2044  for(k = 0; k < nvel; ++k)
2045  {
2046  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConds = fields[k]->GetBndConditions();
2047  Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp;
2048  if(m_HomogeneousType == eHomogeneous1D)
2049  {
2050  bndCondExp = m_fields[k]->GetPlane(2*mode)->GetBndCondExpansions();
2051  }
2052  else
2053  {
2054  bndCondExp = m_fields[k]->GetBndCondExpansions();
2055  }
2056 
2057  for(i = 0; i < bndCondExp.num_elements(); ++i)
2058  {
2059  const Array<OneD, const NekDouble > bndCondCoeffs = bndCondExp[i]->GetCoeffs();
2060  cnt = 0;
2061  for(n = 0; n < nz_loc; ++n)
2062  {
2063  if(bndConds[i]->GetBoundaryConditionType()
2064  == SpatialDomains::eDirichlet)
2065  {
2066  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
2067  {
2068  if (m_equationType == eSteadyNavierStokes && m_initialStep == false)
2069  {
2070  //This condition set all the Dirichlet BC at 0 after
2071  //the initial step of the Newton method
2072  bnd[bndmap[bndcnt++]] = 0;
2073  }
2074  else
2075  {
2076  bnd[bndmap[bndcnt++]] = bndCondCoeffs[cnt++];
2077  }
2078  }
2079  }
2080  else
2081  {
2082  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
2083  {
2084  fh_bnd[bndmap[bndcnt++]]
2085  += bndCondCoeffs[cnt++];
2086  }
2087  }
2088  }
2089  }
2090  }
2091 
2092  m_mat[mode].m_CoupledBndSys->Solve(fh_bnd,bnd,m_locToGloMap[mode]);
2093 
2094  // unpack pressure and velocity boundary systems.
2095  offset = cnt = 0;
2096  int totpcoeffs = pressure->GetNcoeffs();
2097  Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2098  for(i = 0; i < nel; ++i)
2099  {
2100  eid = fields[0]->GetOffset_Elmt_Id(i);
2101  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2102  nint = pressure->GetExp(eid)->GetNcoeffs();
2103 
2104  for(j = 0; j < nvel; ++j)
2105  {
2106  for(k = 0; k < nbnd; ++k)
2107  {
2108  f_bnd[cnt+k] = loctoglosign[offset+j*nbnd+k]*bnd[loctoglomap[offset + j*nbnd + k]];
2109  }
2110  cnt += nbnd;
2111  }
2112  offset += nvel*nbnd + nint*nz_loc;
2113  }
2114 
2115  pressure->SetPhysState(false);
2116 
2117  offset = cnt = cnt1 = 0;
2118  for(i = 0; i < nel; ++i)
2119  {
2120  eid = fields[0]->GetOffset_Elmt_Id(i);
2121  nint = pressure->GetExp(eid)->GetNcoeffs();
2122  nbnd = fields[0]->GetExp(eid)->NumBndryCoeffs();
2123  cnt1 = pressure->GetCoeff_Offset(eid);
2124 
2125  for(n = 0; n < nz_loc; ++n)
2126  {
2127  for(j = 0; j < nint; ++j)
2128  {
2129  p_coeffs[n*totpcoeffs + cnt1+j] =
2130  f_p[cnt+j] = bnd[loctoglomap[offset +
2131  (nvel*nz_loc)*nbnd +
2132  n*nint + j]];
2133  }
2134  cnt += nint;
2135  }
2136  offset += (nvel*nbnd + nint)*nz_loc;
2137  }
2138 
2139  // Back solve first level of static condensation for interior
2140  // velocity space and store in F_int
2141  F_int = F_int + Transpose(*m_mat[mode].m_D_int)*F_p
2142  - Transpose(*m_mat[mode].m_Btilde)*F_bnd;
2143  F_int = (*m_mat[mode].m_Cinv)*F_int;
2144 
2145  // Unpack solution from Bnd and F_int to v_coeffs
2146  cnt = cnt1 = 0;
2147  for(i = 0; i < nel; ++i) // loop over elements
2148  {
2149  eid = fields[m_velocity[0]]->GetOffset_Elmt_Id(i);
2150  fields[0]->GetExp(eid)->GetBoundaryMap(bmap);
2151  fields[0]->GetExp(eid)->GetInteriorMap(imap);
2152  nbnd = bmap.num_elements();
2153  nint = imap.num_elements();
2154  offset = fields[0]->GetCoeff_Offset(eid);
2155 
2156  for(j = 0; j < nvel; ++j) // loop over velocity fields
2157  {
2158  for(n = 0; n < nz_loc; ++n)
2159  {
2160  for(k = 0; k < nbnd; ++k)
2161  {
2162  fields[j]->SetCoeff(n*nplanecoeffs +
2163  offset+bmap[k],
2164  f_bnd[cnt+k]);
2165  }
2166 
2167  for(k = 0; k < nint; ++k)
2168  {
2169  fields[j]->SetCoeff(n*nplanecoeffs +
2170  offset+imap[k],
2171  f_int[cnt1+k]);
2172  }
2173  cnt += nbnd;
2174  cnt1 += nint;
2175  }
2176  }
2177  }
2178 
2179  for(j = 0; j < nvel; ++j)
2180  {
2181  fields[j]->SetPhysState(false);
2182  }
2183  }
2184 
2185  void CoupledLinearNS::v_Output(void)
2186  {
2187  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_fields.num_elements()+1);
2188  std::vector<std::string> variables(m_fields.num_elements()+1);
2189  int i;
2190 
2191  for(i = 0; i < m_fields.num_elements(); ++i)
2192  {
2193  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2194  variables[i] = m_boundaryConditions->GetVariable(i);
2195  }
2196 
2197  fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2198  // project pressure field to velocity space
2199  if(m_singleMode==true)
2200  {
2201  Array<OneD, NekDouble > tmpfieldcoeffs (m_fields[0]->GetNcoeffs()/2);
2202  m_pressure->GetPlane(0)->BwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetCoeffs(), m_pressure->GetPlane(0)->UpdatePhys());
2203  m_pressure->GetPlane(1)->BwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetCoeffs(), m_pressure->GetPlane(1)->UpdatePhys());
2204  m_fields[0]->GetPlane(0)->FwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetPhys(),fieldcoeffs[i]);
2205  m_fields[0]->GetPlane(1)->FwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetPhys(),tmpfieldcoeffs);
2206  for(int e=0; e<m_fields[0]->GetNcoeffs()/2; e++)
2207  {
2208  fieldcoeffs[i][e+m_fields[0]->GetNcoeffs()/2] = tmpfieldcoeffs[e];
2209  }
2210  }
2211  else
2212  {
2213  m_pressure->BwdTrans_IterPerExp(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
2214  m_fields[0]->FwdTrans_IterPerExp(m_pressure->GetPhys(),fieldcoeffs[i]);
2215  }
2216  variables[i] = "p";
2217 
2218  std::string outname = m_sessionName + ".fld";
2219 
2220  WriteFld(outname,m_fields[0],fieldcoeffs,variables);
2221  }
2222 
2223  int CoupledLinearNS::v_GetForceDimension()
2224  {
2225  return m_session->GetVariables().size();
2226  }
2227 }
2228 
2229 /**
2230  * $Log: CoupledLinearNS.cpp,v $
2231  **/
Generated on Wed Mar 4 2015 16:33:13 for Nektar++ by   doxygen 1.8.1.2