complex_harmonic.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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Driver to solve the harmonic equation as a quadratic eigenvalue problem,
31 //by solving two first-order systems with homogeneous Dirichlet boundary
32 //conditions. The sign is chosen to that all the eigenvalues are imaginary
33 //when Mu is small (Mu=0 actually leads to a singular mass matrix so should
34 //be avoided), but pairs merge to become real as Mu increases.
35 
36 // Generic oomph-lib routines
37 #include "generic.h"
38 
39 // Include the mesh
40 #include "meshes/one_d_mesh.h"
41 
42 using namespace std;
43 
44 using namespace oomph;
45 
46 ///Namespace for the shift applied to the eigenproblem
47 namespace EigenproblemShift
48 {
49  //Parameter chosen to that only one complex-conjugate pair has merged
50  double Mu = 6.5; //Will lead to singular mass matrix if set to zero
51 }
52 
53 
54 //===================================================================
55 /// Function-type-object to perform comparison of complex data types
56 /// Needed to sort the complex eigenvalues into order based on the
57 /// size of the real part.
58 //==================================================================
59 template <class T>
61 {
62 public:
63 
64  /// Comparison in terms of magnitude of complex number
65  bool operator()(const complex<T> &x, const complex<T> &y) const
66  {
67  //Return the order in terms of magnitude if they are not equal
68  //Include a tolerance to avoid processor specific ordering
69  if(std::abs(std::abs(x) - std::abs(y)) > 1.0e-10 )
70  {return std::abs(x) < std::abs(y);}
71  //Otherwise sort on real part, again with a tolerance
72  else if(std::abs(x.real() - y.real()) > 1.0e-10)
73  {return x.real() < y.real();}
74  //Otherwise sort on imaginary part
75  else
76  {return x.imag() < y.imag();}
77  }
78 };
79 
80 
81 //=================================================================
82 /// A class for all elements that solve the eigenvalue problem
83 /// \f[
84 /// \frac{\partial w}{\partial x} = \lambda u
85 /// \f]
86 /// \f[
87 /// \frac{\partial u}{\partial x} = (\lambda - \mu) w
88 /// \f]
89 /// This class contains the generic maths. Shape functions, geometric
90 /// mapping etc. must get implemented in derived class.
91 //================================================================
92 class ComplexHarmonicEquations : public virtual FiniteElement
93 {
94 
95 public:
96  /// Empty Constructor
98 
99  /// \short Access function: First eigenfunction value at local node n
100  /// Note that solving the eigenproblem does not assign values
101  /// to this storage space. It is used for output purposes only.
102  virtual inline double u(const unsigned& n) const
103  {return nodal_value(n,0);}
104 
105  /// \short Second eigenfunction value at local node n
106  virtual inline double w(const unsigned& n) const
107  {return nodal_value(n,1);}
108 
109  /// Output the eigenfunction with default number of plot points
110  void output(ostream &outfile)
111  {
112  unsigned nplot=5;
113  output(outfile,nplot);
114  }
115 
116  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
117  /// Nplot plot points
118  void output(ostream &outfile, const unsigned &nplot)
119  {
120  //Vector of local coordinates
121  Vector<double> s(1);
122 
123  // Tecplot header info
124  outfile << tecplot_zone_string(nplot);
125 
126  // Loop over plot points
127  unsigned num_plot_points=nplot_points(nplot);
128  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
129  {
130  // Get local coordinates of plot point
131  get_s_plot(iplot,nplot,s);
132  //Output the coordinate and the eigenfunction
133  outfile << interpolated_x(s,0) << " " << interpolated_u(s)
134  << " " << interpolated_w(s) << std::endl;
135  }
136  // Write tecplot footer (e.g. FE connectivity lists)
137  write_tecplot_zone_footer(outfile,nplot);
138  }
139 
140  /// \short Assemble the contributions to the jacobian and mass
141  /// matrices
143  Vector<double> &residuals,
144  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix)
145  {
146  //Find out how many nodes there are
147  unsigned n_node = nnode();
148 
149  //Set up memory for the shape functions and their derivatives
150  Shape psi(n_node);
151  DShape dpsidx(n_node,1);
152 
153  //Set the number of integration points
154  unsigned n_intpt = integral_pt()->nweight();
155 
156  //Integers to store the local equation and unknown numbers
157  int local_eqn=0, local_unknown=0;
158 
159  //Loop over the integration points
160  for(unsigned ipt=0;ipt<n_intpt;ipt++)
161  {
162  //Get the integral weight
163  double w = integral_pt()->weight(ipt);
164 
165  //Call the derivatives of the shape and test functions
166  double J = dshape_eulerian_at_knot(ipt,psi,dpsidx);
167 
168  //Premultiply the weights and the Jacobian
169  double W = w*J;
170 
171  //Assemble the contributions to the mass matrix
172  //Loop over the test functions
173  for(unsigned l=0;l<n_node;l++)
174  {
175  //Get the local equation number
176  local_eqn = u_local_eqn(l,0);
177  //If it's not a boundary condition
178  if(local_eqn >= 0)
179  {
180  //Loop over the shape functions
181  for(unsigned l2=0;l2<n_node;l2++)
182  {
183  local_unknown = u_local_eqn(l2,0);
184  //If at a non-zero degree of freedom add in the entry
185  if(local_unknown >= 0)
186  {
187  //This corresponds to the top left B block
188  mass_matrix(local_eqn, local_unknown) += psi(l2)*psi(l)*W;
189  }
190  local_unknown = u_local_eqn(l2,1);
191  //If at a non-zero degree of freedom add in the entry
192  if(local_unknown >= 0)
193  {
194  //This corresponds to the top right A block
195  jacobian(local_eqn,local_unknown) += dpsidx(l2,0)*psi(l)*W;
196  }
197  }
198  }
199 
200  //Get the local equation number
201  local_eqn = u_local_eqn(l,1);
202  //IF it's not a boundary condition
203  if(local_eqn >= 0)
204  {
205  //Loop over the shape functions
206  for(unsigned l2=0;l2<n_node;l2++)
207  {
208  local_unknown = u_local_eqn(l2,0);
209  //If at a non-zero degree of freedom add in the entry
210  if(local_unknown >= 0)
211  {
212  //This corresponds to the lower left A block
213  jacobian(local_eqn,local_unknown) += dpsidx(l2,0)*psi(l)*W;
214  }
215  local_unknown = u_local_eqn(l2,1);
216  //If at a non-zero degree of freedom add in the entry
217  if(local_unknown >= 0)
218  {
219  //This corresponds to the lower right B block
220  mass_matrix(local_eqn, local_unknown) += psi(l2)*psi(l)*W;
221  //This corresponds to the lower right \mu B block
222  jacobian(local_eqn,local_unknown) +=
223  EigenproblemShift::Mu*psi(l2)*psi(l)*W;
224  }
225  }
226  }
227  }
228  }
229  } //end_of_fill_in_contribution_to_jacobian_and_mass_matrix
230 
231  /// Return FE representation of function value u(s) at local coordinate s
232  inline double interpolated_u(const Vector<double> &s) const
233  {
234  unsigned n_node = nnode();
235 
236  //Local shape functions
237  Shape psi(n_node);
238 
239  //Find values of basis function
240  this->shape(s,psi);
241 
242  //Initialise value of u
243  double interpolated_u = 0.0;
244 
245  //Loop over the local nodes and sum
246  for(unsigned l=0;l<n_node;l++) {interpolated_u+=u(l)*psi[l];}
247 
248  //Return the interpolated value of the eigenfunction
249  return(interpolated_u);
250  }
251 
252  /// Return FE representation of function value w(s) at local coordinate s
253  inline double interpolated_w(const Vector<double> &s) const
254  {
255  unsigned n_node = nnode();
256 
257  //Local shape functions
258  Shape psi(n_node);
259 
260  //Find values of basis function
261  this->shape(s,psi);
262 
263  //Initialise value of u
264  double interpolated_u = 0.0;
265 
266  //Loop over the local nodes and sum
267  for(unsigned l=0;l<n_node;l++) {interpolated_u+=w(l)*psi[l];}
268 
269  //Return the interpolated value of the eigenfunction
270  return(interpolated_u);
271  }
272 
273 
274 protected:
275 
276  /// \short Shape/test functions and derivs w.r.t. to global coords at
277  /// local coord. s; return Jacobian of mapping
278  virtual double dshape_eulerian(const Vector<double> &s,
279  Shape &psi,
280  DShape &dpsidx) const=0;
281 
282  /// \short Shape/test functions and derivs w.r.t. to global coords at
283  /// integration point ipt; return Jacobian of mapping
284  virtual double dshape_eulerian_at_knot(const unsigned &ipt,
285  Shape &psi,
286  DShape &dpsidx) const=0;
287 
288  /// \short Access function that returns the local equation number
289  /// of the unknown in the problem. Default is to assume that it is the
290  /// first (only) value stored at the node.
291  virtual inline int u_local_eqn(const unsigned &n, const unsigned &i)
292  {return nodal_local_eqn(n,i);}
293 
294  private:
295 
296 };
297 
298 
299 
300 //======================================================================
301 /// QComplexHarmonicElement<NNODE_1D> elements are 1D Elements with
302 /// NNODE_1D nodal points that are used to solve the ComplexHarmonic eigenvalue
303 /// Problem described by ComplexHarmonicEquations.
304 //======================================================================
305 template <unsigned NNODE_1D>
306 class QComplexHarmonicElement : public virtual QElement<1,NNODE_1D>,
308 {
309 
310  public:
311 
312  ///\short Constructor: Call constructors for QElement and
313  /// Poisson equations
314  QComplexHarmonicElement() : QElement<1,NNODE_1D>(),
316 
317  /// \short Required # of `values' (pinned or dofs)
318  /// at node n. Here there are two (u and w)
319  inline unsigned required_nvalue(const unsigned &n) const {return 2;}
320 
321  /// \short Output function overloaded from ComplexHarmonicEquations
322  void output(ostream &outfile)
324 
325  /// \short Output function overloaded from ComplexHarmonicEquations
326  void output(ostream &outfile, const unsigned &Nplot)
327  {ComplexHarmonicEquations::output(outfile,Nplot);}
328 
329 
330 protected:
331 
332 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
333  inline double dshape_eulerian(const Vector<double> &s,
334  Shape &psi,
335  DShape &dpsidx) const
336  {return QElement<1,NNODE_1D>::dshape_eulerian(s,psi,dpsidx);}
337 
338 
339  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
340  /// integration point ipt. Return Jacobian.
341  inline double dshape_eulerian_at_knot(const unsigned& ipt,
342  Shape &psi,
343  DShape &dpsidx) const
344  {return QElement<1,NNODE_1D>::dshape_eulerian_at_knot(ipt,psi,dpsidx);}
345 
346 }; //end_of_QComplexHarmonic_class_definition
347 
348 
349 //==start_of_problem_class============================================
350 /// 1D ComplexHarmonic problem in unit interval.
351 //====================================================================
352 template<class ELEMENT,class EIGEN_SOLVER>
353 class ComplexHarmonicProblem : public Problem
354 {
355 public:
356 
357  /// Constructor: Pass number of elements and pointer to source function
358  ComplexHarmonicProblem(const unsigned& n_element);
359 
360  /// Destructor. Clean up the mesh and solver
361  ~ComplexHarmonicProblem(){delete this->mesh_pt();
362  delete this->eigen_solver_pt();}
363 
364  /// Solve the problem
365  void solve(const unsigned &label);
366 
367  /// \short Doc the solution, pass the number of the case considered,
368  /// so that output files can be distinguished.
369  void doc_solution(const unsigned& label);
370 
371 }; // end of problem class
372 
373 
374 
375 //=====start_of_constructor===============================================
376 /// \short Constructor for 1D ComplexHarmonic problem in unit interval.
377 /// Discretise the 1D domain with n_element elements of type ELEMENT.
378 /// Specify function pointer to source function.
379 //========================================================================
380 template<class ELEMENT,class EIGEN_SOLVER>
382  const unsigned& n_element)
383 {
384  //Create the eigen solver
385  this->eigen_solver_pt() = new EIGEN_SOLVER;
386 
387  //Get the positive eigenvalues, shift is zero by default
388  static_cast<EIGEN_SOLVER*>(eigen_solver_pt())
389  ->get_eigenvalues_right_of_shift();
390 
391  //Set domain length
392  double L=1.0;
393 
394  // Build mesh and store pointer in Problem
395  Problem::mesh_pt() = new OneDMesh<ELEMENT>(n_element,L);
396 
397  // Set the boundary conditions for this problem: By default, all nodal
398  // values are free -- we only need to pin the ones that have
399  // Dirichlet conditions.
400 
401  // Pin the single nodal value at the single node on mesh
402  // boundary 0 (= the left domain boundary at x=0)
403  mesh_pt()->boundary_node_pt(0,0)->pin(0);
404 
405  // Pin the single nodal value at the single node on mesh
406  // boundary 1 (= the right domain boundary at x=1)
407  mesh_pt()->boundary_node_pt(1,0)->pin(0);
408 
409  // Setup equation numbering scheme
410  assign_eqn_numbers();
411 
412 } // end of constructor
413 
414 
415 
416 //===start_of_doc=========================================================
417 /// Doc the solution in tecplot format. Label files with label.
418 //========================================================================
419 template<class ELEMENT,class EIGEN_SOLVER>
421  const unsigned& label)
422 {
423 
424  ofstream some_file;
425  char filename[100];
426 
427  // Number of plot points
428  unsigned npts;
429  npts=5;
430 
431  // Output solution with specified number of plot points per element
432  sprintf(filename,"soln%i.dat",label);
433  some_file.open(filename);
434  mesh_pt()->output(some_file,npts);
435  some_file.close();
436 
437 } // end of doc
438 
439 //=======================start_of_solve==============================
440 /// Solve the eigenproblem
441 //===================================================================
442 template<class ELEMENT,class EIGEN_SOLVER>
444 solve(const unsigned& label)
445 {
446  //Set external storage for the eigenvalues
447  Vector<complex<double> > eigenvalues;
448  //Set external storage for the eigenvectors
449  Vector<DoubleVector> eigenvectors;
450  //Desired number eigenvalues enough to include the first two real
451  //eigenvalues and then the complex conjugate pair. If the number is set
452  //to four then ARPACK prefers the degenerate (real) eigenvalues Mu, Mu.
453  unsigned n_eval=7;
454 
455  //Solve the eigenproblem
456  this->solve_eigenproblem(n_eval,eigenvalues,eigenvectors);
457 
458  //We now need to sort the output based on the size of the real part
459  //of the eigenvalues.
460  //This is because the solver does not necessarily sort the eigenvalues
461  Vector<complex<double> > sorted_eigenvalues = eigenvalues;
462  sort(sorted_eigenvalues.begin(),sorted_eigenvalues.end(),
464 
465  //Read out the smallest eigenvalue
466  complex<double> temp_evalue = sorted_eigenvalues[0];
467  unsigned smallest_index=0;
468  //Loop over the unsorted eigenvalues and find the entry that corresponds
469  //to our second smallest eigenvalue.
470  for(unsigned i=0;i<eigenvalues.size();i++)
471  {
472  //Note that equality tests for doubles are bad, but it was just
473  //sorted data, so should be fine
474  if(eigenvalues[i] == temp_evalue) {smallest_index=i; break;}
475  }
476 
477  //Normalise the eigenvector
478  {
479  //Get the dimension of the eigenvector
480  unsigned dim = eigenvectors[smallest_index].nrow();
481  double length=0.0;
482  //Loop over all the entries
483  for(unsigned i=0;i<dim;i++)
484  {
485  //Add the contribution to the length
486  length += std::pow(eigenvectors[smallest_index][i],2.0);
487  }
488  //Now take the magnitude
489  length = sqrt(length);
490  //Fix the sign
491  if(eigenvectors[smallest_index][0] < 0) {length *= -1.0;}
492  //Finally normalise
493  for(unsigned i=0;i<dim;i++)
494  {
495  eigenvectors[smallest_index][i] /= length;
496  }
497  }
498 
499  //Now assign the second eigenvector to the dofs of the problem
500  this->assign_eigenvector_to_dofs(eigenvectors[smallest_index]);
501  //Output solution for this case (label output files with "1")
502  this->doc_solution(label);
503 
504  char filename[100];
505  sprintf(filename,"eigenvalues%i.dat",label);
506 
507  //Open an output file for the sorted eigenvalues
508  ofstream evalues(filename);
509  for(unsigned i=0;i<n_eval;i++)
510  {
511  //Print to screen
512  cout << sorted_eigenvalues[i].real() << " "
513  << sorted_eigenvalues[i].imag() << std::endl;
514  //Send to file
515  evalues << sorted_eigenvalues[i].real() << " "
516  << sorted_eigenvalues[i].imag() << std::endl;
517  }
518 
519  evalues.close();
520 } //end_of_solve
521 
522 
523 ////////////////////////////////////////////////////////////////////////
524 ////////////////////////////////////////////////////////////////////////
525 ////////////////////////////////////////////////////////////////////////
526 
527 
528 //======start_of_main==================================================
529 /// Driver for 1D Poisson problem
530 //=====================================================================
531 int main(int argc, char **argv)
532 {
533 //Want to test Trilinos if we have it, so we must initialise MPI
534 //if we have compiled with it
535 #ifdef OOMPH_HAS_MPI
536  MPI_Helpers::init(argc,argv);
537 #endif
538 
539  // Set up the problem:
540  unsigned n_element=100; //Number of elements
541 
542  clock_t t_start1 = clock();
543  //Solve with ARPACK
544  {
546  problem(n_element);
547 
548  std::cout << "Matrix size " << problem.ndof() << std::endl;
549 
550  problem.solve(1);
551  }
552  clock_t t_end1 = clock();
553 
554  clock_t t_start2 = clock();
555  //Solve with LAPACK_QZ
556  {
558  problem(n_element);
559 
560  problem.solve(2);
561  }
562  clock_t t_end2 = clock();
563 
564 #ifdef OOMPH_HAS_TRILINOS
565  clock_t t_start3 = clock();
566 //Solve with Anasazi
567  {
568  ComplexHarmonicProblem<QComplexHarmonicElement<3>,ANASAZI> problem(n_element);
569  problem.solve(3);
570  }
571  clock_t t_end3 = clock();
572 #endif
573 
574  std::cout << "ARPACK TIME: " << (double)(t_end1 - t_start1)/CLOCKS_PER_SEC
575  << std::endl;
576 
577  std::cout << "LAPACK TIME: " << (double)(t_end2 - t_start2)/CLOCKS_PER_SEC
578  << std::endl;
579 
580 #ifdef OOMPH_HAS_TRILINOS
581  std::cout << "ANASAZI TIME: " << (double)(t_end3 - t_start3)/CLOCKS_PER_SEC
582  << std::endl;
583 #endif
584 
585 #ifdef OOMPH_HAS_MPI
586  MPI_Helpers::finalize();
587 #endif
588 
589 } // end of main
590 
591 
592 
593 
594 
595 
596 
597 
598 
599 
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n. Here there are two (u and w)
ComplexHarmonicProblem(const unsigned &n_element)
Constructor: Pass number of elements and pointer to source function.
double interpolated_u(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
virtual double w(const unsigned &n) const
Second eigenfunction value at local node n.
1D ComplexHarmonic problem in unit interval.
double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
void solve(const unsigned &label)
Solve the problem.
~ComplexHarmonicProblem()
Destructor. Clean up the mesh and solver.
double interpolated_w(const Vector< double > &s) const
Return FE representation of function value w(s) at local coordinate s.
int main(int argc, char **argv)
Driver for 1D Poisson problem.
ComplexHarmonicEquations()
Empty Constructor.
void output(ostream &outfile)
Output function overloaded from ComplexHarmonicEquations.
bool operator()(const complex< T > &x, const complex< T > &y) const
Comparison in terms of magnitude of complex number.
void doc_solution(const unsigned &label)
Doc the solution, pass the number of the case considered, so that output files can be distinguished...
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void output(ostream &outfile)
Output the eigenfunction with default number of plot points.
virtual double u(const unsigned &n) const
Access function: First eigenfunction value at local node n Note that solving the eigenproblem does no...
virtual int u_local_eqn(const unsigned &n, const unsigned &i)
Access function that returns the local equation number of the unknown in the problem. Default is to assume that it is the first (only) value stored at the node.
QComplexHarmonicElement()
Constructor: Call constructors for QElement and Poisson equations.
void output(ostream &outfile, const unsigned &Nplot)
Output function overloaded from ComplexHarmonicEquations.
void output(ostream &outfile, const unsigned &nplot)
Output FE representation of soln: x,y,u or x,y,z,u at Nplot plot points.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Assemble the contributions to the jacobian and mass matrices.