#include "generic.h"
#include "poisson.h"
#include "meshes/simple_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::get_source (const Vector< double > &x, double &source)
	Source function to make it an exact solution. More...

void	run (const string &dir_name, LinearSolver *linear_solver_pt, const unsigned nel_1d, bool mess_up_order)

int	main ()
	Driver code for 2D Poisson problem – compare different solvers. More...

Function Documentation

◆ main()

int main ( )

Driver code for 2D Poisson problem – compare different solvers.

Use SuperLU with compressed row storage

Use compressed row storage

Switch on full doc

Use SuperLU with compressed column storage

Use compressed row storage

Switch on full doc

Use dense matrix LU decomposition

 {
  
  // Pointer to linear solver
  LinearSolver* linear_solver_pt;
  
  // Result directory
  string dir_name;
  
  // Cpu start/end times
  clock_t cpu_start,cpu_end;
  
  // Storage for cpu times with and without messed up order
  Vector<double> superlu_cr_cpu(2);
  Vector<double> superlu_cc_cpu(2);
  Vector<double> hsl_ma42_cpu(2);
  Vector<double> hsl_ma42_reordered_cpu(2);
  Vector<double> dense_lu_cpu(2);
  Vector<double> fd_lu_cpu(2);
  
  // Flag to indicate if order should be messed up:
  bool mess_up_order;
  
  // Number of elements along 1D edge of mesh; total number of elements
  // is nel_1d x nel_1d
  unsigned nel_1d;
  
  // Do run with and without messed up elements
  for (unsigned mess=0;mess<2;mess++)
   {
  
    // Mess up order?
    if (mess==0)
     {
      mess_up_order=true;
     }
    else
     {
      mess_up_order=false;
     }
  
    //-----------------------------------------
    
    cout << std::endl;
    cout << " Use SuperLU with compressed row storage: " << std::endl;
    cout << "========================================= " << std::endl;
    cout << std::endl;
  
    // Start cpu clock
    cpu_start=clock();
    
    // Build linear solver
    linear_solver_pt = new SuperLUSolver;
    
    static_cast<SuperLUSolver*>(linear_solver_pt)
     ->use_compressed_row_for_superlu_serial();
    
    static_cast<SuperLUSolver*>(linear_solver_pt)->enable_doc_stats();
    
    // Choose result directory
    dir_name="RESLT_cr";
  
    // Number of elements along 1D edge of mesh; total number of elements
    // is nel_1d x nel_1d
    nel_1d=20;
  
    // Run it
    run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
    
    // Note: solver does not have to be deleted -- it's killed automatically
    // in the problem destructor.
    
    // End cpu clock
    cpu_end=clock();
    
    // Timing
    superlu_cr_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
    
    
    
    
    //--------------------------------------------
  
  
    cout << std::endl;
    cout << " Use SuperLU with compressed column storage: " << std::endl;
    cout << "============================================ " << std::endl;
    cout << std::endl;
  
    
    // Start cpu clock
    cpu_start=clock();
    
    // Build linear solver
    linear_solver_pt = new SuperLUSolver;
    
    static_cast<SuperLUSolver*>(linear_solver_pt)
     ->use_compressed_column_for_superlu_serial();
    
    static_cast<SuperLUSolver*>(linear_solver_pt)->enable_doc_stats();
  
    // Choose result directory
    dir_name="RESLT_cc";
  
    // Number of elements along 1D edge of mesh; total number of elements
    // is nel_1d x nel_1d
    nel_1d=20;
  
    // Run it
    run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
    
    // Note: solver does not have to be deleted -- it's killed automatically
    // in the problem destructor.
    
    // End cpu clock
    cpu_end=clock();
    
    // Timing
    superlu_cc_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
    
    
 #ifdef HAVE_HSL_SOURCES
    
    //---------------------------------------------------------
    
    cout << std::endl;
    cout << " Use HSL frontal solver MA42 without element re-ordering: " << std::endl;
    cout << "========================================================= " << std::endl;
    cout << std::endl;
  
    // Start cpu clock
    cpu_start=clock();
    
    // Build linear solver
    linear_solver_pt = new HSL_MA42;
    
    static_cast<HSL_MA42*>(linear_solver_pt)->enable_doc_stats();
    
    // Choose result directory
    dir_name="RESLT_frontal";
  
  
    // Number of elements along 1D edge of mesh; total number of elements
    // is nel_1d x nel_1d
    nel_1d=20;
  
    // Run it
    run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
    
    // Note: solver does not have to be deleted -- it's killed automatically
    // in the problem destructor.
    
    
    // End cpu clock
    cpu_end=clock();
    
    // Timing
    hsl_ma42_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
    
    
    
    //--------------------------------------------------------------------------
    
    cout << std::endl;
    cout << " Use HSL frontal solver MA42 with element re-ordering: " << std::endl;
    cout << "====================================================== " << std::endl;
    cout << std::endl;
  
    // Start cpu clock
    cpu_start=clock();
    
    // Build linear solver
    linear_solver_pt = new HSL_MA42;
    
    static_cast<HSL_MA42*>(linear_solver_pt)->enable_doc_stats();
    
    
    static_cast<HSL_MA42*>(linear_solver_pt)->enable_reordering();
    
    // Choose result directory
    dir_name="RESLT_frontal_reordered";
  
  
    // Number of elements along 1D edge of mesh; total number of elements
    // is nel_1d x nel_1d
    nel_1d=20;
  
    // Run it
    run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
    
    // Note: solver does not have to be deleted -- it's killed automatically
    // in the problem destructor.
    
    
    // End cpu clock
    cpu_end=clock();
    
    // Timing
    hsl_ma42_reordered_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
    
  
 #endif
  
    //-----------------------------------
    
    cout << std::endl;
    cout << " Use dense matrix LU decomposition: " << std::endl;
    cout << "=================================== " << std::endl;
    cout << std::endl;
  
    // Start cpu clock
    cpu_start=clock();
    
    // Build linear solver
    linear_solver_pt = new DenseLU;
    
    // Choose result directory
    dir_name="RESLT_dense_LU";
  
    // Number of elements along 1D edge of mesh; total number of elements
    // is nel_1d x nel_1d
    nel_1d=4;
  
    // Run it
    run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
    
    // Note: solver does not have to be deleted -- it's killed automatically
    // in the problem destructor.
    
    
    // End cpu clock
    cpu_end=clock();
    
    // Timing
    dense_lu_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
    
  
  
  
    //-----------------------------------
    
    cout << std::endl;
    cout << " Use dense FD-ed matrix LU decomposition: " << std::endl;
    cout << "========================================= " << std::endl;
    cout << std::endl;
  
    // Start cpu clock
    cpu_start=clock();
    
    // Build linear solver
    linear_solver_pt = new FD_LU;
    
    // Choose result directory
    dir_name="RESLT_FD_LU";
  
    // Number of elements along 1D edge of mesh; total number of elements
    // is nel_1d x nel_1d
    nel_1d=4;
  
    // Run it
    run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
    
    // Note: solver does not have to be deleted -- it's killed automatically
    // in the problem destructor.
       
    // End cpu clock
    cpu_end=clock();
    
    // Timing
    fd_lu_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
    
   }
  
  
  // Doc timings with and without messed up elements
  for (unsigned mess=0;mess<2;mess++)
   {
    cout << std::endl   << std::endl   << std::endl ;
    if (mess==0)
     {
      cout << "TIMINGS WITH MESSED UP ORDERING OF ELEMENTS: " << std::endl;
      cout << "============================================ " << std::endl;
  
     }
    else
     {
      cout << "TIMINGS WITHOUT MESSED UP ORDERING OF ELEMENTS: " << std::endl;
      cout << "=============================================== " << std::endl;
     }
  
    cout << "CPU time with SuperLU compressed row                 : " 
         <<  superlu_cr_cpu[mess] << std::endl;
    cout << "CPU time with SuperLU compressed col                 : " 
         <<  superlu_cc_cpu[mess] << std::endl;
 #ifdef HAVE_HSL_SOURCES
    cout << "CPU time with MA42 frontal solver                    : " 
         <<  hsl_ma42_cpu[mess] << std::endl;
    cout << "CPU time with MA42 frontal solver (incl. reordering) : " 
         <<  hsl_ma42_reordered_cpu[mess] << std::endl;
 #endif
    cout << "CPU time with dense LU solver (fewer elements!)      : " 
         <<  dense_lu_cpu[mess] << std::endl;  
    cout << "CPU time with dense LU solver & FD (fewer elements!) : " 
         <<  fd_lu_cpu[mess] << std::endl;  
    cout << std::endl   << std::endl   << std::endl ;
   }
  
  
  
  
 } //end of main

References run().

◆ run()

void run	(	const string &	dir_name,
		LinearSolver *	linear_solver_pt,
		const unsigned	nel_1d,
		bool	mess_up_order
	)

Build and run problem with specified linear solver. Also pass flag to specify if order of elements should be messed up. nel_1d is the number of elements along the 1D mesh edge. Total number of elements is nel_1d x nel_1d.

Set linear solver:

 {
  
  //Set up the problem
  //------------------
  
  // Create the problem with 2D nine-node elements from the
  // QPoissonElement family. Pass pointer to source function
  // and flag to specify if order of elements should be messed up.
  PoissonProblem<QPoissonElement<2,3> > 
   problem(&TanhSolnForPoisson::get_source,nel_1d,mess_up_order);
  
  
  //--------------------
  problem.linear_solver_pt() = linear_solver_pt;
  
  
  // Create label for output
  //------------------------
  DocInfo doc_info;
  
  // Set output directory
  doc_info.set_directory(dir_name);
  
  // Step number
  doc_info.number()=0;
  
  // Check if we're ready to go:
  //----------------------------
  cout << "\n\n\nProblem self-test:";
  if (problem.self_test()==0) 
   {
    cout << "passed: Problem can be solved." << std::endl;
   }
  else 
   {
    throw OomphLibError("Self test failed",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  
  
  // 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; 
  
  // Solve the problem
  problem.newton_solve();
  
  //Output the solution
  problem.doc_solution(doc_info);
  
  
 } //end of run

References TanhSolnForPoisson::Alpha, TanhSolnForPoisson::get_source(), oomph::DocInfo::number(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, problem, oomph::DocInfo::set_directory(), and TanhSolnForPoisson::TanPhi.

Classes

Namespaces

Functions

Function Documentation

◆ main()

◆ run()