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