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