HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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  ~CE_Multigrid() override;
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  /// If struts is true, ghost cells that are outside two or more of the axes
47  /// are also updated (like + in the following diagram, where * represent the
48  /// ghost cells that are filled in regardless of the value of struts).
49  /// +*************+
50  /// * *
51  /// * ACTUAL GRID *
52  /// * *
53  /// +*************+
54  void updateBoundaryConditions(bool struts = false) const;
55 
56  /// Calculates the maximum coarse grid elements in the coarsest grid, based
57  /// on the minPerAxis parameter.
58  int getMaxCoarse(int minPerAxis) const;
59 
60  /// Returns a grid's level within the multigrid hierarchy.
61  int getLevel() const { return myLevel; }
62 
63  /// Returns a vector containing the odd/even parity for each axis.
65  {
66  UT_Vector3I parity;
67  for (int i = 0; i < 3; i++)
68  parity[i] = getRes(i) % 2;
69  return parity;
70  }
71 
72  /// Initializes this to be the same dimensions, boundary conditions,
73  /// etc, of the given grid.
74  void match(const CE_Multigrid &src);
75 
76  /// Computes rhe residual r = b - A * x, where b == *this
77  void subtractApplyLaplacian(CE_Multigrid& x, CE_Multigrid &r) const;
78 
79 
80  /// Computes the inf and 2-norm of the residual r = b - A * x,
81  /// where b == *this.
82  void computeResidualNorms(const CE_Multigrid &x,
83  fpreal64 &norminf, fpreal64 &norm2) const;
84 
85  /// Performs red-black Gauss-Seidel smoothing, where b
86  /// == *this, and the implicit matrix is the Laplacian operator. If
87  /// useSmall is true, this use callsmoothLaplacianGaussSeidelSmall() if
88  /// possible.
89  void smoothLaplacianGaussSeidel(CE_Multigrid &x, int iterations = 1,
90  bool useSmall=true) const;
91 
92  /// Performs red-black Gauss-Seidel smoothing, using a kernel that only
93  /// works when the entire grid can fit in one workgroup on the OpenCL
94  /// device. Returns false if the criteria to run are not met.
95  bool smoothLaplacianGaussSeidelSmall(CE_Multigrid &x,
96  int iterations = 1) const;
97 
98  /// Directly solve the Poisson equation that this grid represents. This
99  /// function should only be called for small grids. At the moment this
100  /// always calls smoothLaplacianGaussSeidelSmall() with a large
101  /// number of iterations.
102  void directSolve(CE_Multigrid &x) const;
103 
104  /// Performs coarsening using full-weighting along one axis.
105  /// On output @code uh.getReas(axis) == getRes(axis) / 2 @endcode
106  void coarsenAlongAxis(CE_Multigrid& uh, int axis) const;
107 
108  /// Coarsen this grid into a lower resolution version of the grid.
109  UT_Vector3I coarsen(int minPerAxis, CE_Multigrid &uh) const;
110 
111  /// Interpolate this array into u using inverse of full-weighting along one
112  /// axis. On output
113  /// @code u.getRes(axis) == getRes(axis) * 2 + parity[axis] @endcode
114  void interpolateAlongAxis(CE_Multigrid& u, int axis) const;
115 
116 
117  /// Interpolate this grid into a fine version of the grid.
118  /// coarsenAxis is a vector of booleans specifying which axes to
119  /// interpolate
120  /// Where interpolateAxis[axis] is true, ensure that
121  /// @code u[axis] = getRes(axis) * 2 + parity[axis] @endcode
122  void interpolate(const UT_Vector3I &coarsenAxis, const UT_Vector3I &parity,
123  CE_Multigrid &u) const;
124 
125  /// Perform recursive V-cycle for multigrid algorithm, until reaching
126  /// numElements() < maxcoarse,
127  /// at which point a direct solve is done using directSolve.. At
128  /// each level, nSmoothDown and nSmoothUp iterations of Gauss-Seidel
129  /// smoothing are done, unless smoothTopLevelDown is false, in which
130  /// case no smoothing is done for the top-level grid on the down-cycle. This
131  /// can be used as an optimization when calling vcycle iteratively; the
132  /// previous vcycle call will have smoothed the top-level grid as its last
133  /// step, so you can often skip smoothing that same grid on the next vcycle
134  /// without loss of convergence.
135  void vcycle(int minPerAxis, int nSmoothDown, int nSmoothUp,
136  CE_Multigrid &x, bool smoothTopLevelDown = true) const;
137 
138  /// Performs recursive Full Multigrid to generate an initial guess for the
139  /// solution of the Poisson equation represented by this array. At
140  /// each level, nSmoothDown and nSmoothUp iterations of smoothing are done.
141  void fullMultigrid(int minPerAxis, int nSmoothDown,
142  int nSmoothUp, CE_Multigrid &x) const;
143 
144  /// Iteratively solve the Poisson equation, assuming *this is the right
145  /// hand side. This function will do at least miniter iterations, then
146  /// continue to iterate until the 2-norm of the residual has
147  /// been reduced by reltol, the inf-norm of the residual is below abstol,
148  /// or maxiters is reached.
149  /// The optional resnorm0 and resnorm2 parameters will contain the
150  /// 0- and 2-norms of the residuals for each iteration.
151  int solvePoisson(fpreal64 abstol, fpreal64 reltol, int miniter, int maxiter,
152  CE_Multigrid &x,
153  UT_ValArray<fpreal64> *resnorminf,
154  UT_ValArray<fpreal64> *resnorm2,
155  bool finishBetweenIters = false) const;
156 
157 protected:
158  /// Coarsen this grid into a lower resolution version of the grid, by explicitly
159  /// coarsening each axis sequentially.
160  UT_Vector3I coarsenByAxis(int minPerAxis, CE_Multigrid &uh) const;
161  /// Interpolate this grid into a fine version of the grid. This version explicitly
162  /// interpolates each axis sequentially.
163  void interpolateByAxis(const UT_Vector3I &coarsenAxis, const UT_Vector3I &parity,
164  CE_Multigrid &u) const;
165 
166 protected:
167  void initLaplacian();
168 
169  fpreal32 myDiag, myOmega;
171  UT_Vector3I myBoundariesNeg, myBoundariesPos;
172  UT_Vector3 myBoundModNeg, myBoundModPos;
173  UT_Vector3 myBoundAdjModNeg, myBoundAdjModPos;
174  UT_Vector3 myDiagModNeg, myDiagModPos;
175  int myLevel;
176  bool myAllOpen, myAllClosed;
178 };
179 
180 
181 
182 #endif
183 #endif
#define CE_API
Definition: CE_API.h:10
UT_Vector3 myBoundModPos
Definition: CE_Multigrid.h:172
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:169
UT_Vector3I getParity() const
Returns a vector containing the odd/even parity for each axis.
Definition: CE_Multigrid.h:64
UT_Vector3T< int64 > UT_Vector3I
GLenum src
Definition: glcorearb.h:1792
float fpreal32
Definition: SYS_Types.h:200
GLint GLenum GLint x
Definition: glcorearb.h:408
double fpreal64
Definition: SYS_Types.h:201
int getLevel() const
Returns a grid's level within the multigrid hierarchy.
Definition: CE_Multigrid.h:61
UT_Vector3I myOddCoarsenings
Definition: CE_Multigrid.h:177
UT_Vector3 myDiagModPos
Definition: CE_Multigrid.h:174
UT_Vector3I myBoundariesPos
Definition: CE_Multigrid.h:171
GLuint res
Definition: glew.h:11549
UT_Vector3 mySpacing
Definition: CE_Multigrid.h:170
UT_Vector3 myBoundAdjModPos
Definition: CE_Multigrid.h:173
GLfloat f
Definition: glcorearb.h:1925
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