1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
|
// DO NOT EDIT !
// This file is generated using the MantaFlow preprocessor (prep generate).
/******************************************************************************
*
* MantaFlow fluid solver framework
* Copyright 2011 Tobias Pfaff, Nils Thuerey
*
* This program is free software, distributed under the terms of the
* Apache License, Version 2.0
* http://www.apache.org/licenses/LICENSE-2.0
*
* Multigrid solver by Florian Ferstl (florian.ferstl.ff@gmail.com)
*
* This is an implementation of the solver developed by Dick et al. [1]
* without topology awareness (= vertex duplication on coarser levels). This
* simplification allows us to use regular grids for all levels of the multigrid
* hierarchy and works well for moderately complex domains.
*
* [1] Solving the Fluid Pressure Poisson Equation Using Multigrid-Evaluation
* and Improvements, C. Dick, M. Rogowsky, R. Westermann, IEEE TVCG 2015
*
******************************************************************************/
#ifndef _MULTIGRID_H
#define _MULTIGRID_H
#include "vectorbase.h"
#include "grid.h"
namespace Manta {
//! Multigrid solver
class GridMg {
public:
//! constructor: preallocates most of required memory for multigrid hierarchy
GridMg(const Vec3i &gridSize);
~GridMg(){};
//! update system matrix A from symmetric 7-point stencil
void setA(const Grid<Real> *pA0,
const Grid<Real> *pAi,
const Grid<Real> *pAj,
const Grid<Real> *pAk);
//! set right-hand side after setting A
void setRhs(const Grid<Real> &rhs);
bool isASet() const
{
return mIsASet;
}
bool isRhsSet() const
{
return mIsRhsSet;
}
//! perform VCycle iteration
// - if src is null, then a zero vector is used instead
// - returns norm of residual after VCylcle
Real doVCycle(Grid<Real> &dst, const Grid<Real> *src = nullptr);
// access
void setCoarsestLevelAccuracy(Real accuracy)
{
mCoarsestLevelAccuracy = accuracy;
}
Real getCoarsestLevelAccuracy() const
{
return mCoarsestLevelAccuracy;
}
void setSmoothing(int numPreSmooth, int numPostSmooth)
{
mNumPreSmooth = numPreSmooth;
mNumPostSmooth = numPostSmooth;
}
int getNumPreSmooth() const
{
return mNumPreSmooth;
}
int getNumPostSmooth() const
{
return mNumPostSmooth;
}
//! Set factor for automated downscaling of trivial equations:
// 1*x_i = b_i ---> trivialEquationScale*x_i = trivialEquationScale*b_i
// Factor should be significantly smaller than the scale of the entries in A.
// Info: Trivial equations of the form x_i = b_i can have a negative
// effect on the coarse grid operators of the multigrid hierarchy (due
// to scaling mismatches), which can lead to slow multigrid convergence.
// To avoid this, the solver checks for such equations when updating A
// (and rhs) and scales these equations by a fixed factor < 1.
void setTrivialEquationScale(Real scale)
{
mTrivialEquationScale = scale;
}
private:
Vec3i vecIdx(int v, int l) const
{
return Vec3i(v % mSize[l].x,
(v % (mSize[l].x * mSize[l].y)) / mSize[l].x,
v / (mSize[l].x * mSize[l].y));
}
int linIdx(Vec3i V, int l) const
{
return V.x + V.y * mPitch[l].y + V.z * mPitch[l].z;
}
bool inGrid(Vec3i V, int l) const
{
return V.x >= 0 && V.y >= 0 && V.z >= 0 && V.x < mSize[l].x && V.y < mSize[l].y &&
V.z < mSize[l].z;
}
void analyzeStencil(int v, bool is3D, bool &isStencilSumNonZero, bool &isEquationTrivial) const;
void genCoarseGrid(int l);
void genCoraseGridOperator(int l);
void smoothGS(int l, bool reversedOrder);
void calcResidual(int l);
Real calcResidualNorm(int l);
void solveCG(int l);
void restrict(int l_dst, const std::vector<Real> &src, std::vector<Real> &dst) const;
void interpolate(int l_dst, const std::vector<Real> &src, std::vector<Real> &dst) const;
private:
enum VertexType : char {
vtInactive = 0,
vtActive = 1,
vtActiveTrivial = 2, // only on finest level 0
vtRemoved = 3, //-+
vtZero = 4, // +-- only during coarse grid generation
vtFree = 5 //-+
};
struct CoarseningPath {
Vec3i U, W, N;
int sc, sf;
Real rw, iw;
bool inUStencil;
};
int mNumPreSmooth;
int mNumPostSmooth;
Real mCoarsestLevelAccuracy;
Real mTrivialEquationScale;
std::vector<std::vector<Real>> mA;
std::vector<std::vector<Real>> mx;
std::vector<std::vector<Real>> mb;
std::vector<std::vector<Real>> mr;
std::vector<std::vector<VertexType>> mType;
std::vector<std::vector<double>> mCGtmp1, mCGtmp2, mCGtmp3, mCGtmp4;
std::vector<Vec3i> mSize, mPitch;
std::vector<CoarseningPath> mCoarseningPaths0;
bool mIs3D;
int mDim;
int mStencilSize;
int mStencilSize0;
Vec3i mStencilMin;
Vec3i mStencilMax;
bool mIsASet;
bool mIsRhsSet;
// provide kernels with access
friend struct knActivateVertices;
friend struct knActivateCoarseVertices;
friend struct knSetRhs;
friend struct knGenCoarseGridOperator;
friend struct knSmoothColor;
friend struct knCalcResidual;
friend struct knResidualNormSumSqr;
friend struct knRestrict;
friend struct knInterpolate;
}; // GridMg
} // namespace Manta
#endif
|
