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