HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UT_MultigridArray.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_MultigridArray.h ( UT Library, C++)
7  *
8  * COMMENTS:
9  * This provides an array that can solve the 2D or 3D Poisson equation
10  * using the multigrid algorithm.
11  * There are only explicit template instantations for fpreal32 and
12  * fpreal64
13  *
14  * The created array has elements indexed from 0, ie: [0..xdiv-1].
15  */
16 
17 #ifndef __UT_MultigridArray__
18 #define __UT_MultigridArray__
19 
20 #include "UT_API.h"
21 #include "UT_Vector.h"
22 #include "UT_Vector3.h"
23 #include "UT_ThreadedAlgorithm.h"
24 #include "UT_ValArray.h"
25 #include "UT_VoxelArray.h"
26 #include "UT_Matrix.h"
27 
28 
30 {
33 };
34 
35 ///
36 /// UT_MultigridArrayT
37 /// This provides an array that can solve the 2D or 3D Poisson equation
38 /// using the multigrid algorithm.
39 /// There are currently only explicit template instantations for fpreal32 and
40 /// fpreal64.
41 ///
42 /// The created array has elements indexed from 0, ie: [0..xdiv-1].
43 template<typename T>
45 {
46 public:
47  /// Default constructor, array will not be initialized after this, i.e.
48  /// isInit() == false
50  virtual ~UT_MultigridArrayT();
51 
52  /// Initialize array to given size and spacing
53  /// NB: For the 2d case, UT_MultigridArrays MUST have the shape (nx, 1, ny).
54  void init( const UT_Vector3I &res,
55  const UT_Vector3T<T> &spacing,
56  const UT_Vector3I &boundariesNeg,
57  const UT_Vector3I &boundariesPos,
58  int level = 0,
59  UT_Vector3I oddCoarsenings = UT_Vector3I(0,0,0));
60 
61  /// Initializes this array from the provided UT_VoxelArray, assuming the
62  /// provided grid spacing. After this function returns, isInit() will
63  /// return true.
64  /// If this array has already been initialized, it will be resized as
65  /// needed to match the UT_VoxelArray.
66  /// NB: In the 2d case, UT_MultigridArrays MUST have the shape (nx, 1, ny).
67  /// This function ensures that will be the case, but if a provided 2D
68  /// UT_VoxelArray has a different shape, e.g. (nx, ny, 1), then the shapes
69  /// of the two arrays will not match after this function, although the
70  /// number of voxels will. copyToVoxels will only resize if the number of
71  /// voxels differs, so initFromVoxels and copyToVoxels with the same
72  /// UT_VoxelArray will work fine in the 2D and 3D case.
73  void initFromVoxels(const UT_VoxelArray<T> &f,
74  const UT_Vector3T<T> &spacing,
75  const UT_Vector3I &boundariesNeg,
76  const UT_Vector3I &boundariesPos);
77 
78  /// Copy the array data to the provided UT_VoxelArray. The UT_VoxelArray
79  /// will be resized if the number of voxels doesn't match the number of
80  /// elements in this array.
81  void copyToVoxels(UT_VoxelArray<T> &f) const;
82 
83  /// Initializes this to be the same dimensions, boundary conditions,
84  /// etc, of the given array.
85  /// The values of this may be reset to zero.
86  void match( const UT_MultigridArrayT &src);
87 
88  /// Whether this array has been initialized - default constructor does not.
89  bool isInit()
90  {
91  return myVec.isInit();
92  }
93 
94  // 3D operator access
95  T operator()(int i, int j, int k) const
96  {
97  return myVec(i + j * myYStride + k * myZStride);
98  }
99 
100  T &operator()(int i, int j, int k)
101  {
102  return myVec(i + j * myYStride + k * myZStride);
103  };
104 
105  T operator()(const UT_Vector3I &idx) const
106  {
107  return myVec(idx[0] + idx[1] * myYStride + idx[2] * myZStride);
108  }
109 
110  T &operator()(const UT_Vector3I &idx)
111  {
112  return myVec(idx[0] + idx[1] * myYStride + idx[2] * myZStride);
113  };
114 
115  /// Returns the ghost value for the negative side boundaries, given the
116  /// values at idx and idxadj, which should represent the index for a
117  /// border cell and the one adjacent to it along the provided axis.
118  T ghostValueNeg(int axis, const UT_Vector3I &idx,
119  const UT_Vector3I &idxadj) const
120  {
121  return myBoundModNeg[axis] * (*this)(idx) +
122  myBoundAdjModNeg[axis] * (*this)(idxadj);
123  }
124 
125  /// Returns the ghost value for the positive side boundaries, given the
126  /// values at idx and idxadj, which should represent the index for a
127  /// border cell and the one adjacent to it along the provided axis.
128  T ghostValuePos(int axis, const UT_Vector3I &idx,
129  const UT_Vector3I &idxadj) const
130  {
131  return myBoundModPos[axis] * (*this)(idx) +
132  myBoundAdjModPos[axis] * (*this)(idxadj);
133  }
134 
135  /// Returns the ghost value for the negative side boundaries, given the
136  /// provided values U and Uadj, which should be the values at the border
137  /// cell and the one adjacent to it along the provided axis.
138  T ghostValueNeg(int axis, T U, T Uadj) const
139  {
140  return myBoundModNeg[axis] * U + myBoundAdjModNeg[axis] * Uadj;
141  }
142 
143  /// Returns the ghost value for the positive side boundaries, given the
144  /// provided values U and Uadj, which should be the values at the border
145  /// cell and the one adjacent to it along the provided axis.
146  T ghostValuePos(int axis, T U, T Uadj) const
147  {
148  return myBoundModPos[axis] * U + myBoundAdjModPos[axis] * Uadj;
149  }
150 
151  UT_MultigridArrayT<T> &operator =(const UT_MultigridArrayT &a);
152 
153  UT_MultigridArrayT<T> &operator +=(const UT_MultigridArrayT &a);
154 
155  void zero()
156  {
157  myVec.zero();
158  }
159 
160  T norm(int n) const
161  {
162  return myVec.norm(n);
163  }
164 
165  exint getRes(int axis) const
166  {
167  return myRes[axis];
168  }
169 
171  {
172  return myRes;
173  }
174 
175  const UT_Vector3T<T> &getSpacing() const
176  {
177  return mySpacing;
178  }
179 
181  {
182  return myBoundariesNeg;
183  }
184 
186  {
187  return myBoundariesPos;
188  }
189 
190  int getLevel() const
191  {
192  return myLevel;
193  }
194 
196  {
197  UT_Vector3I parity;
198  for (int i=0; i < 3; i++)
199  parity[i] = getRes(i) % 2;
200  return parity;
201  }
202 
203  /// Returns the total number of elements in the array.
205  {
206  return myRes[0] * myRes[1] * myRes[2];
207  }
208 
209  /// Calculates the maximum coarse grid elements in the coarsest grid, based
210  /// on the minPerAxis parameter.
211  exint getMaxCoarse(exint minPerAxis) const;
212 
213  /// Returns whether this is a 1-D array.
214  bool is1D() const
215  {
216  int ndim = 0;
217  for (int i=0; i < 3; i++)
218  if (getRes(i) > 1)
219  ndim++;
220  return ndim < 2;
221  }
222 
223  bool shouldMultiThread() const
224  {
225  return numElements() > 1000;
226  }
227 
228  const T *data() const
229  {
230  return myVec.getData();
231  }
232 
233  T *data()
234  {
235  return myVec.getData();
236  }
237 
238  /// Computes the 0 and 2-norm of this in one pass
239  void computeNorms(fpreal64 &norminf, fpreal64 &norm2) const;
240 
241  /// Computes the inf and 2-norm of the residual r = b - A * x,
242  /// where b == *this
243  void computeResidualNorms(const UT_MultigridArrayT &x,
244  fpreal64 &norminf, fpreal64 &norm2) const;
245 
246  /// Computes rhe residual r = b - A * x, where b == *this
247  THREADED_METHOD2_CONST(UT_MultigridArrayT, shouldMultiThread(),
248  subtractApplyLaplacian,
249  const UT_MultigridArrayT&, x,
251  void subtractApplyLaplacianPartial(const UT_MultigridArrayT& x,
253  const UT_JobInfo& info) const;
254 
255  /// Performs red-black Gauss-Seidel smoothing according to parity, where b
256  /// == *this, and the implicit matrix is the Laplacian operator.
257  /// This should be called twice, with parity=0 and 1, for a full smoothing
258  /// cycle.
259  THREADED_METHOD2_CONST(UT_MultigridArrayT, shouldMultiThread(),
260  smoothLaplacianGaussSeidel,
261  int, parity,
262  UT_MultigridArrayT&, x)
263  void smoothLaplacianGaussSeidelPartial(int parity, UT_MultigridArrayT& x,
264  const UT_JobInfo& info) const;
265 
266  /// Performs coarsening using full-weighting along one axis.
267  /// On output uh.shape[axis] == this->shape[axis] / 2
268  THREADED_METHOD2_CONST(UT_MultigridArrayT, shouldMultiThread(),
269  coarsenAlongAxis,
270  UT_MultigridArrayT&, uh,
271  int, axis)
272  void coarsenAlongAxisPartial(UT_MultigridArrayT& uh, int axis,
273  const UT_JobInfo& info) const;
274 
275  /// Interpolate this array into u using inverse of full-weighting along one
276  /// axis.
277  /// On output u.shape[axis] == this->shape[axis] * 2 + parity.
278  THREADED_METHOD2_CONST(UT_MultigridArrayT, shouldMultiThread(),
279  interpolateAlongAxis,
280  UT_MultigridArrayT&, u,
281  int, axis)
282  void interpolateAlongAxisPartial(UT_MultigridArrayT& u, int axis,
283  const UT_JobInfo& info) const;
284 
285 
286  /// Build the dense matrix in A that corresponds to this grid's Laplacian
287  /// operator.
288  template <typename S>
289  void buildMatrix(UT_MatrixT<S> &A) const;
290 
291  /// Directly solve the Poisson equation that this array represents, using
292  /// either Cholesky, LU, or SVD decomposition, depending on the nature of
293  /// the Poisson problem. Should only be called for very coarse grids.
294  void directSolve(UT_MultigridArrayT &x) const;
295 
296  /// Coarsen this array into a lower resolution version of the grid.
297  /// coarsenAxis is a vector of booleans specifying which axes to coarsen.
298  /// Where coarsenAxis[axis] is true, ensure that uh[axis] = getRes(axis) / 2
299  UT_Vector3I coarsen(exint minPerAxis, UT_MultigridArrayT &uh) const;
300 
301  /// Interpolate this array into a fine version of the grid.
302  /// interpolateAxis is a vector of booleans specifying which axes to
303  /// interpolate
304  /// Where interpolateAxis[axis] is true, ensure that
305  /// u[axis] = getRes(axis) * 2 + parity[axis].
306  void interpolate(const UT_Vector3I &coarsenAxis, const UT_Vector3I &parity,
307  UT_MultigridArrayT &u) const;
308 
309  /// Perform recursive V-cycle for multigrid algorithm, until reaching
310  /// numElements() < maxcoarse,
311  /// at which point a direct solve is done using directSolve.. At
312  /// each level, nSmoothDown and nSmoothUp iterations of Gauss-Seidel
313  /// smoothing are done, unless smoothTopLevelDown is false, in which
314  /// case no smoothing is done for the top-level grid on the down-cycle. This
315  /// can be used as an optimization when calling vcycle iteratively; the
316  /// previous vcycle call will have smoothed the top-level grid as its last
317  /// step, so you can often skip smoothing that same grid on the next vcycle
318  /// without loss of convergence.
319  void vcycle(exint minPerAxis, int nSmoothDown, int nSmoothUp,
320  UT_MultigridArrayT &x, bool smoothTopLevelDown = true) const;
321 
322  /// Performs recursive Full Multigrid to generate an initial guess for the
323  /// solution of the Poisson equation represented by this array. At
324  /// each level, nSmoothDown and nSmoothUp iterations of smoothing are done.
325  void fullMultigrid(exint minPerAxis, int nSmoothDown,
326  int nSmoothUp, UT_MultigridArrayT &x) const;
327 
328  /// Iteratively solve the Poisson equation, assuming *this is the right
329  /// hand side. This function will do at least miniter iterations, then
330  /// continue to iterate until the 2-norm of the residual has
331  /// been reduced by reltol, the inf-norm of the residual is below abstol,
332  /// or maxiters is reached.
333  /// The optional resnorm0 and resnorm2 parameters will contain the
334  /// 0- and 2-norms of the residuals for each iteration.
335  int solvePoisson(fpreal64 abstol, fpreal64 reltol, int miniter, int maxiter,
336  UT_MultigridArrayT &x,
337  UT_ValArray<fpreal64> *resnorminf,
338  UT_ValArray<fpreal64> *resnorm2) const;
339 
340 protected:
341  /// Directly solves the Poisson equation for this right-hand side, storing
342  /// the result in x, using UT_MatrixSolver::choleskyDecomp and
343  /// UT_MatrixSolver::choleskyBackSub.
344  /// This should only be called for relatively small N when the coarsest
345  /// grid level has been reached, and only when the result of buildMatrix
346  /// is positive (semi-definite)
347  bool solvePoissonCholesky(UT_MultigridArrayT &x) const;
348 
349  /// Directly solves the Poisson equation for this right-hand side, using
350  /// LU decomposition.
351  /// This should only be called for relatively small N when the coarsest
352  /// grid level has been reached, and only when the result of buildMatrix
353  /// is non-singular.
354  bool solvePoissonLU(UT_MultigridArrayT &x) const;
355 
356  /// Directly solves the Poisson equation for this right-hand side, using
357  /// SVD decomposition. This can be used when the result of buildMatrix is
358  /// singular. This should only be called for relatively small N when the
359  /// coarsest grid level has been reached.
360  bool solvePoissonSVD(UT_MultigridArrayT &x) const;
361 
362  void initStorage();
363  void initLaplacian();
364 
365  THREADED_METHOD2_CONST(UT_MultigridArrayT, shouldMultiThread(),
366  computeNormsInternal,
367  fpreal64*, norminf,
368  fpreal64*, norm2)
369  void computeNormsInternalPartial(fpreal64 *norminf, fpreal64 *norm2,
370  const UT_JobInfo& info) const;
371 
372  THREADED_METHOD3_CONST(UT_MultigridArrayT, shouldMultiThread(),
373  computeResidualNormsInternal,
374  const UT_MultigridArrayT&, x,
375  fpreal64*, norminf,
376  fpreal64*, norm2)
377  void computeResidualNormsInternalPartial(const UT_MultigridArrayT& x,
378  fpreal64 *norminf, fpreal64 *norm2,
379  const UT_JobInfo& info) const;
380 
381  /// Since the multi-threading threshold is 1000,
382  /// pick a block size less than that.
383  static const exint PARALLEL_BLOCK_SIZE_VEC = 512;
384 
385  /// This is for parallelism over just z, so a z section of size just 8
386  /// can be huge if x and y are large.
387  static const exint PARALLEL_BLOCK_SIZE_Z = 8;
388 
389  /// UT_VectorT that holds the array data.
390  UT_VectorT<T> myVec;
391 
392  /// Number of elements in each dimension.
393  UT_Vector3I myRes;
394 
395  /// Grid spacing in each dimension.
396  UT_Vector3T<T> mySpacing;
397 
398  /// Array stride in y and z axes.
399  exint myYStride, myZStride;
400 
401  /// Laplacian operator values calculated based on grid spacing.
402  T myInvdx2, myInvdy2, myInvdz2, myDiag, myOmega;
403  UT_Vector3I myBoundariesNeg, myBoundariesPos;
404  UT_Vector3T<T> myBoundModNeg, myBoundModPos;
405  UT_Vector3T<T> myBoundAdjModNeg, myBoundAdjModPos;
406  UT_Vector3T<T> myDiagModNeg, myDiagModPos;
407  int myLevel;
408  bool myAllOpen, myAllClosed;
409  UT_Vector3I myOddCoarsenings;
410 };
411 
412 typedef UT_MultigridArrayT<fpreal> UT_MultigridArrayR;
413 typedef UT_MultigridArrayT<fpreal32> UT_MultigridArrayF;
414 typedef UT_MultigridArrayT<fpreal64> UT_MultigridArrayD;
415 
416 #endif
T operator()(const UT_Vector3I &idx) const
UT_Vector3I getRes() const
bool shouldMultiThread() const
const UT_Vector3I & getBoundariesPos() const
T ghostValueNeg(int axis, const UT_Vector3I &idx, const UT_Vector3I &idxadj) const
exint numElements() const
Returns the total number of elements in the array.
exint getRes(int axis) const
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLint level
Definition: glcorearb.h:107
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
UT_Vector3I getParity() const
#define UT_API
Definition: UT_API.h:12
UT_Vector3T< int64 > UT_Vector3I
bool isInit()
Whether this array has been initialized - default constructor does not.
png_uint_32 i
Definition: png.h:2877
GLdouble n
Definition: glcorearb.h:2007
const UT_Vector3I & getBoundariesNeg() const
GLfloat f
Definition: glcorearb.h:1925
T & operator()(const UT_Vector3I &idx)
int64 exint
Definition: SYS_Types.h:109
double fpreal64
Definition: SYS_Types.h:185
T & operator()(int i, int j, int k)
T norm(int n) const
const T * data() const
T ghostValueNeg(int axis, T U, T Uadj) const
const UT_Vector3T< T > & getSpacing() const
#define THREADED_METHOD2_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2)
double fpreal
Definition: SYS_Types.h:263
bool is1D() const
Returns whether this is a 1-D array.
GLint GLenum GLint x
Definition: glcorearb.h:408
UT_BoundaryCondition
GLboolean r
Definition: glcorearb.h:1221
#define const
Definition: zconf.h:214
T operator()(int i, int j, int k) const
#define THREADED_METHOD3_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3)
float fpreal32
Definition: SYS_Types.h:184
T ghostValuePos(int axis, const UT_Vector3I &idx, const UT_Vector3I &idxadj) const
T ghostValuePos(int axis, T U, T Uadj) const
GLenum src
Definition: glcorearb.h:1792