Chombo + EB: BiCGStabSolver.H Source File

00001 #ifdef CH_LANG_CC
00002 /*
00003 *      _______              __
00004 *     / ___/ /  ___  __ _  / /  ___
00005 *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
00006 *    \___/_//_/\___/_/_/_/_.__/\___/
00007 *    Please refer to Copyright.txt, in Chombo's root directory.
00008 */
00009 #endif
00010 
00011 // BVS, June 18, 2003
00012 
00013 #ifndef _BICGSTABSOLVER_H_
00014 #define _BICGSTABSOLVER_H_
00015 
00016 #include "LinearSolver.H"
00017 #include "parstream.H"
00018 #include "CH_Timer.H"
00019 #include "NamespaceHeader.H"
00020 
00022 
00025 template <class T>
00026 class BiCGStabSolver : public LinearSolver<T>
00027 {
00028 public:
00029 
00030   BiCGStabSolver();
00031 
00032   virtual ~BiCGStabSolver();
00033 
00034   virtual void setHomogeneous(bool a_homogeneous)
00035      {m_homogeneous = a_homogeneous;}
00036 
00038 
00043   virtual void define(LinearOp<T>* a_op, bool a_homogeneous);
00044 
00046   virtual void solve(T& a_phi, const T& a_rhs);
00047 
00049   virtual void setConvergenceMetrics(Real a_metric,
00050                                      Real a_tolerance);
00051 
00053 
00057   bool m_homogeneous;
00058 
00060 
00063   LinearOp<T>* m_op;
00064 
00066 
00069   int m_imax;
00070 
00072 
00076   int m_verbosity;
00077 
00079 
00082   Real m_eps;
00083 
00085 
00088   Real m_reps;
00089 
00091 
00095   Real m_convergenceMetric;
00096 
00098 
00101   Real m_hang;
00102 
00104 
00111   int m_exitStatus;
00112 
00114 
00117   Real m_small;
00118 
00120 
00123   int m_numRestarts;
00124 
00126 
00130   int m_normType;
00131 
00132 };
00133 
00134 // *******************************************************
00135 // BiCGStabSolver Implementation
00136 // *******************************************************
00137 
00138 // For large elliptic problem, bottom smoother needed 67 iterations.
00139 // Bumping imax to 80.  (ndk)
00140 template <class T>
00141 BiCGStabSolver<T>::BiCGStabSolver()
00142   :m_homogeneous(false), m_op(NULL), m_imax(80), m_verbosity(3), m_eps(1.0E-6),
00143    m_reps(1.0E-12), m_convergenceMetric(-1.0), m_hang(1E-8), m_exitStatus(-1),
00144    m_small(1.0E-30), m_numRestarts(5), m_normType(2){;}
00145 
00146 template <class T>
00147 BiCGStabSolver<T>::~BiCGStabSolver()
00148 {
00149   m_op = NULL;
00150 }
00151 
00152 template <class T>
00153 void BiCGStabSolver<T>::define(LinearOp<T>* a_operator, bool a_homogeneous)
00154 {
00155   m_homogeneous = a_homogeneous;
00156   m_op = a_operator;
00157 }
00158 
00159 template <class T>
00160 void BiCGStabSolver<T>::solve(T& a_phi, const T& a_rhs)
00161 {
00162   CH_TIMERS("BiCGStabSolver::solve");
00163 
00164   CH_TIMER("BiCGStabSolver::solve::Initialize",timeInitialize);
00165   CH_TIMER("BiCGStabSolver::solve::MainLoop",timeMainLoop);
00166   CH_TIMER("BiCGStabSolver::solve::Cleanup",timeCleanup);
00167 
00168   CH_START(timeInitialize);
00169 
00170   T r, r_tilde, e, p, p_tilde, s_tilde, t, v;
00171 
00172   m_op->create(r,       a_rhs);
00173   m_op->create(r_tilde, a_rhs);
00174   m_op->create(e,       a_phi);
00175   m_op->create(p,       a_rhs);
00176   m_op->create(p_tilde, a_phi);
00177   m_op->create(s_tilde, a_phi);
00178   m_op->create(t,       a_rhs);
00179   m_op->create(v,       a_rhs);
00180 
00181   int recount = 0;
00182 
00183   CH_assert(m_op != NULL);
00184 
00185   m_op->residual(r, a_phi, a_rhs, m_homogeneous);
00186 
00187   m_op->assignLocal(r_tilde, r);
00188   m_op->setToZero(e);
00189   // (DFM 2/1/07) these next two need to be set to zero to prevent
00190   // problems in the multilevel case
00191   m_op->setToZero(p_tilde);
00192   m_op->setToZero(s_tilde);
00193 
00194   int i = 0;
00195 
00196   Real rho[4]={0,0,0,0}; // rho[0] = r_i , rho[1] = r_(i-1), etc.
00197   Real norm[2];
00198   norm[0] = m_op->norm(r, m_normType);
00199   Real initial_norm = norm[0];
00200   Real initial_rnorm = norm[0];
00201   norm[1] = norm[0];
00202 
00203   Real alpha[2] = {0,0};
00204   Real  beta[2] = {0,0};
00205   Real omega[2] = {0,0};
00206 
00207   bool init = true;
00208   int restarts = 0;
00209 
00210   if (m_verbosity >= 3)
00211     {
00212       pout() << "      BiCGStab:: initial Residual norm = "
00213              << initial_norm << "\n";
00214     }
00215 
00216   // if a convergence metric has been supplied, replace initial residual
00217   // with the supplied convergence metric...
00218   if (m_convergenceMetric > 0)
00219     {
00220       initial_norm = m_convergenceMetric;
00221     }
00222 
00223   CH_STOP(timeInitialize);
00224 
00225   CH_START(timeMainLoop);
00226   while ((i<m_imax && norm[0] > m_eps*norm[1]) && (norm[1] > 0))
00227     {
00228       i++;
00229 
00230       norm[1] = norm[0];
00231       alpha[1]= alpha[0];
00232       beta[1] = beta[0];
00233       omega[1]= omega[0];
00234 
00235       if (m_verbosity >= 5)
00236         {
00237           pout() << "      BiCGStab::       norm[0]  = "  << norm[0] << ", "
00238                  <<                        "norm[1]  = "  << norm[1]
00239                  << "\n";
00240           pout() << "      BiCGStab::       alpha[0] = "  << alpha[0] << ", "
00241                  <<                        "alpha[1] = "  << alpha[1]
00242                  << "\n";
00243           pout() << "      BiCGStab::       beta[0]  = "  << beta[0] << ", "
00244                  <<                        "beta[1]  = "  << beta[1]
00245                  << "\n";
00246           pout() << "      BiCGStab::       omega[0] = "  << omega[0] << ", "
00247                  <<                        "omega[1] = "  << omega[1]
00248                  << "\n";
00249         }
00250 
00251       rho[3] = rho[2];
00252       rho[2] = rho[1];
00253       rho[1] = m_op->dotProduct(r_tilde, r);
00254 
00255       if (m_verbosity >= 5)
00256         {
00257           pout() << "      BiCGStab::       rho[1] = "  << rho[1] << ", "
00258                  <<                        "rho[2] = "  << rho[2] << ", "
00259                  <<                        "rho[3] = "  << rho[3]
00260                  << "\n";
00261         }
00262 
00263       if (rho[1] == 0.0)
00264         {
00265           // we are finished, we will not converge anymore
00266           m_op->incr(a_phi, e, 1.0);
00267 
00268           if (m_verbosity >= 3)
00269             {
00270               pout() << "      BiCGStab:: rho = 0, returning"
00271                      << " -- Residual norm = "
00272                      << norm[0] << "\n";
00273             }
00274 
00275           CH_STOP(timeMainLoop);
00276 
00277           CH_START(timeCleanup);
00278 
00279           m_exitStatus = 2;
00280           // need to call clear just in case
00281           // this is not ideal -- maybe change the return
00282           // to exit
00283           m_op->clear(r);
00284           m_op->clear(r_tilde);
00285           m_op->clear(e);
00286           m_op->clear(p);
00287           m_op->clear(p_tilde);
00288           m_op->clear(s_tilde);
00289           m_op->clear(t);
00290           m_op->clear(v);
00291 
00292           CH_STOP(timeCleanup);
00293 
00294           return;
00295         }
00296 
00297       if (init)
00298         {
00299           m_op->assignLocal(p, r);
00300           init = false;
00301         }
00302       else
00303         {
00304           beta[1] = (rho[1]/rho[2])*(alpha[1]/omega[1]);
00305           m_op->scale(p, beta[1]);
00306           m_op->incr(p, v, -beta[1]*omega[1]);
00307           m_op->incr(p, r, 1.0);
00308         }
00309 
00310       if (m_verbosity >= 5)
00311         {
00312           pout() << "      BiCGStab::       beta[1]  = "  << beta[1]
00313                  << "\n";
00314         }
00315 
00316       m_op->preCond(p_tilde, p);
00317       m_op->applyOp(v, p_tilde, true);
00318       Real m = m_op->dotProduct(r_tilde, v);
00319       alpha[0] = rho[1]/m;
00320 
00321       if (m_verbosity >= 5)
00322         {
00323           pout() << "      BiCGStab::       rho[1] = "  << rho[1]   << ", "
00324                  <<                             "m = "  << m        << ", "
00325                  <<                      "alpha[0] = "  << alpha[0]
00326                  << "\n";
00327         }
00328 
00329       if (Abs(m) > m_small*Abs(rho[1]))
00330         {
00331           m_op->incr(r, v, -alpha[0]);
00332           norm[0] = m_op->norm(r, m_normType);
00333           m_op->incr(e, p_tilde, alpha[0]);
00334         }
00335       else
00336         {
00337           m_op->setToZero(r);
00338           norm[0] = 0.0;
00339         }
00340 
00341       if (norm[0] > m_eps*initial_norm && norm[0] > m_reps*initial_rnorm)
00342         {
00343           m_op->preCond(s_tilde, r);
00344           m_op->applyOp(t, s_tilde, true);
00345           omega[0] = m_op->dotProduct(t, r)/m_op->dotProduct(t, t);
00346           m_op->incr(e, s_tilde, omega[0]);
00347           m_op->incr(r, t,      -omega[0]);
00348           norm[0] = m_op->norm(r, m_normType);
00349         }
00350 
00351       if (m_verbosity >= 4)
00352         {
00353           pout() << "      BiCGStab::     iteration = "  << i        << ", error norm = " << norm[0] << ", rate = " << norm[1]/norm[0] << "\n";
00354         }
00355 
00356       if (norm[0] <= m_eps*initial_norm || norm[0] <= m_reps*initial_rnorm)
00357         {
00358           // converged to tolerance
00359           m_exitStatus = 1;
00360           break;
00361         }
00362 
00363       if (omega[0] == 0.0 || norm[0] > (1-m_hang)*norm[1])
00364         {
00365           if (recount == 0)
00366             {
00367               recount = 1;
00368             }
00369           else
00370             {
00371               recount = 0;
00372               m_op->incr(a_phi, e, 1.0);
00373 
00374               if (restarts == m_numRestarts)
00375                 {
00376                   if (m_verbosity >= 3)
00377                     {
00378                       pout() << "      BiCGStab:: max number of restarts reached"
00379                              << " -- Residual norm = "
00380                              << norm[0] << "\n";
00381                     }
00382 
00383                   CH_STOP(timeMainLoop);
00384 
00385                   CH_START(timeCleanup);
00386 
00387                   m_exitStatus = 3;
00388 
00389                   // need to call clear just in case
00390                   // this is not ideal -- maybe change the return
00391                   // to exit
00392                   m_op->clear(r);
00393                   m_op->clear(r_tilde);
00394                   m_op->clear(e);
00395                   m_op->clear(p);
00396                   m_op->clear(p_tilde);
00397                   m_op->clear(s_tilde);
00398                   m_op->clear(t);
00399                   m_op->clear(v);
00400 
00401                   CH_STOP(timeCleanup);
00402 
00403                   return;
00404                 }
00405 
00406               {
00407                 CH_TIME("BiCGStabSolver::solve::Restart");
00408 
00409                 m_op->residual(r, a_phi, a_rhs, m_homogeneous);
00410                 norm[0] = m_op->norm(r, m_normType);
00411                 rho[1]=0.0; rho[1]=0.0; rho[2]=0.0; rho[3]=0.0;
00412                 alpha[0]=0; beta[0]=0; omega[0]=0;
00413                 m_op->assignLocal(r_tilde, r);
00414                 m_op->setToZero(e);
00415 
00416                 restarts++;
00417               }
00418 
00419               if (m_verbosity >= 4)
00420                 {
00421                   pout() << "      BiCGStab::   restart =  " << restarts << "\n";
00422                 }
00423 
00424               init = true;
00425             }
00426         }
00427     }
00428   CH_STOP(timeMainLoop);
00429 
00430   CH_START(timeCleanup);
00431 
00432   if (m_verbosity >= 3)
00433     {
00434       pout() << "      BiCGStab:: " << i << " iterations, final Residual norm = "
00435              << norm[0] << "\n";
00436     }
00437 
00438   m_op->incr(a_phi, e, 1.0);
00439 
00440   m_op->clear(r);
00441   m_op->clear(r_tilde);
00442   m_op->clear(e);
00443   m_op->clear(p);
00444   m_op->clear(p_tilde);
00445   m_op->clear(s_tilde);
00446   m_op->clear(t);
00447   m_op->clear(v);
00448 
00449   CH_STOP(timeCleanup);
00450 }
00451 
00452 template <class T>
00453 void BiCGStabSolver<T>::setConvergenceMetrics(Real a_metric,
00454                                               Real a_tolerance)
00455 {
00456   m_convergenceMetric = a_metric;
00457   m_eps = a_tolerance;
00458 }
00459 
00460 #include "NamespaceFooter.H"
00461 #endif /*_BICGSTABSOLVER_H_*/