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