linear_solvers/TrilinosSolver_test.cc File Reference

#include "generic.h"
#include "poisson.h"
#include "meshes/rectangular_quadmesh.h"

Classes
class	PoissonProblem< ELEMENT >
	Micky mouse Poisson problem. More...

Namespaces
	TanhSolnForPoisson
	Namespace for exact solution for Poisson equation with "sharp step".

Functions
void	TanhSolnForPoisson::get_exact_u (const Vector< double > &x, Vector< double > &u)
	Exact solution as a Vector. More...
void	TanhSolnForPoisson::source_function (const Vector< double > &x, double &source)
	Source function required to make the solution above an exact solution. More...
int	main (int argc, char **argv)
	Driver code for 2D Poisson problem. More...

Function Documentation

main()

int main	(	int argc	,
		char **	argv
	)

Driver code for 2D Poisson problem.

{
  // Note: there is a memory leak in trilinos itself as 
  // demonstrated by this bit of code!
  if (false)
   {
    // set up a parameter list
    Teuchos::ParameterList ml_parameters;
    ML_Epetra::SetDefaults("SA", ml_parameters);
    
    std::cout << "hello" << std::endl;
    exit(0);
   }
 
 // Set the orientation of the "step" to 45 degrees
 TanhSolnForPoisson::TanPhi=1.0;
 
 // Initial value for the steepness of the "step"
 TanhSolnForPoisson::Alpha=1.0;
  
 // Create the problem with 2D nine-node elements from the
 // QPoissonElement family. Pass pointer to source function.
 PoissonProblem<RefineableQPoissonElement<2,2> >
  problem(&TanhSolnForPoisson::source_function);   
 
 
 DocInfo doc_info;
 doc_info.set_directory("RESLT");
 
 // File to report the number of Newton iterations
 ofstream conv_file;
 char filename[100];
 sprintf(filename,"%s/conv.dat",doc_info.directory().c_str());
 conv_file.open(filename);
 
 
 // Trilinos solvers with the ML preconditioner
 //--------------------------------------------
 {
 
  cout << "==================================="
       << "=================================" << endl;
  cout << "Testing Trilinos solver preconditioned with "
       << "TrilinosMLPreconditioner" << endl;
  cout << "======================================="
       << "=============================" << endl;
 
  // Create a Trilinos Solver 
  TrilinosAztecOOSolver* linear_solver_pt = new TrilinosAztecOOSolver;
  
  // Create the Trilinos ML preconditioner
  TrilinosMLPreconditioner* preconditioner_pt = new TrilinosMLPreconditioner;
 
  // Set the preconditioner pointer
  linear_solver_pt->preconditioner_pt() = preconditioner_pt;
 
  // Set linear solver
  problem.linear_solver_pt() = linear_solver_pt;
  
  linear_solver_pt->enable_doc_time();
  
  // Loop over specific built-in Trilinos solvers
  for (unsigned s=0; s<=2; s++)
   {
    // Select the linear solver
    switch(s)
     {
     case 0: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::CG;
      cout << "Using CG as solver" << endl;
      cout << "==================" << endl;
      break;
     case 1: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::GMRES;
      cout << "Using GMRES as solver" << endl;
      cout << "=====================" << endl;
      break;
     case 2: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::BiCGStab;
      cout << "Using BiCGStab as solver" << endl;
      cout << "========================" << endl;
      break;
     default: cout << "Error selecting solver" << endl;
     } 
    
    // Solve the problem
     problem.newton_solve(); 
    
     // Doc
     problem.doc_solution(doc_info);
     doc_info.number()++;
     oomph_info << "Number of Newton iterations: " << problem.newton_count() 
                << std::endl << std::endl;
     conv_file << problem.newton_count() << std::endl;
 
   } // end of loop over preconditioner
 
  // delete the Trilinos solver
  delete linear_solver_pt;
  delete preconditioner_pt;  
 }
 
 
 
 // Trilinos solvers with the IFPACK preconditioner
 //------------------------------------------------
 {
 
  cout << "================================="
       << "=======================================" << endl;
  cout << "Testing Trilinos solver preconditioned" 
       << " with TrilinosIFPACKPreconditioner" << endl;
  cout << "==========================================="
       << "=============================" << endl;
 
    // Create a Trilinos Solver 
  TrilinosAztecOOSolver* linear_solver_pt = new TrilinosAztecOOSolver;
  
  // create the Trilinos IFPACK preconditioner
  TrilinosIFPACKPreconditioner* preconditioner_pt = 
   new TrilinosIFPACKPreconditioner;
 
  // set the preconditioner pointer
  linear_solver_pt->preconditioner_pt() = preconditioner_pt;
  
  // Set linear solver
  problem.linear_solver_pt() = linear_solver_pt;
  
  linear_solver_pt->enable_doc_time();
  
  // Loop over specific built-in Trilinos solvers
  for (unsigned s=0; s<=2; s++)
   {
    // Select the linear solver
    switch(s)
     {
     case 0: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::CG;
      cout << "Using CG as solver" << endl;
      cout << "==================" << endl;
      break;
     case 1: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::GMRES;
      cout << "Using GMRES as solver" << endl;
      cout << "=====================" << endl;
      break;
     case 2: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::BiCGStab;
      cout << "Using BiCGStab as solver" << endl;
      cout << "========================" << endl;
      break;
     default: cout << "Error selecting solver" << endl;
     } 
    
    // Solve the problem
     problem.newton_solve();
 
     // Doc
     problem.doc_solution(doc_info);
     doc_info.number()++;
     oomph_info << "Number of Newton iterations: " << problem.newton_count() 
                << std::endl << std::endl;
     conv_file << problem.newton_count() << std::endl;
    
   } // end of loop over preconditioner
 
  // delete the Trilinos solver
  delete linear_solver_pt;
  delete preconditioner_pt;
 }
 
 
 
 
 
 // Test the oomph-lib iterative linear solvers with Trilinos preconditioners
 //--------------------------------------------------------------------------
 
 cout
  << "======================================" 
  << "==================================" 
  << endl;
 cout 
  << "Testing oomph-lib iterative linear solvers " 
  << "with Trilinos preconditioners"
  << endl;
 cout
  << "=======================================" 
  << "=================================" 
  << endl  << endl  << endl;
 
 
 
 // Pointer to linear solver
 LinearSolver* linear_solver_pt=0;
 
 // Loop over two Trilonos preconditioners
 for (unsigned i_prec=0;i_prec<2;i_prec++)
  {
 
   // Pointer to preconditioner
   Preconditioner* preconditioner_pt=0;
 
   if (i_prec==0)
    {
     // Create a IFPACK preconditioner
     preconditioner_pt = new TrilinosIFPACKPreconditioner;
     cout << "Using  TrilinosIFPACKPreconditioner as preconditioner" << endl;
     cout << "=====================================================" << endl;
    }
   else
    {
     // Create a ML preconditioner
     preconditioner_pt = new TrilinosMLPreconditioner;
     cout << "Using  TrilinosMLPreconditioner as preconditioner" << endl;
     cout << "=================================================" << endl;
    }
 
 
   // Test CG
   //--------
   cout << "CG iterative solver" << endl;
   cout << "-------------------" << endl;
   linear_solver_pt = new CG<CRDoubleMatrix>;
   static_cast<IterativeLinearSolver*>(linear_solver_pt)->tolerance()=1.0e-10;
   problem.linear_solver_pt() = linear_solver_pt;
   
   // set the preconditioner in the iterative solver
   static_cast<IterativeLinearSolver*>(linear_solver_pt)->
    preconditioner_pt()= preconditioner_pt;
   
   // Solve the problem
   problem.newton_solve(); 
  
 
   // Doc
   problem.doc_solution(doc_info);
   doc_info.number()++;
   oomph_info << "Number of Newton iterations: " << problem.newton_count() 
              << std::endl << std::endl;
   conv_file << problem.newton_count() << std::endl;
 
   delete linear_solver_pt;
   
   
   // Test BiCGStab
   //--------------
   cout << "BiCGStab iterative solver" << endl;
   cout << "-------------------------" << endl;
   linear_solver_pt = new BiCGStab<CRDoubleMatrix>;
   static_cast<IterativeLinearSolver*>(linear_solver_pt)->tolerance()=1.0e-10;
   problem.linear_solver_pt() = linear_solver_pt;
   
   // set the preconditioner in the iterative solver
   static_cast<IterativeLinearSolver*>(linear_solver_pt)->preconditioner_pt()
    = preconditioner_pt;
   
   // Solve the problem
   problem.newton_solve(); 
  
   // Doc
   problem.doc_solution(doc_info);
   doc_info.number()++;
   oomph_info << "Number of Newton iterations: " << problem.newton_count() 
              << std::endl << std::endl;
   conv_file << problem.newton_count() << std::endl;
 
   delete linear_solver_pt;
 
   
   // Test GMRES
   //-----------
   cout << "GMRES iterative solver" << endl;
   cout << "----------------------" << endl;
   linear_solver_pt = new GMRES<CRDoubleMatrix>;
   static_cast<IterativeLinearSolver*>(linear_solver_pt)->tolerance()=1.0e-10;
   problem.linear_solver_pt() = linear_solver_pt;
     
   // set the preconditioner in the iterative solver
   static_cast<IterativeLinearSolver*>(linear_solver_pt)->preconditioner_pt()
    = preconditioner_pt;
  
   // Solve the problem
   problem.newton_solve();
 
   // Doc
   problem.doc_solution(doc_info);
   doc_info.number()++;
   oomph_info << "Number of Newton iterations: " << problem.newton_count() 
              << std::endl << std::endl;
   conv_file << problem.newton_count() << std::endl;
 
   delete linear_solver_pt;
 
   // Kill preconditioner
   delete preconditioner_pt;
   preconditioner_pt=0;
 
  } // end of loop over preconditioners
 
 
 // Trilinos solvers with the OOMPH-lib
 //------------------------------------------------
 {
 
  cout << "==================================" 
       << "=======================================" << endl;
  cout << "Testing Trilinos solver preconditioned with" 
       << " MatrixBasedDiagPreconditioner" << endl;
  cout << "=======================================" 
       << "==================================" << endl;
 
  // Create a Trilinos Solver 
  TrilinosAztecOOSolver* linear_solver_pt = new TrilinosAztecOOSolver;
  
  // create the oomphlib diagonal preconditioner
  Preconditioner* preconditioner_pt = 
   new MatrixBasedDiagPreconditioner;
 
  // set the preconditioner pointer
  linear_solver_pt->preconditioner_pt() = preconditioner_pt;
 
  // Set linear solver
  problem.linear_solver_pt() = linear_solver_pt;
  
  linear_solver_pt->enable_doc_time();
  
  // Loop over specific built-in Trilinos solvers
  for (unsigned s=0; s<=2; s++)
   {
        // Select the linear solver
    switch(s)
     {
     case 0: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::CG;
      cout << "Using CG as solver" << endl;
      cout << "==================" << endl;
      break;
     case 1: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::GMRES;
      cout << "Using GMRES as solver" << endl;
      cout << "=====================" << endl;
      break;
     case 2: 
      linear_solver_pt->solver_type() = TrilinosAztecOOSolver::BiCGStab;
      cout << "Using BiCGStab as solver" << endl;
      cout << "========================" << endl;
      break;
     default: cout << "Error selecting solver" << endl;
     } 
    
    // Solve the problem
    problem.newton_solve();
    
    // Doc
    problem.doc_solution(doc_info);
    doc_info.number()++;
    oomph_info << "Number of Newton iterations: " << problem.newton_count() 
               << std::endl << std::endl;
    conv_file << problem.newton_count() << std::endl;
    
   } // end of loop over preconditioner
 
  // delete the Trilinos solver
  delete linear_solver_pt;
  
  // Delete oomph-lib preconditioner
  delete preconditioner_pt;
 }
 
 conv_file.close();
 
 // This code is not used, just compiled -- don't delete as it's
 // extracted by doxygen for a tutorial
 bool never_get_here=true;
 if (!never_get_here)
  {
   // Create oomph-lib linear solver
   IterativeLinearSolver* linear_solver_pt=new GMRES<CRDoubleMatrix>;
 
   // Create Trilinos IFPACK preconditioner as oomph-lib Preconditioner
   Preconditioner* preconditioner_pt=new TrilinosIFPACKPreconditioner;
 
   // Pass pointer to preconditioner to oomph-lib IterativeLinearSolver
   linear_solver_pt->preconditioner_pt()=preconditioner_pt;
  }
 
} //end of main

References TanhSolnForPoisson::Alpha, oomph::DocInfo::directory(), oomph::LinearSolver::enable_doc_time(), MergeRestartFiles::filename, oomph::DocInfo::number(), oomph::oomph_info, oomph::IterativeLinearSolver::preconditioner_pt(), problem, s, oomph::DocInfo::set_directory(), oomph::TrilinosAztecOOSolver::solver_type(), TanhSolnForPoisson::source_function(), and TanhSolnForPoisson::TanPhi.