two_d_poisson_compare_solvers.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 //Driver for a simple 2D poisson problem -- compare various linear solvers
27 
28 //Generic routines
29 #include "generic.h"
30 
31 // The Poisson equations
32 #include "poisson.h"
33 
34 // The mesh
35 #include "meshes/simple_rectangular_quadmesh.h"
36 
37 using namespace std;
38 
39 using namespace oomph;
40 
41 //===== start_of_namespace=============================================
42 /// Namespace for exact solution for Poisson equation with "sharp step"
43 //=====================================================================
44 namespace TanhSolnForPoisson
45 {
46 
47  /// Parameter for steepness of "step"
48  double Alpha=1.0;
49 
50  /// Parameter for angle Phi of "step"
51  double TanPhi=0.0;
52 
53  /// Exact solution as a Vector
54  void get_exact_u(const Vector<double>& x, Vector<double>& u)
55  {
56  u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
57  }
58 
59  /// Source function required to make the solution above an exact solution
60  void get_source(const Vector<double>& x, double& source)
61  {
62  source = 2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
63  (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*
64  Alpha*Alpha*TanPhi*TanPhi+2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
65  (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*Alpha;
66  }
67 
68 } // end of namespace
69 
70 
71 
72 
73 
74 
75 
76 
77 //====== start_of_problem_class=======================================
78 /// 2D Poisson problem on rectangular domain, discretised with
79 /// 2D QPoisson elements. The specific type of element is
80 /// specified via the template parameter.
81 //====================================================================
82 template<class ELEMENT>
83 class PoissonProblem : public Problem
84 {
85 
86 public:
87 
88  /// Constructor: Pass pointer to source function and flag to indicate
89  /// if order of elements should be messed up. nel_1d is the number of
90  /// elements along the 1d mesh edges -- total number of elements is
91  /// nel_1d x nel_1d.
92  PoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt,
93  const unsigned& nel_1d,
94  const bool& mess_up_order);
95 
96  /// Destructor (empty)
98 
99  /// Update the problem specs before solve: Reset boundary conditions
100  /// to the values from the exact solution.
102 
103  /// Update the problem after solve (empty)
105 
106  /// Doc the solution. DocInfo object stores flags/labels for where the
107  /// output gets written to
108  void doc_solution(DocInfo& doc_info);
109 
110 private:
111 
112  /// Pointer to source function
113  PoissonEquations<2>::PoissonSourceFctPt Source_fct_pt;
114 
115 }; // end of problem class
116 
117 
118 
119 
120 //=====start_of_constructor===============================================
121 /// Constructor for Poisson problem: Pass pointer to source function
122 /// and a flag that specifies if the order of the elements should
123 /// be messed up. nel_1d is the number of elements along
124 /// the 1d mesh edges -- total number of elements is nel_1d x nel_1d.
125 //========================================================================
126 template<class ELEMENT>
128 PoissonProblem(PoissonEquations<2>::PoissonSourceFctPt source_fct_pt,
129  const unsigned& nel_1d, const bool& mess_up_order)
130  : Source_fct_pt(source_fct_pt)
131 {
132 
133  // Setup mesh
134 
135  // # of elements in x-direction
136  unsigned n_x=nel_1d;
137 
138  // # of elements in y-direction
139  unsigned n_y=nel_1d;
140 
141  // Domain length in x-direction
142  double l_x=1.0;
143 
144  // Domain length in y-direction
145  double l_y=2.0;
146 
147  // Build and assign mesh
148  Problem::mesh_pt() = new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
149 
150 
151 
152  // Mess up order of elements in the Problem's mesh:
153  //-------------------------------------------------
154  if (mess_up_order)
155  {
156  unsigned n_element=mesh_pt()->nelement();
157 
158  // Make copy
159  Vector<ELEMENT*> tmp_element_pt(n_element);
160  for (unsigned e=0;e<n_element;e++)
161  {
162  tmp_element_pt[e]=dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
163  }
164 
165  // Reorder
166  unsigned e_half=unsigned(0.5*double(n_element));
167  unsigned e_lo=e_half-3;
168  unsigned e_up=n_element-e_lo;
169  unsigned count=0;
170  for (unsigned e=0;e<e_lo;e++)
171  {
172  mesh_pt()->element_pt(count)=tmp_element_pt[e];
173  count++;
174  mesh_pt()->element_pt(count)=tmp_element_pt[n_element-e-1];
175  count++;
176  }
177  for (unsigned e=e_lo;e<e_up;e++)
178  {
179  mesh_pt()->element_pt(count)=tmp_element_pt[e];
180  count++;
181  }
182  }
183 
184  // Set the boundary conditions for this problem: All nodes are
185  // free by default -- only need to pin the ones that have Dirichlet conditions
186  // here.
187  unsigned num_bound = mesh_pt()->nboundary();
188  for(unsigned ibound=0;ibound<num_bound;ibound++)
189  {
190  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
191  for (unsigned inod=0;inod<num_nod;inod++)
192  {
193  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
194  }
195  }
196 
197  // Complete the build of all elements so they are fully functional
198 
199  // Loop over the elements to set up element-specific
200  // things that cannot be handled by the (argument-free!) ELEMENT
201  // constructor: Pass pointer to source function
202  unsigned n_element = mesh_pt()->nelement();
203  for(unsigned i=0;i<n_element;i++)
204  {
205  // Upcast from GeneralsedElement to the present element
206  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
207 
208  //Set the source function pointer
209  el_pt->source_fct_pt() = Source_fct_pt;
210  }
211 
212 
213  // Setup equation numbering scheme
214  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
215 
216 } // end of constructor
217 
218 
219 
220 
221 //========================================start_of_actions_before_newton_solve===
222 /// Update the problem specs before solve: (Re-)set boundary conditions
223 /// to the values from the exact solution.
224 //========================================================================
225 template<class ELEMENT>
227 {
228  // How many boundaries are there?
229  unsigned num_bound = mesh_pt()->nboundary();
230 
231  //Loop over the boundaries
232  for(unsigned ibound=0;ibound<num_bound;ibound++)
233  {
234  // How many nodes are there on this boundary?
235  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
236 
237  // Loop over the nodes on boundary
238  for (unsigned inod=0;inod<num_nod;inod++)
239  {
240  // Get pointer to node
241  Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
242 
243  // Extract nodal coordinates from node:
244  Vector<double> x(2);
245  x[0]=nod_pt->x(0);
246  x[1]=nod_pt->x(1);
247 
248  // Compute the value of the exact solution at the nodal point
249  Vector<double> u(1);
251 
252  // Assign the value to the one (and only) nodal value at this node
253  nod_pt->set_value(0,u[0]);
254  }
255  }
256 } // end of actions before solve
257 
258 
259 
260 //===============start_of_doc=============================================
261 /// Doc the solution: doc_info contains labels/output directory etc.
262 //========================================================================
263 template<class ELEMENT>
264 void PoissonProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
265 {
266 
267  ofstream some_file;
268  char filename[100];
269 
270  // Number of plot points: npts x npts
271  unsigned npts=5;
272 
273  // Output solution
274  //-----------------
275  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
276  doc_info.number());
277  some_file.open(filename);
278  mesh_pt()->output(some_file,npts);
279  some_file.close();
280 
281 
282  // Output exact solution
283  //----------------------
284  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
285  doc_info.number());
286  some_file.open(filename);
287  mesh_pt()->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
288  some_file.close();
289 
290  // Doc error and return of the square of the L2 error
291  //---------------------------------------------------
292  double error,norm;
293  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
294  doc_info.number());
295  some_file.open(filename);
296  mesh_pt()->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
297  error,norm);
298  some_file.close();
299 
300  // Doc L2 error and norm of solution
301  cout << "\nNorm of error : " << sqrt(error) << std::endl;
302  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
303 
304 } // end of doc
305 
306 
307 
308 
309 
310 
311 //===== start_of_run======================================================
312 /// Build and run problem with specified linear solver. Also pass
313 /// flag to specify if order of elements should be messed up.
314 /// nel_1d is the number of elements along the 1D mesh edge.
315 /// Total number of elements is nel_1d x nel_1d.
316 //========================================================================
317 void run(const string& dir_name,
318  LinearSolver* linear_solver_pt,
319  const unsigned nel_1d,
320  bool mess_up_order)
321 {
322 
323  //Set up the problem
324  //------------------
325 
326  // Create the problem with 2D nine-node elements from the
327  // QPoissonElement family. Pass pointer to source function
328  // and flag to specify if order of elements should be messed up.
330  problem(&TanhSolnForPoisson::get_source,nel_1d,mess_up_order);
331 
332 
333  /// Set linear solver:
334  //--------------------
335  problem.linear_solver_pt() = linear_solver_pt;
336 
337 
338  // Create label for output
339  //------------------------
340  DocInfo doc_info;
341 
342  // Set output directory
343  doc_info.set_directory(dir_name);
344 
345  // Step number
346  doc_info.number()=0;
347 
348  // Check if we're ready to go:
349  //----------------------------
350  cout << "\n\n\nProblem self-test:";
351  if (problem.self_test()==0)
352  {
353  cout << "passed: Problem can be solved." << std::endl;
354  }
355  else
356  {
357  throw OomphLibError("Self test failed",
358  OOMPH_CURRENT_FUNCTION,
359  OOMPH_EXCEPTION_LOCATION);
360  }
361 
362 
363  // Set the orientation of the "step" to 45 degrees
365 
366  // Initial value for the steepness of the "step"
368 
369  // Solve the problem
370  problem.newton_solve();
371 
372  //Output the solution
373  problem.doc_solution(doc_info);
374 
375 
376 } //end of run
377 
378 
379 
380 
381 
382 //===== start_of_main=====================================================
383 /// Driver code for 2D Poisson problem -- compare different solvers
384 //========================================================================
385 int main()
386 {
387 
388  // Pointer to linear solver
389  LinearSolver* linear_solver_pt;
390 
391  // Result directory
392  string dir_name;
393 
394  // Cpu start/end times
395  clock_t cpu_start,cpu_end;
396 
397  // Storage for cpu times with and without messed up order
398  Vector<double> superlu_cr_cpu(2);
399  Vector<double> superlu_cc_cpu(2);
400  Vector<double> hsl_ma42_cpu(2);
401  Vector<double> hsl_ma42_reordered_cpu(2);
402  Vector<double> dense_lu_cpu(2);
403  Vector<double> fd_lu_cpu(2);
404 
405  // Flag to indicate if order should be messed up:
406  bool mess_up_order;
407 
408  // Number of elements along 1D edge of mesh; total number of elements
409  // is nel_1d x nel_1d
410  unsigned nel_1d;
411 
412  // Do run with and without messed up elements
413  for (unsigned mess=0;mess<2;mess++)
414  {
415 
416  // Mess up order?
417  if (mess==0)
418  {
419  mess_up_order=true;
420  }
421  else
422  {
423  mess_up_order=false;
424  }
425 
426  /// Use SuperLU with compressed row storage
427  //-----------------------------------------
428 
429  cout << std::endl;
430  cout << " Use SuperLU with compressed row storage: " << std::endl;
431  cout << "========================================= " << std::endl;
432  cout << std::endl;
433 
434  // Start cpu clock
435  cpu_start=clock();
436 
437  // Build linear solver
438  linear_solver_pt = new SuperLUSolver;
439 
440  /// Use compressed row storage
441  static_cast<SuperLUSolver*>(linear_solver_pt)
442  ->use_compressed_row_for_superlu_serial();
443 
444  /// Switch on full doc
445  static_cast<SuperLUSolver*>(linear_solver_pt)->enable_doc_stats();
446 
447  // Choose result directory
448  dir_name="RESLT_cr";
449 
450  // Number of elements along 1D edge of mesh; total number of elements
451  // is nel_1d x nel_1d
452  nel_1d=20;
453 
454  // Run it
455  run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
456 
457  // Note: solver does not have to be deleted -- it's killed automatically
458  // in the problem destructor.
459 
460  // End cpu clock
461  cpu_end=clock();
462 
463  // Timing
464  superlu_cr_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
465 
466 
467 
468 
469  /// Use SuperLU with compressed column storage
470  //--------------------------------------------
471 
472 
473  cout << std::endl;
474  cout << " Use SuperLU with compressed column storage: " << std::endl;
475  cout << "============================================ " << std::endl;
476  cout << std::endl;
477 
478 
479  // Start cpu clock
480  cpu_start=clock();
481 
482  // Build linear solver
483  linear_solver_pt = new SuperLUSolver;
484 
485  /// Use compressed row storage
486  static_cast<SuperLUSolver*>(linear_solver_pt)
487  ->use_compressed_column_for_superlu_serial();
488 
489  /// Switch on full doc
490  static_cast<SuperLUSolver*>(linear_solver_pt)->enable_doc_stats();
491 
492  // Choose result directory
493  dir_name="RESLT_cc";
494 
495  // Number of elements along 1D edge of mesh; total number of elements
496  // is nel_1d x nel_1d
497  nel_1d=20;
498 
499  // Run it
500  run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
501 
502  // Note: solver does not have to be deleted -- it's killed automatically
503  // in the problem destructor.
504 
505  // End cpu clock
506  cpu_end=clock();
507 
508  // Timing
509  superlu_cc_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
510 
511 
512 #ifdef HAVE_HSL_SOURCES
513 
514  /// Use HSL frontal solver MA42 without element re-ordering
515  //---------------------------------------------------------
516 
517  cout << std::endl;
518  cout << " Use HSL frontal solver MA42 without element re-ordering: " << std::endl;
519  cout << "========================================================= " << std::endl;
520  cout << std::endl;
521 
522  // Start cpu clock
523  cpu_start=clock();
524 
525  // Build linear solver
526  linear_solver_pt = new HSL_MA42;
527 
528  /// Switch on full doc
529  static_cast<HSL_MA42*>(linear_solver_pt)->enable_doc_stats();
530 
531  // Choose result directory
532  dir_name="RESLT_frontal";
533 
534 
535  // Number of elements along 1D edge of mesh; total number of elements
536  // is nel_1d x nel_1d
537  nel_1d=20;
538 
539  // Run it
540  run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
541 
542  // Note: solver does not have to be deleted -- it's killed automatically
543  // in the problem destructor.
544 
545 
546  // End cpu clock
547  cpu_end=clock();
548 
549  // Timing
550  hsl_ma42_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
551 
552 
553 
554  /// Use HSL frontal solver MA42 with element reordering by Sloan's algorithm
555  //--------------------------------------------------------------------------
556 
557  cout << std::endl;
558  cout << " Use HSL frontal solver MA42 with element re-ordering: " << std::endl;
559  cout << "====================================================== " << std::endl;
560  cout << std::endl;
561 
562  // Start cpu clock
563  cpu_start=clock();
564 
565  // Build linear solver
566  linear_solver_pt = new HSL_MA42;
567 
568  /// Switch on full doc
569  static_cast<HSL_MA42*>(linear_solver_pt)->enable_doc_stats();
570 
571 
572  /// Switch on re-ordering
573  static_cast<HSL_MA42*>(linear_solver_pt)->enable_reordering();
574 
575  // Choose result directory
576  dir_name="RESLT_frontal_reordered";
577 
578 
579  // Number of elements along 1D edge of mesh; total number of elements
580  // is nel_1d x nel_1d
581  nel_1d=20;
582 
583  // Run it
584  run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
585 
586  // Note: solver does not have to be deleted -- it's killed automatically
587  // in the problem destructor.
588 
589 
590  // End cpu clock
591  cpu_end=clock();
592 
593  // Timing
594  hsl_ma42_reordered_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
595 
596 
597 #endif
598 
599  /// Use dense matrix LU decomposition
600  //-----------------------------------
601 
602  cout << std::endl;
603  cout << " Use dense matrix LU decomposition: " << std::endl;
604  cout << "=================================== " << std::endl;
605  cout << std::endl;
606 
607  // Start cpu clock
608  cpu_start=clock();
609 
610  // Build linear solver
611  linear_solver_pt = new DenseLU;
612 
613  // Choose result directory
614  dir_name="RESLT_dense_LU";
615 
616  // Number of elements along 1D edge of mesh; total number of elements
617  // is nel_1d x nel_1d
618  nel_1d=4;
619 
620  // Run it
621  run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
622 
623  // Note: solver does not have to be deleted -- it's killed automatically
624  // in the problem destructor.
625 
626 
627  // End cpu clock
628  cpu_end=clock();
629 
630  // Timing
631  dense_lu_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
632 
633 
634 
635 
636  /// Use dense matrix LU decomposition
637  //-----------------------------------
638 
639  cout << std::endl;
640  cout << " Use dense FD-ed matrix LU decomposition: " << std::endl;
641  cout << "========================================= " << std::endl;
642  cout << std::endl;
643 
644  // Start cpu clock
645  cpu_start=clock();
646 
647  // Build linear solver
648  linear_solver_pt = new FD_LU;
649 
650  // Choose result directory
651  dir_name="RESLT_FD_LU";
652 
653  // Number of elements along 1D edge of mesh; total number of elements
654  // is nel_1d x nel_1d
655  nel_1d=4;
656 
657  // Run it
658  run(dir_name,linear_solver_pt,nel_1d,mess_up_order);
659 
660  // Note: solver does not have to be deleted -- it's killed automatically
661  // in the problem destructor.
662 
663  // End cpu clock
664  cpu_end=clock();
665 
666  // Timing
667  fd_lu_cpu[mess]=double(cpu_end-cpu_start)/CLOCKS_PER_SEC;
668 
669  }
670 
671 
672  // Doc timings with and without messed up elements
673  for (unsigned mess=0;mess<2;mess++)
674  {
675  cout << std::endl << std::endl << std::endl ;
676  if (mess==0)
677  {
678  cout << "TIMINGS WITH MESSED UP ORDERING OF ELEMENTS: " << std::endl;
679  cout << "============================================ " << std::endl;
680 
681  }
682  else
683  {
684  cout << "TIMINGS WITHOUT MESSED UP ORDERING OF ELEMENTS: " << std::endl;
685  cout << "=============================================== " << std::endl;
686  }
687 
688  cout << "CPU time with SuperLU compressed row : "
689  << superlu_cr_cpu[mess] << std::endl;
690  cout << "CPU time with SuperLU compressed col : "
691  << superlu_cc_cpu[mess] << std::endl;
692 #ifdef HAVE_HSL_SOURCES
693  cout << "CPU time with MA42 frontal solver : "
694  << hsl_ma42_cpu[mess] << std::endl;
695  cout << "CPU time with MA42 frontal solver (incl. reordering) : "
696  << hsl_ma42_reordered_cpu[mess] << std::endl;
697 #endif
698  cout << "CPU time with dense LU solver (fewer elements!) : "
699  << dense_lu_cpu[mess] << std::endl;
700  cout << "CPU time with dense LU solver & FD (fewer elements!) : "
701  << fd_lu_cpu[mess] << std::endl;
702  cout << std::endl << std::endl << std::endl ;
703  }
704 
705 
706 
707 
708 } //end of main
709 
710 
711 
712 
713 
2D Poisson problem on rectangular domain, discretised with 2D QPoisson elements. The specific type of...
PoissonEquations< 2 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the exact solutio...
void actions_after_newton_solve()
Update the problem after solve (empty)
PoissonProblem(PoissonEquations< 2 >::PoissonSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function.
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to.
~PoissonProblem()
Destructor (empty)
Namespace for exact solution for Poisson equation with "sharp step".
double TanPhi
Parameter for angle Phi of "step".
void get_source(const Vector< double > &x, double &source)
Source function required to make the solution above an exact solution.
double Alpha
Parameter for steepness of "step".
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
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 sh...
int main()
Driver code for 2D Poisson problem – compare different solvers.