adv_diff_iterative_linear_solver_tester.cc File Reference

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

Classes
class	TwoMeshFluxAdvectionDiffusionProblem< ELEMENT >

Namespaces
	TanhSolnForAdvectionDiffusion

Functions
void	TanhSolnForAdvectionDiffusion::get_exact_u (const Vector< double > &x, Vector< double > &u)
	Exact solution as a Vector. More...
void	TanhSolnForAdvectionDiffusion::get_exact_u (const Vector< double > &x, double &u)
	Exact solution as a scalar. More...
void	TanhSolnForAdvectionDiffusion::source_function (const Vector< double > &x_vect, double &source)
	Source function required to make the solution above an exact solution. More...
void	TanhSolnForAdvectionDiffusion::prescribed_flux_on_fixed_x_boundary (const Vector< double > &x, double &flux)
	Flux required by the exact solution on a boundary on which x is fixed. More...
void	TanhSolnForAdvectionDiffusion::wind_function (const Vector< double > &x, Vector< double > &wind)
	Wind. More...
void	run_it (LinearSolver *linear_solver_pt)
int	main (int argc, char *argv[])

Function Documentation

main()

int main	(	int argc	,
		char *	argv[]
	)

Check out various iterative linear solvers on an advection diffusion problem. Command line arg specifies the solver.

Possible endless loop for memory leak testing

{
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 //Check
 if (CommandLineArgs::Argc!=2)
  {
   throw OomphLibError(
    "Code needs one argument to specify which solver to use",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
 
 // Default value for Peclet number -- can be reset to zero for
 // solvers that require symmetric matrices
 TanhSolnForAdvectionDiffusion::Peclet=200.0;
 
 // Tolerance for iterative linear solvers
 double tol=1.0e-12;
 
 // Pointer to linear solver
 LinearSolver* linear_solver_pt=new SuperLUSolver;
 
 // Pointer to iterative linear solver
 IterativeLinearSolver* it_linear_solver_pt=0;
 
 // Preconditioner
 Preconditioner* prec_pt=new ILUZeroPreconditioner<CRDoubleMatrix>;
 
 
 bool doit=true;
 while (doit)
  {
   
   // Switch based on one and only command line argument
   int solver_flag=atoi(CommandLineArgs::Argv[1]);
   switch (solver_flag)
    {
     
     
     // SuperLU
     //--------
    case 0:
     
     // This solver can handle non-symmetric matrices
     TanhSolnForAdvectionDiffusion::Peclet=200.0;
     
     std::cout << "Solving with SuperLU direct solver" 
               << std::endl;  
     break;
 
     
     // GMRES
     //------
    case 1:
     
     // This solver can handle non-symmetric matrices
     TanhSolnForAdvectionDiffusion::Peclet=200.0;
     
     // Build solver
     it_linear_solver_pt = new GMRES<CRDoubleMatrix>;
     
     // Build/pass preconditioner
     it_linear_solver_pt->preconditioner_pt() = prec_pt;
     
     // Switch on doc of convergence history and adjust tolerance
     it_linear_solver_pt->tolerance()=tol;
     it_linear_solver_pt->
      open_convergence_history_file_stream("RESLT/convergence.dat");
     
     // Set solver flag
     linear_solver_pt=it_linear_solver_pt;
     
     // doc
     std::cout << "Solving with GMRES" << std::endl;  
     
     break;
     
     
     
     
     // BiCGStab
     //---------
    case 2:
   
     // This solver can handle non-symmetric matrices
     TanhSolnForAdvectionDiffusion::Peclet=200.0;
 
     // Build solver
     it_linear_solver_pt = new BiCGStab<CRDoubleMatrix>;
   
     // Build/pass preconditioner
     it_linear_solver_pt->preconditioner_pt() = prec_pt;
   
     // Switch on doc of convergence history and adjust tolerance
     it_linear_solver_pt->tolerance()=tol;
     it_linear_solver_pt->
      open_convergence_history_file_stream("RESLT/convergence.dat");
   
     // Set solver flag
     linear_solver_pt=it_linear_solver_pt;
   
     // doc
     std::cout << "Solving with BiCGStab" << std::endl;  
  
     break;
 
 
 
     // CG
     //----
    case 3:
   
     // This solver can't handle non-symmetric matrices
     TanhSolnForAdvectionDiffusion::Peclet=0.0;
 
     // Build solver
     it_linear_solver_pt = new CG<CRDoubleMatrix>;
   
     // Build/pass preconditioner
     it_linear_solver_pt->preconditioner_pt() = prec_pt;
   
     // Switch on doc of convergence history and adjust tolerance
     it_linear_solver_pt->tolerance()=tol;
     it_linear_solver_pt->
      open_convergence_history_file_stream("RESLT/convergence.dat");
   
     // Set solver flag
     linear_solver_pt=it_linear_solver_pt;
   
     // doc
     std::cout << "Solving with CG" << std::endl;  
  
     break;
 
 
     // GS
     //----
    case 4:
      
     // This solver can't handle non-symmetric matrices
     TanhSolnForAdvectionDiffusion::Peclet=0.0;
 
     // Build solver 
     it_linear_solver_pt = new GS<CRDoubleMatrix>;
   
     // Note: Doesn't have a preconditioner
   
     // Switch on doc of convergence history and adjust tolerance
     it_linear_solver_pt->tolerance()=tol;
     it_linear_solver_pt->
      open_convergence_history_file_stream("RESLT/convergence.dat");
   
     // Set solver flag
     linear_solver_pt=it_linear_solver_pt;
   
     // doc
     std::cout << "Solving with GS" << std::endl;  
  
     break;
     
     // Damped Jacobi
     //--------------
    case 5:
      
     // This solver can't handle non-symmetric matrices
     TanhSolnForAdvectionDiffusion::Peclet=0.0;
 
     // Build solver 
     it_linear_solver_pt = new DampedJacobi<CRDoubleMatrix>;
   
     // Note: Doesn't have a preconditioner
   
     // Switch on doc of convergence history and adjust tolerance
     it_linear_solver_pt->tolerance()=tol;
     it_linear_solver_pt->
      open_convergence_history_file_stream("RESLT/convergence.dat");
   
     // Set solver flag
     linear_solver_pt=it_linear_solver_pt;
   
     // doc
     std::cout << "Solving with damped Jacobi" << std::endl;  
  
     break;
 
 
     // Wrong
     //------
    default:
 
     throw OomphLibError("Invalid command line argument: solver flag: ",
             OOMPH_CURRENT_FUNCTION,
             OOMPH_EXCEPTION_LOCATION);
     break;
   
    }
 
   // Build problem and solve it with specified linear solver
   run_it(linear_solver_pt);
 
   delete it_linear_solver_pt;
 
   // Comment this out if you want an endless loop to test
   // for memory leaks.
   doit=false;
 
  }
 
} //end of main

References oomph::CommandLineArgs::Argc, oomph::CommandLineArgs::Argv, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, TanhSolnForAdvectionDiffusion::Peclet, oomph::IterativeLinearSolver::preconditioner_pt(), run_it(), Flag_definition::setup(), and oomph::IterativeLinearSolver::tolerance().

run_it()

void run_it

(

LinearSolver *

linear_solver_pt

)

Run 2D AdvectionDiffusion problem with flux boundary conditions, using the iterative linear solver and preconditioner specified via the pointers.

{
 
 //Set up the problem
 //------------------
 
 //Set up the problem with 2D nine-node elements from the
 //RefineableQuadAdvectionDiffusionElement family. 
 //Pass pointer to source function. 
 TwoMeshFluxAdvectionDiffusionProblem<RefineableQAdvectionDiffusionElement<2,
  3> > problem(&TanhSolnForAdvectionDiffusion::source_function,
               &TanhSolnForAdvectionDiffusion::wind_function);
 
 // Set linear solver
 problem.linear_solver_pt()=linear_solver_pt;
 
 // Change max. number of iterations (some of the solvers used here
 // are rubbish!
 IterativeLinearSolver* iter_solver_pt=
  dynamic_cast<IterativeLinearSolver*>(linear_solver_pt);
 if (iter_solver_pt!=0)
  {
   iter_solver_pt->max_iter()=1000;
  }
 
 // Get Jacobian matrix and create dummy rhs for re-solve
 // (solution from re-solve should be a vector of ones)
 unsigned n_dof=problem.ndof();
 CRDoubleMatrix matrix;
 DoubleVector residual;
 DoubleVector other_rhs;
 DoubleVector one;
 DoubleVector solution;
 problem.get_jacobian(residual,matrix);
 one.build(residual.distribution_pt(),0.0);
 other_rhs.build(residual.distribution_pt(),0.0);
 one.initialise(1.0);
 matrix.multiply(one,other_rhs);
 linear_solver_pt->enable_resolve();
 
 
 // Normal solve
 //-------------
 
 // Create label for output
 DocInfo doc_info;
 
 // Set output directory
 doc_info.set_directory("RESLT");
 
 // Solve the problem
 problem.newton_solve();
 
 //Output solution
 problem.doc_solution(doc_info);
 
 
 // Re-solve
 //---------
 
 // Open convergence history for resolve
 IterativeLinearSolver* it_solver_pt=
  dynamic_cast<IterativeLinearSolver*>(linear_solver_pt);
 if (it_solver_pt!=0)
  {
   it_solver_pt->
    open_convergence_history_file_stream("RESLT/resolve_convergence.dat");
  }
 
 // Resolve
 linear_solver_pt->resolve(other_rhs,solution);
 
 // Get error
 double error=0.0;
 for (unsigned i=0;i<n_dof;i++)
  {
   error+=pow((solution[i]-1.0),2);
  }
 
 // Doc
 ofstream error_file("RESLT/resolve_error.dat");
 error_file << "Error for resolve: " 
            << sqrt(error/double(n_dof)) 
            << std::endl;
 std::cout  << "Error for resolve: " 
            << sqrt(error/double(n_dof)) 
            << std::endl;
 error_file.close();
 
 // Switch of ability to resolve
 linear_solver_pt->disable_resolve();
 
 
 
 // Linear-algebra solve
 //---------------------
 
 if (it_solver_pt!=0)
  {
   it_solver_pt->
    open_convergence_history_file_stream("RESLT/la_solve_convergence.dat");
  }
 
 DoubleVector solution2;
 linear_solver_pt->solve(&matrix,other_rhs,solution2);
 
 // Get error
 error=0.0;
 for (unsigned i=0;i<n_dof;i++)
  {
   error+=pow((solution2[i]-1.0),2);
  }
 
 // Doc
 ofstream error_file2("RESLT/la_solve_error.dat");
 error_file2 << "Error for LA solve: " 
             << sqrt(error/double(n_dof)) 
             << std::endl;
 std::cout << "Error for LA solve: " 
           << sqrt(error/double(n_dof)) 
           << std::endl;
 error_file2.close();
 
} //end of run_it

References oomph::DoubleVector::build(), oomph::LinearSolver::disable_resolve(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::enable_resolve(), calibrate::error, i, oomph::DoubleVector::initialise(), matrix(), oomph::IterativeLinearSolver::max_iter(), Eigen::ArrayBase< Derived >::pow(), problem, oomph::LinearSolver::resolve(), oomph::DocInfo::set_directory(), BiharmonicTestFunctions1::solution(), oomph::LinearSolver::solve(), TanhSolnForAdvectionDiffusion::source_function(), sqrt(), and TanhSolnForAdvectionDiffusion::wind_function().

Referenced by main().