tminres
SerialExample/ex2.cpp
Go to the documentation of this file.
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 }
 All Classes Files Functions Friends