tminres
|
00001 // tminres is free software; you can redistribute it and/or modify it under the 00002 // terms of the GNU Lesser General Public License (as published by the Free 00003 // Software Foundation) version 2.1 dated February 1999. 00004 // 00005 // Authors: 00006 // - Umberto Villa, Emory University - uvilla@emory.edu 00007 // - Michael Saunders, Stanford University 00008 // - Santiago Akle, Stanford University 00009 00037 #include <tminres.hpp> 00038 #include "SimpleVector.hpp" 00039 #include <cmath> 00040 #include <cassert> 00041 00043 00049 class StokesOperator { 00050 public: 00051 00053 00056 StokesOperator(int n_ ); 00057 00059 virtual ~StokesOperator(); 00060 00062 void Apply(const SimpleVector & X, SimpleVector & Y) const; 00063 00064 private: 00065 00066 inline int getUIndex(int i, int j) const 00067 { 00068 if(i == -1) 00069 i = n; 00070 00071 if(i==n+1) 00072 i = 0; 00073 00074 if(j== -1) 00075 j = n-1; 00076 00077 if(j==n) 00078 j = 0; 00079 00080 return offsetU + i + (n+1)*j; 00081 } 00082 00083 inline int getVIndex(int i, int j) const 00084 { 00085 if(i == -1) 00086 i = n-1; 00087 00088 if(i==n) 00089 i = 0; 00090 00091 if(j== -1) 00092 j = n; 00093 00094 if(j==n+1) 00095 j = 0; 00096 00097 return offsetV + i + n*j; 00098 } 00099 00100 inline int getPIndex(int i, int j) const 00101 { 00102 if(i == -1 || i == n || j == -1 || j == n) 00103 return -1; 00104 00105 return offsetP + i + n*j; 00106 } 00107 00108 00109 int n; 00110 double h; 00111 int offsetU; 00112 int offsetV; 00113 int offsetP; 00114 }; 00115 00117 00120 class StokesDiagonalScaling 00121 { 00122 public: 00124 00127 StokesDiagonalScaling(int n_); 00128 00129 virtual ~StokesDiagonalScaling(); 00130 00132 void Apply(const SimpleVector & X, SimpleVector & Y) const; 00133 00134 private: 00135 int n; 00136 double h; 00137 int offsetU; 00138 int offsetV; 00139 int offsetP; 00140 00141 }; 00142 00143 00144 00145 int main() 00146 { 00147 //(1) Discretization parameters and problem size 00148 int n(4); //number of cell for each edge 00149 int dim(n*n + 2*(n+1)*n); //total number of unknows (n^2 pressure and (n+1)n for each velocity component) 00150 //(2) Define the linear operator "op" we want to solve. 00151 StokesOperator op(n); 00152 //(3) Define the exact solution (at random) 00153 SimpleVector sol(dim); 00154 sol.Randomize( 1 ); 00155 00156 //(4) Define the "rhs" as "rhs = op*sol" 00157 SimpleVector rhs(dim); 00158 op.Apply(sol, rhs); 00159 double rhsNorm( sqrt(InnerProduct(rhs,rhs)) ); 00160 std::cout << "|| rhs || = " << rhsNorm << "\n"; 00161 00162 //(5) Define the preconditioner 00163 StokesDiagonalScaling prec(n); 00164 00165 //(6) Use an identically zero initial guess 00166 SimpleVector x(dim); 00167 x = 0; 00168 00169 //(7) Set the minres parameters 00170 double shift(0); 00171 int max_iter(100); 00172 double tol(1e-6); 00173 bool show(true); 00174 00175 //(8) Solve the problem with minres 00176 MINRES(op, x, rhs, &prec, shift, max_iter, tol, show); 00177 00178 //(9) Compute the error || x_ex - x_minres ||_2 00179 subtract(x, sol, x); 00180 double err2 = InnerProduct(x,x); 00181 00182 std::cout<< "|| x_ex - x_n || = " << sqrt(err2) << "\n"; 00183 00184 return 0; 00185 } 00186 00187 00188 StokesOperator::StokesOperator(int n_): 00189 n(n_), 00190 h(1./static_cast<double>(n_)), 00191 offsetU(0), 00192 offsetV( (n+1)*n ), 00193 offsetP( (n+1)*n + n*(n+1) ) 00194 { 00195 assert(n > 0); 00196 } 00197 00198 StokesOperator::~StokesOperator() 00199 { 00200 // TODO Auto-generated destructor stub 00201 } 00202 00203 void StokesOperator::Apply(const SimpleVector & X, SimpleVector & Y) const 00204 { 00205 Y = 0.; 00206 for(int i(0); i<n+1; ++i) 00207 for(int j(0); j<n+1; ++j) 00208 { 00209 double uij ( X[getUIndex(i ,j ) ] ); 00210 double uimj( X[getUIndex(i-1,j ) ] ); 00211 double uipj( X[getUIndex(i+1,j ) ] ); 00212 double uijm( X[getUIndex(i ,j-1) ] ); 00213 double uijp( X[getUIndex(i ,j+1) ] ); 00214 00215 double vij ( X[getVIndex(i ,j ) ] ); 00216 double vimj( X[getVIndex(i-1,j ) ] ); 00217 double vipj( X[getVIndex(i+1,j ) ] ); 00218 double vijm( X[getVIndex(i ,j-1) ] ); 00219 double vijp( X[getVIndex(i ,j+1) ] ); 00220 00221 double pij ( X.at(getPIndex(i ,j ) ) ); 00222 double pimj( X.at(getPIndex(i-1,j ) ) ); 00223 double pijm( X.at(getPIndex(i ,j-1) ) ); 00224 00225 double laplacianu( 4.*uij - uimj - uipj - uijm - uijp ); 00226 laplacianu *= (1/h/h); 00227 double laplacianv( 4.*vij - vimj - vipj - vijm - vijp ); 00228 laplacianv *= (1/h/h); 00229 00230 double dxp( pij - pimj); 00231 dxp *= 1./h; 00232 double dyp( pijm - pij); 00233 dyp *= 1./h; 00234 00235 double div( uipj - uij + vij - vijp); 00236 div *= -1./h; 00237 00238 if ( j != n ) 00239 Y[getUIndex(i,j)] += laplacianu + dxp + uij; 00240 00241 if ( i !=n ) 00242 Y[getVIndex(i,j)] += laplacianv + dyp + vij; 00243 00244 if( i != n && j != n ) 00245 Y[getPIndex(i,j)] += div; 00246 00247 00248 } 00249 } 00250 00251 00252 StokesDiagonalScaling::StokesDiagonalScaling(int n_): 00253 n(n_), 00254 h(1./static_cast<double>(n_)), 00255 offsetU(0), 00256 offsetV( (n+1)*n ), 00257 offsetP( (n+1)*n + n*(n+1) ) 00258 { 00259 assert(n > 0); 00260 } 00261 00262 00263 StokesDiagonalScaling::~StokesDiagonalScaling() 00264 { 00265 00266 } 00267 00268 void StokesDiagonalScaling::Apply(const SimpleVector & X, SimpleVector & Y) const 00269 { 00270 int i(0); 00271 00272 double diagValue(4./h/h + 1.); 00273 for( ; i < offsetP; ++i) 00274 Y[i] = X[i]/diagValue; 00275 00276 for( ; i < offsetP + n*n; ++i) 00277 Y[i] = X[i]; 00278 }