HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CE_Multigrid.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: CE_Multigrid.h ( CE Library, C++)
7  *
8  * COMMENTS: Compute Engine Multigrid solver
9  */
10 
11 #ifndef __CE_Multigrid__
12 #define __CE_Multigrid__
13 
14 #include "CE_Grid.h"
15 
16 #ifdef CE_ENABLED
17 
18 #include <UT/UT_Vector3.h>
19 
20 /// This class provides a specializaion of CE_Grid that can solve the 2D or
21 /// 3D Poisson equation using the multigrid algorithm, running on an OpenCL
22 /// device. It's interface closely matches that of UT_MultigridArray.
24 {
25 public:
26  CE_Multigrid();
27  virtual ~CE_Multigrid();
28 
29  /// Initialize array to given size, spacing, and boundary conditions.
30  void init( const UT_Vector3I &res,
31  const UT_Vector3 &spacing,
32  const UT_Vector3I &boundariesNeg,
33  const UT_Vector3I &boundariesPos,
34  int level = 0,
35  UT_Vector3I oddCoarsenings = UT_Vector3I(0, 0, 0));
36 
37  /// Initializes this grid from the provided UT_VoxelArray, assuming the
38  /// provided grid spacing.
39  void initFromVoxels(const UT_VoxelArrayF &f,
40  const UT_Vector3 &spacing,
41  const UT_Vector3I &boundariesNeg,
42  const UT_Vector3I &boundariesPos);
43 
44  /// Update the grid's ghost cells with values that enforce the boundary
45  /// conditions at its particular level in the multigrid hierarchy.
46  void updateBoundaryConditions() const;
47 
48  /// Calculates the maximum coarse grid elements in the coarsest grid, based
49  /// on the minPerAxis parameter.
50  int getMaxCoarse(int minPerAxis) const;
51 
52  /// Returns a grid's level within the multigrid hierarchy.
53  int getLevel() const { return myLevel; }
54 
55  /// Returns a vector containing the odd/even parity for each axis.
57  {
58  UT_Vector3I parity;
59  for (int i = 0; i < 3; i++)
60  parity[i] = getRes(i) % 2;
61  return parity;
62  }
63 
64  /// Initializes this to be the same dimensions, boundary conditions,
65  /// etc, of the given grid.
66  void match(const CE_Multigrid &src);
67 
68  /// Computes rhe residual r = b - A * x, where b == *this
69  void subtractApplyLaplacian(CE_Multigrid& x, CE_Multigrid &r) const;
70 
71 
72  /// Computes the inf and 2-norm of the residual r = b - A * x,
73  /// where b == *this.
74  void computeResidualNorms(const CE_Multigrid &x,
75  fpreal64 &norminf, fpreal64 &norm2) const;
76 
77  /// Performs red-black Gauss-Seidel smoothing, where b
78  /// == *this, and the implicit matrix is the Laplacian operator. If
79  /// useSmall is true, this use callsmoothLaplacianGaussSeidelSmall() if
80  /// possible.
81  void smoothLaplacianGaussSeidel(CE_Multigrid &x, int iterations = 1,
82  bool useSmall=true) const;
83 
84  /// Performs red-black Gauss-Seidel smoothing, using a kernel that only
85  /// works when the entire grid can fit in one workgroup on the OpenCL
86  /// device. Returns false if the criteria to run are not met.
87  bool smoothLaplacianGaussSeidelSmall(CE_Multigrid &x,
88  int iterations = 1) const;
89 
90  /// Directly solve the Poisson equation that this grid represents. This
91  /// function should only be called for small grids. At the moment this
92  /// always calls smoothLaplacianGaussSeidelSmall() with a large
93  /// number of iterations.
94  void directSolve(CE_Multigrid &x) const;
95 
96  /// Performs coarsening using full-weighting along one axis.
97  /// On output @code uh.getReas(axis) == getRes(axis) / 2 @endcode
98  void coarsenAlongAxis(CE_Multigrid& uh, int axis) const;
99 
100  /// Coarsen this grid into a lower resolution version of the grid.
101  UT_Vector3I coarsen(int minPerAxis, CE_Multigrid &uh) const;
102 
103  /// Interpolate this array into u using inverse of full-weighting along one
104  /// axis. On output
105  /// @code u.getRes(axis) == getRes(axis) * 2 + parity[axis] @endcode
106  void interpolateAlongAxis(CE_Multigrid& u, int axis) const;
107 
108 
109  /// Interpolate this grid into a fine version of the grid.
110  /// coarsenAxis is a vector of booleans specifying which axes to
111  /// interpolate
112  /// Where interpolateAxis[axis] is true, ensure that
113  /// @code u[axis] = getRes(axis) * 2 + parity[axis] @endcode
114  void interpolate(const UT_Vector3I &coarsenAxis, const UT_Vector3I &parity,
115  CE_Multigrid &u) const;
116 
117  /// Perform recursive V-cycle for multigrid algorithm, until reaching
118  /// numElements() < maxcoarse,
119  /// at which point a direct solve is done using directSolve.. At
120  /// each level, nSmoothDown and nSmoothUp iterations of Gauss-Seidel
121  /// smoothing are done, unless smoothTopLevelDown is false, in which
122  /// case no smoothing is done for the top-level grid on the down-cycle. This
123  /// can be used as an optimization when calling vcycle iteratively; the
124  /// previous vcycle call will have smoothed the top-level grid as its last
125  /// step, so you can often skip smoothing that same grid on the next vcycle
126  /// without loss of convergence.
127  void vcycle(int minPerAxis, int nSmoothDown, int nSmoothUp,
128  CE_Multigrid &x, bool smoothTopLevelDown = true) const;
129 
130  /// Performs recursive Full Multigrid to generate an initial guess for the
131  /// solution of the Poisson equation represented by this array. At
132  /// each level, nSmoothDown and nSmoothUp iterations of smoothing are done.
133  void fullMultigrid(int minPerAxis, int nSmoothDown,
134  int nSmoothUp, CE_Multigrid &x) const;
135 
136  /// Iteratively solve the Poisson equation, assuming *this is the right
137  /// hand side. This function will do at least miniter iterations, then
138  /// continue to iterate until the 2-norm of the residual has
139  /// been reduced by reltol, the inf-norm of the residual is below abstol,
140  /// or maxiters is reached.
141  /// The optional resnorm0 and resnorm2 parameters will contain the
142  /// 0- and 2-norms of the residuals for each iteration.
143  int solvePoisson(fpreal64 abstol, fpreal64 reltol, int miniter, int maxiter,
144  CE_Multigrid &x,
145  UT_ValArray<fpreal64> *resnorminf,
146  UT_ValArray<fpreal64> *resnorm2,
147  bool finishBetweenIters = false) const;
148 
149 protected:
150  void initLaplacian();
151 
152  fpreal32 myDiag, myOmega;
154  UT_Vector3I myBoundariesNeg, myBoundariesPos;
155  UT_Vector3 myBoundModNeg, myBoundModPos;
156  UT_Vector3 myBoundAdjModNeg, myBoundAdjModPos;
157  UT_Vector3 myDiagModNeg, myDiagModPos;
158  int myLevel;
159  bool myAllOpen, myAllClosed;
161 };
162 
163 
164 
165 #endif
166 #endif
#define CE_API
Definition: CE_API.h:10
UT_Vector3 myBoundModPos
Definition: CE_Multigrid.h:155
GLint level
Definition: glcorearb.h:107
void match(const CE_Grid &src)
Match the src CE_Grid in terms of size and border conditions.
fpreal32 myOmega
Definition: CE_Multigrid.h:152
UT_Vector3I getParity() const
Returns a vector containing the odd/even parity for each axis.
Definition: CE_Multigrid.h:56
UT_Vector3T< int64 > UT_Vector3I
png_uint_32 i
Definition: png.h:2877
int getLevel() const
Returns a grid's level within the multigrid hierarchy.
Definition: CE_Multigrid.h:53
GLfloat f
Definition: glcorearb.h:1925
UT_Vector3I myOddCoarsenings
Definition: CE_Multigrid.h:160
double fpreal64
Definition: SYS_Types.h:185
UT_Vector3 myDiagModPos
Definition: CE_Multigrid.h:157
UT_Vector3I myBoundariesPos
Definition: CE_Multigrid.h:154
UT_Vector3 mySpacing
Definition: CE_Multigrid.h:153
UT_Vector3 myBoundAdjModPos
Definition: CE_Multigrid.h:156
GLint GLenum GLint x
Definition: glcorearb.h:408
void initFromVoxels(const UT_VoxelArrayF &src, int xghost=1, int yghost=1, int zghost=1, int xpad=1, int ypad=1, int zpad=1)
GLboolean r
Definition: glcorearb.h:1221
UT_Vector3I getRes() const
Definition: CE_Grid.h:101
float fpreal32
Definition: SYS_Types.h:184
GLenum src
Definition: glcorearb.h:1792