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