Nektar++
GlobalLinSysXxtFull.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File GlobalLinSysXxtFull.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: GlobalLinSysXxtFull definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 #include <MultiRegions/ExpList.h>
39 
40 namespace Nektar
41 {
42  namespace MultiRegions
43  {
44  /**
45  * @class GlobalLinSysXxtFull
46  */
47 
48  /**
49  * Registers the class with the Factory.
50  */
53  "XxtFull",
55  "Xxt Full Matrix.");
56 
57 
58  /// Constructor for full direct matrix solve.
60  const GlobalLinSysKey &pLinSysKey,
61  const boost::weak_ptr<ExpList> &pExp,
62  const boost::shared_ptr<AssemblyMap>
63  &pLocToGloMap)
64  : GlobalLinSys (pLinSysKey, pExp, pLocToGloMap),
65  GlobalLinSysXxt(pLinSysKey, pExp, pLocToGloMap)
66  {
67 
69  "This routine should only be used when using a Full XXT"
70  " matrix solve");
71 
72  CreateMap(pLocToGloMap);
73  AssembleMatrixArrays(pLocToGloMap);
74  }
75 
76 
78  {
79 
80  }
81 
82 
83  /**
84  * Solve the linear system using a full global matrix system.
85  */
87  const Array<OneD, const NekDouble> &pInput,
88  Array<OneD, NekDouble> &pOutput,
89  const AssemblyMapSharedPtr &pLocToGloMap,
90  const Array<OneD, const NekDouble> &pDirForcing)
91  {
92  bool dirForcCalculated = (bool) pDirForcing.num_elements();
93  int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
94  int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
95 
96  Array<OneD, NekDouble> tmp (nGlobDofs);
97  Array<OneD, NekDouble> tmp2(nGlobDofs);
98  Array<OneD, NekDouble> tmp3 = pOutput + nDirDofs;
99 
100  if(nDirDofs)
101  {
102  // calculate the dirichlet forcing
103  if(dirForcCalculated)
104  {
105  Vmath::Vsub(nGlobDofs, pInput.get(), 1,
106  pDirForcing.get(), 1,
107  tmp.get(), 1);
108  }
109  else
110  {
111  // Calculate the dirichlet forcing and substract it
112  // from the rhs
113  //int nLocDofs = pLocToGloMap->GetNumLocalCoeffs();
114 
115  m_expList.lock()->GeneralMatrixOp(
116  m_linSysKey,
117  pOutput, tmp, eGlobal);
118 
119  Vmath::Vsub( nGlobDofs, pInput.get(),1,
120  tmp.get(), 1,
121  tmp.get(), 1);
122  }
123 
124  SolveLinearSystem(pLocToGloMap->GetNumLocalCoeffs(),
125  tmp, tmp2, pLocToGloMap);
126 
127  // Enforce the Dirichlet boundary conditions on the solution
128  // array as XXT discards them.
129  Vmath::Vcopy(nDirDofs, pOutput, 1,
130  tmp2, 1);
131  }
132  else
133  {
134  Vmath::Vcopy(nGlobDofs, pInput, 1, tmp, 1);
135  SolveLinearSystem(pLocToGloMap->GetNumLocalCoeffs(),
136  tmp,tmp2, pLocToGloMap);
137  }
138 
139  // Perturb the output array (previous solution) by the result of
140  // this solve to get full solution.
141  Vmath::Vadd(nGlobDofs - nDirDofs,
142  tmp2 + nDirDofs, 1, tmp3, 1, tmp3, 1);
143 
144  }
145 
146 
147  /**
148  * Create the inverse multiplicity map.
149  * @param locToGloMap Local to global mapping information.
150  */
152  const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
153  {
154  const Array<OneD, const int> &vMap
155  = pLocToGloMap->GetLocalToGlobalMap();
156  unsigned int nGlo = pLocToGloMap->GetNumGlobalCoeffs();
157  unsigned int nEntries = pLocToGloMap->GetNumLocalCoeffs();
158  unsigned int i;
159 
160  // Count the multiplicity of each global DOF on this process
161  Array<OneD, NekDouble> vCounts(nGlo, 0.0);
162  for (i = 0; i < nEntries; ++i)
163  {
164  vCounts[vMap[i]] += 1.0;
165  }
166 
167  // Get universal multiplicity by globally assembling counts
168  pLocToGloMap->UniversalAssemble(vCounts);
169 
170  // Construct a map of 1/multiplicity for use in XXT solve
172  for (i = 0; i < nEntries; ++i)
173  {
174  m_locToGloSignMult[i] = 1.0/vCounts[vMap[i]];
175  }
176 
177  m_map = pLocToGloMap->GetLocalToGlobalMap();
178  }
179 
180 
181  /**
182  * Construct the local matrix row index, column index and value index
183  * arrays and initialize the XXT data structure with this information.
184  * @param locToGloMap Local to global mapping information.
185  */
187  const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
188  {
189  ExpListSharedPtr vExp = m_expList.lock();
190  unsigned int nElmt = vExp->GetNumElmts();
191  DNekScalMatSharedPtr loc_mat;
192  unsigned int iCount = 0;
193  unsigned int rCount = 0;
194  unsigned int nRows = 0;
195  unsigned int nEntries = 0;
196  unsigned int numDirBnd = pLocToGloMap->GetNumGlobalDirBndCoeffs();
197  unsigned int nLocal = pLocToGloMap->GetNumLocalCoeffs();
198  const Array<OneD, NekDouble> &vMapSign
199  = pLocToGloMap->GetLocalToGlobalSign();
200  bool doSign = pLocToGloMap->GetSignChange();
201  unsigned int i = 0, j = 0, k = 0, n = 0;
202  int gid1;
203  Array<OneD, unsigned int> vSizes(nElmt);
204 
205  // First construct a map of the number of local DOFs in each block
206  // and the number of matrix entries for each block
207 
208  // Dimension of matrix is just the linear vertex space
211  {
212  for (n = 0; n < nElmt; ++n)
213  {
214  i = vExp->GetOffset_Elmt_Id(n);
215  vSizes[n] = vExp->GetExp(i)->GetNverts();
216  nEntries += vSizes[n]*vSizes[n];
217  }
218  }
219  else
220  {
221  for (n = 0; n < nElmt; ++n)
222  {
223  i = vExp->GetOffset_Elmt_Id(n);
224  vSizes[n] = vExp->GetExp(i)->GetNcoeffs();
225  nEntries += vSizes[n]*vSizes[n];
226  }
227  }
228 
229  // Set up i-index, j-index and value arrays
230  m_Ai = Array<OneD, unsigned int>(nEntries);
231  m_Aj = Array<OneD, unsigned int>(nEntries);
232  m_Ar = Array<OneD, double>(nEntries, 0.0);
233 
234  // Set up the universal ID array for XXT
235  Array<OneD, unsigned long> vId(nLocal);
236 
237  // Loop over each elemental block, extract matrix indices and value
238  // and set the universal ID array
239  for(n = iCount = 0; n < nElmt; ++n)
240  {
241  loc_mat = GetBlock(vExp->GetOffset_Elmt_Id(n));
242  nRows = loc_mat->GetRows();
243 
244  for(i = 0; i < nRows; ++i)
245  {
246  gid1 = pLocToGloMap->GetLocalToGlobalMap(iCount + i);
247  for(j = 0; j < nRows; ++j)
248  {
249  k = rCount + i*vSizes[n] + j;
250  m_Ai[k] = iCount + i;
251  m_Aj[k] = iCount + j;
252  m_Ar[k] = (*loc_mat)(i,j);
253  if (doSign)
254  {
255  m_Ar[k] *= vMapSign[iCount+i]*vMapSign[iCount+j];
256  }
257  }
258 
259  // Dirichlet DOFs are not included in the solve, so we set
260  // these to the special XXT id=0.
261  if (gid1 < numDirBnd)
262  {
263  vId[iCount + i] = 0;
264  }
265  else
266  {
267  vId[iCount + i]
268  = pLocToGloMap->GetGlobalToUniversalMap(gid1);
269  }
270  }
271  iCount += vSizes[n];
272  rCount += vSizes[n]*vSizes[n];
273  }
274 
275  // Set up XXT and output some stats
276  LibUtilities::CommSharedPtr vComm = pLocToGloMap->GetComm();
277  m_crsData = Xxt::Init(nLocal, vId, m_Ai, m_Aj, m_Ar, vComm);
279  }
280  }
281 }
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
static std::string className
Name of class.
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
void CreateMap(const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
void SolveLinearSystem(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir=0)
Solve the linear system for given input and output vectors.
Definition: GlobalLinSys.h:193
GlobalLinSysXxtFull(const GlobalLinSysKey &pLinSysKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Constructor for full direct matrix solve.
virtual void v_Solve(const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Solve the linear system for given input and output vectors using a specified local to global map...
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Global coefficients.
DNekScalMatSharedPtr GetBlock(unsigned int n)
Definition: GlobalLinSys.h:214
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
static GlobalLinSysSharedPtr create(const GlobalLinSysKey &pLinSysKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
Describe a linear system.
static struct crs_data * Init(unsigned int pRank, const Nektar::Array< OneD, unsigned long > pId, const Nektar::Array< OneD, unsigned int > pAi, const Nektar::Array< OneD, unsigned int > pAj, const Nektar::Array< OneD, NekDouble > pAr, const LibUtilities::CommSharedPtr &pComm)
Initialise the matrix-solve.
Definition: Xxt.hpp:157
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
Definition: GlobalLinSys.h:123
void AssembleMatrixArrays(const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
Array< OneD, unsigned int > m_Aj
A global linear system.
Definition: GlobalLinSys.h:70
Array< OneD, unsigned int > m_Ai
GlobalLinSysFactory & GetGlobalLinSysFactory()
void nektar_crs_stats(struct crs_data *data)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Array< OneD, NekDouble > m_locToGloSignMult
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:125