Table of contents | |
---|---|
This program was contributed by Toby D. Young and Wolfgang Bangerth.
The problem we want to solve in this example is an eigenspectrum problem. Eigenvalue problems appear in a wide context of problems, for example in the computation of electromagnetic standing waves in cavities, vibration modes of drum membranes, or oscillations of lakes and estuaries. One of the most enigmatic applications is probably the computation of stationary or quasi-static wave functions in quantum mechanics. The latter application is what we would like to investigate here, though the general techniques outlined in this program are of course equally applicable to the other applications above.
Eigenspectrum problems have the general form
where the Dirichlet boundary condition on could also be replaced by Neumann or Robin conditions;
is an operator that generally also contains differential operators.
Under suitable conditions, the above equations have a set of solutions ,
, where
can be a finite or infinite set (and in the latter case it may be a discrete or sometimes at least in part a continuous set). In either case, let us note that there is no longer just a single solution, but a set of solutions (the various eigenfunctions and corresponding eigenvalues) that we want to compute. The problem of numerically finding all eigenvalues (eigenfunctions) of such eigenvalue problems is a formidable challenge. In fact, if the set
is infinite, the challenge is of course intractable. Most of the time however we are really only interested in a small subset of these values (functions); and fortunately, the interface to the SLEPc library that we will use for this tutorial program allows us to select which portion of the eigenspectrum and how many solutions we want to solve for.
In this program, the eigenspectrum solvers we use are classes provided by deal.II that wrap around the linear algebra implementation of the SLEPc library; SLEPc itself builds on the PETSc library for linear algebra contents.
The basic equation of stationary quantum mechanics is the Schrödinger equation which models the motion of particles in an external potential . The particle is described by a wave function
that satisfies a relation of the (non-dimensionalized) form
As a consequence, this particle can only exist in a certain number of eigenstates that correspond to the energy eigenvalues admitted as solutions of this equation. The orthodox (Copenhagen) interpretation of quantum mechanics posits that, if a particle has energy
then the probability of finding it at location
is proportional to
where
is the eigenfunction that corresponds to this eigenvalue.
In order to numerically find solutions to this equation, i.e. a set of pairs of eigenvalues/eigenfunctions, we use the usual finite element approach of multiplying the equation from the left with testfunctions, integrating by parts, and searching for solutions in finite dimensional spaces by approximating , where
is a vector of expansion coefficients. We then immediately arrive at the following equation that discretizes the continuous eigenvalue problem:
In matrix and vector notation, this equation then reads:
where is the stiffness matrix arising from the differential operator
, and
is the mass matrix. The solution to the eigenvalue problem is an eigenspectrum
, with associated eigenfunctions
.
It is this form of the eigenvalue problem that involves both matrices and
that we will solve in the current tutorial program. We will want to solve it for the lowermost few eigenvalue/eigenfunction pairs.
The program below is essentially just a slightly modified version of step-4. The things that are different are the following:
The main class (now named EigenvalueProblem
) now no longer has a single solution vector, but a whole set of vectors for the various eigenfunctions we want to compute.
We use PETSc matrices and vectors as in step-17 and step-18 since that is what the SLEPc eigenvalue solvers require.
The function EigenvalueProblem::solve
is entirely different from anything seen so far in the tutorial, as it does not just solve a linear system but actually solves the eigenvalue problem. It is built on the SLEPc library, and more immediately on the deal.II SLEPc wrappers in the class SLEPcWrappers::SolverKrylovSchur.
We use the ParameterHandler class to describe a few input parameters, such as the exact form of the potential , the number of global refinement steps of the mesh, or the number of eigenvalues we want to solve for. We could go much further with this but stop at making only a few of the things that one could select at run time actual input file parameters. In order to see what could be done in this regard, take a look at step-29, step-33, and in particular step-19.
We use the FunctionParser class to make the potential a run-time parameter that can be specified in the input file as a formula.
The rest of the program follows in a pretty straightforward way from step-4.
As mentioned in the introduction, this program is essentially only a slightly revised version of step-4. As a consequence, most of the following include files are as used there, or at least as used already in previous tutorial programs:
#include <base/logstream.h> #include <base/quadrature_lib.h> #include <base/function.h> #include <base/function_parser.h> #include <base/parameter_handler.h> #include <base/utilities.h> #include <grid/tria.h> #include <grid/grid_generator.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <dofs/dof_handler.h> #include <dofs/dof_accessor.h> #include <dofs/dof_tools.h> #include <fe/fe_q.h> #include <fe/fe_values.h> #include <numerics/vectors.h> #include <numerics/matrices.h> #include <numerics/data_out.h> #include <lac/full_matrix.h>
PETSc appears here because SLEPc depends on this library:
#include <lac/petsc_sparse_matrix.h> #include <lac/petsc_vector.h>
And then we need to actually import the interfaces for solvers that SLEPc provides:
#include <lac/slepc_solver.h>
We also need some standard C++:
#include <fstream> #include <iostream>
Finally, as in previous programs, we import all the deal.II class and function names into the global namespace:
using namespace dealii;
EigenvalueProblem
class templateFollowing is the class declaration for the main class template. It looks pretty much exactly like what has already been shown in step-4:
template <int dim> class EigenvalueProblem { public: EigenvalueProblem (const std::string &prm_file); void run (); private: void make_grid_and_dofs (); void assemble_system (); void solve (); void output_results () const; Triangulation<dim> triangulation; FE_Q<dim> fe; DoFHandler<dim> dof_handler;
With these exceptions: For our eigenvalue problem, we need both a stiffness matrix for the left hand side as well as a mass matrix for the right hand side. We also need not just one solution function, but a whole set of these for the eigenfunctions we want to compute, along with the corresponding eigenvalues:
PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; std::vector<PETScWrappers::Vector> eigenfunctions; std::vector<double> eigenvalues;
And then we need an object that will store several run-time parameters that we will specify in an input file :
ParameterHandler parameters;
Finally, we will have an object that contains "constraints" on our degrees of freedom. This could include hanging node constraints if we had adaptively refined meshes (which we don't have in the current program). Here, we will store the constraints for boundary nodes .
ConstraintMatrix constraints; };
EigenvalueProblem
classFirst up, the constructor. The main new part is handling the run-time input parameters. We need to declare their existence first, and then read their values from the input file whose name is specified as an argument to this function:
template <int dim> EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file) : fe (1), dof_handler (triangulation) { parameters.declare_entry ("Global mesh refinement steps", "5", Patterns::Integer (0, 20), "The number of times the 1-cell coarse mesh should " "be refined globally for our computations."); parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "5", Patterns::Integer (0, 100), "The number of eigenvalues/eigenfunctions " "to be computed."); parameters.declare_entry ("Potential", "0", Patterns::Anything(), "A functional description of the potential."); parameters.read_input (prm_file); }
The next function creates a mesh on the domain , refines it as many times as the input file calls for, and then attaches a DoFHandler to it and initializes the matrices and vectors to their correct sizes. We also build the constraints that correspond to the boundary values
.
For the matrices, we use the PETSc wrappers. These have the ability to allocate memory as necessary as non-zero entries are added. This seems inefficient: we could as well first compute the sparsity pattern, initialize the matrices with it, and as we then insert entries we can be sure that we do not need to re-allocate memory and free the one used previously. One way to do that would be to use code like this:
CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(), dof_handler.n_dofs()); DoFTools::make_sparsity_pattern (dof_handler, csp); csp.compress (); stiffness_matrix.reinit (csp); mass_matrix.reinit (csp);
instead of the two reinit()
calls for the stiffness and mass matrices below.
This doesn't quite work, unfortunately. The code above may lead to a few entries in the non-zero pattern to which we only ever write zero entries; most notably, this holds true for off-diagonal entries for those rows and columns that belong to boundary nodes. This shouldn't be a problem, but for whatever reason, PETSc's ILU preconditioner, which we use to solve linear systems in the eigenvalue solver, doesn't like these extra entries and aborts with an error message.
In the absense of any obvious way to avoid this, we simply settle for the second best option, which is have PETSc allocate memory as necessary. That said, since this is not a time critical part, this whole affair is of no further importance.
template <int dim> void EigenvalueProblem<dim>::make_grid_and_dofs () { GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (parameters.get_integer ("Global mesh refinement steps")); dof_handler.distribute_dofs (fe); DoFTools::make_zero_boundary_constraints (dof_handler, constraints); constraints.close (); stiffness_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); mass_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs());
The next step is to take care of the eigenspectrum. In this case, the outputs are eigenvalues and eigenfunctions, so we set the size of the list of eigenfunctions and eigenvalues to be as large as we asked for in the input file :
eigenfunctions .resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions")); for (unsigned int i=0; i<eigenfunctions.size (); ++i) eigenfunctions[i].reinit (dof_handler.n_dofs ()); eigenvalues.resize (eigenfunctions.size ()); }
Here, we assemble the global stiffness and mass matrices from local contributions and
respectively. This function should be immediately familiar if you've seen previous tutorial programs. The only thing new would be setting up an object that described the potential
using the expression that we got from the input file. We then need to evaluate this object at the quadrature points on each cell. If you've seen how to evaluate function objects (see, for example the coefficient in step-5), the code here will also look rather familiar.
template <int dim> void EigenvalueProblem<dim>::assemble_system () { QGauss<dim> quadrature_formula(2); FEValues<dim> fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell); std::vector<unsigned int> local_dof_indices (dofs_per_cell); FunctionParser<dim> potential; potential.initialize (FunctionParser<dim>::default_variable_names (), parameters.get ("Potential"), typename FunctionParser<dim>::ConstMap()); std::vector<double> potential_values (n_q_points); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active (), endc = dof_handler.end (); for (; cell!=endc; ++cell) { fe_values.reinit (cell); cell_stiffness_matrix = 0; cell_mass_matrix = 0; potential.value_list (fe_values.get_quadrature_points(), potential_values); for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) { cell_stiffness_matrix (i, j) += (fe_values.shape_grad (i, q_point) * fe_values.shape_grad (j, q_point) + potential_values[q_point] * fe_values.shape_value (i, q_point) * fe_values.shape_value (j, q_point) ) * fe_values.JxW (q_point); cell_mass_matrix (i, j) += (fe_values.shape_value (i, q_point) * fe_values.shape_value (j, q_point) ) * fe_values.JxW (q_point); }
Now that we have the local matrix contributions, we transfer them into the global objects and take care of zero boundary constraints:
cell->get_dof_indices (local_dof_indices); constraints .distribute_local_to_global (cell_stiffness_matrix, local_dof_indices, stiffness_matrix); constraints .distribute_local_to_global (cell_mass_matrix, local_dof_indices, mass_matrix); }
At the end of the function, we tell PETSc that the matrices have now been fully assembled and that the sparse matrix representation can now be compressed as no more entries will be added:
stiffness_matrix.compress (); mass_matrix.compress (); }
This is the key new functionality of the program. Now that the system is set up, here is a good time to actually solve the problem: As with other examples this is done using a "solve" routine. Essentially, it works as in other programs: you set up a SolverControl object that describes the accuracy to which we want to solve the linear systems, and then we select the kind of solver we want. Here we choose the Krylov-Schur solver of SLEPc, a pretty fast and robust choice for this kind of problem:
template <int dim> void EigenvalueProblem<dim>::solve () {
We start here, as we normally do, by assigning convergence control we want:
SolverControl solver_control (dof_handler.n_dofs(), 1e-9); SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
Before we actually solve for the eigenfunctions and -values, we have to also select which set of eigenvalues to solve for. Lets select those eigenvalues and corresponding eigenfunctions with the smallest real part (in fact, the problem we solve here is symmetric and so the eigenvalues are purely real). After that, we can actually let SLEPc do its work:
eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); eigensolver.solve (stiffness_matrix, mass_matrix, eigenvalues, eigenfunctions, eigenfunctions.size());
The output of the call above is a set of vectors and values. In eigenvalue problems, the eigenfunctions are only determined up to a constant that can be fixed pretty arbitrarily. Knowing nothing about the origin of the eigenvalue problem, SLEPc has no other choice than to normalize the eigenvectors to one in the (vector) norm. Unfortunately this norm has little to do with any norm we may be interested from a eigenfunction perspective: the
norm, or maybe the
norm.
Let us choose the latter and rescale eigenfunctions so that they have instead of
(where
is the
th eigenfunction and
the corresponding vector of nodal values). For the
elements chosen here, we know that the maximum of the function
is attained at one of the nodes, so
, making the normalization in the
norm trivial. Note that this doesn't work as easily if we had chosen
elements with
: there, the maximum of a function does not necessarily have to be attained at a node, and so
(although the equality is usually nearly true).
for (unsigned int i=0; i<eigenfunctions.size(); ++i) eigenfunctions[i] /= eigenfunctions[i].linfty_norm (); }
This is the last significant function of this program. It uses the DataOut class to generate graphical output from the eigenfunctions for later visualization. It works as in many of the other tutorial programs.
The whole collection of functions is then output as a single VTK file.
template <int dim> void EigenvalueProblem<dim>::output_results () const { DataOut<dim> data_out; data_out.attach_dof_handler (dof_handler); for (unsigned int i=0; i<eigenfunctions.size(); ++i) data_out.add_data_vector (eigenfunctions[i], std::string("eigenfunction_") + Utilities::int_to_string(i));
The only thing worth discussing may be that because the potential is specified as a function expression in the input file, it would be nice to also have it as a graphical representation along with the eigenfunctions. The process to achieve this is relatively straightforward: we build an object that represents and then we interpolate this continuous function onto the finite element space. The result we also attach to the DataOut object for visualization.
Vector<double> projected_potential (dof_handler.n_dofs()); { FunctionParser<dim> potential; potential.initialize (FunctionParser<dim>::default_variable_names (), parameters.get ("Potential"), typename FunctionParser<dim>::ConstMap()); VectorTools::interpolate (dof_handler, potential, projected_potential); } data_out.add_data_vector (projected_potential, "interpolated_potential"); data_out.build_patches (); std::ofstream output ("eigenvectors.vtk"); data_out.write_vtk (output); }
This is the function which has the top-level control over everything. It is almost exactly the same as in step-4:
template <int dim> void EigenvalueProblem<dim>::run () { make_grid_and_dofs (); std::cout << " Number of active cells: " << triangulation.n_active_cells () << std::endl << " Number of degrees of freedom: " << dof_handler.n_dofs () << std::endl << std::endl; assemble_system (); solve (); output_results (); for (unsigned int i=0; i<eigenvalues.size(); ++i) std::cout << " Eigenvalue " << i << " : " << eigenvalues[i] << std::endl; }
main
functionint main (int argc, char **argv) { try {
Here is another difference from other steps: We initialize the SLEPc work space which inherently initializes the PETSc work space, run the whole program, ...
SlepcInitialize (&argc, &argv, 0, 0);
{
deallog.depth_console (0);
EigenvalueProblem<2> problem ("step-36.prm");
problem.run ();
}
...and then unitialize the SLEPc work space when the job is done:
SlepcFinalize (); }
All the while, we are watching out if any exceptions should have been generated. If that is so, we panic...
catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; }
If no exceptions are thrown, then we tell the program to stop monkeying around and exit nicely:
std::cout << std::endl << "Job done." << std::endl; return 0; }
The problem's input is parameterized by an input file step-36.prm
which could, for example, contain the following text:
set Global mesh refinement steps = 5 set Number of eigenvalues/eigenfunctions = 5 set Potential = 0
Here, the potential is zero inside the domain, and we know that the eigenvalues are given by where
. Eigenfunctions are sines and cosines with
and
periods in
and
directions. This matches the output our program generates:
examples/@ref step_36 "step-36"> make run ============================ Running @ref step_36 "step-36" Number of active cells: 1024 Number of degrees of freedom: 1089 Eigenvalue 0 : 4.93877 Eigenvalue 1 : 12.3707 Eigenvalue 2 : 12.3707 Eigenvalue 3 : 19.8027 Eigenvalue 4 : 24.837 Job done.
These eigenvalues are exactly the ones that correspond to pairs ,
and
,
, and
. A visualization of the corresponding eigenfunctions would look like this:
![]() | ![]() |
![]() | ![]() |
![]() |
It is always worth playing a few games in the playground! So here goes with a few suggestions:
The potential used above (called the infinite well because it is a flat potential surrounded by infinitely high walls) is interesting because it allows for analytically known solutions. Apart from that, it is rather boring, however. That said, it is trivial to play around with the potential by just setting it to something different in the input file. For example, let us assume that we wanted to work with the following potential in 2d:
In other words, the potential is -100 in two sectors of a circle of radius 0.75, -5 in the other two sectors, and zero outside the circle. We can achieve this by using the following in the input file :
set Potential = if (x^2 + y^2 < 0.75^2, if (x*y > 0, -100, -5), 0)
If in addition we also increase the mesh refinement by one level, we get the following results:
examples/@ref step_36 "step-36"> make run ============================ Running @ref step_36 "step-36" Number of active cells: 4096 Number of degrees of freedom: 4225 Eigenvalue 0 : -74.2562 Eigenvalue 1 : -72.7322 Eigenvalue 2 : -42.7406 Eigenvalue 3 : -42.2232 Eigenvalue 4 : -37.0744
The output file also contains an interpolated version of the potential, which looks like this (note that as expected the lowest few eigenmodes have probability densities that are significant only where the potential is the lowest, i.e. in the top right and bottom left sector of inner circle of the potential):
The first five eigenfunctions are now like this:
![]() | ![]() |
![]() | ![]() |
![]() |
In our derivation of the problem we have assumed that the particle is confined to a domain and that at the boundary of this domain its probability
of being is zero. This is equivalent to solving the eigenvalue problem on all of
and assuming that the energy potential is finite only inside a region
and infinite outside. It is relatively easy to show that
at all locations
where
. So the question is what happens if our potential is not of this form, i.e. there is no bounded domain outside of which the potential is infinite? In that case, it may be worth to just consider a very large domain at the boundary of which
is at least very large, if not infinite. Play around with a few cases like this and explore how the spectrum and eigenfunctions change as we make the computational region larger and larger.
What happens if we investigate the simple harmonic oscillator problem ? This potential is exactly of the form discussed in the previous paragraph and has hyper spherical symmetry. One may want to use a large spherical domain with a large outer radius, to approximate the whole-space problem (say, by invoking GridGenerator::hyper_ball).
The plots above show the wave function , but the physical quantity of interest is actually the probability density
for the particle to be at location
. Some visualization programs can compute derived quantities from the data in an input file, but we can also do so right away when creating the output file. The facility to do that is the DataPostprocessor class that can be used in conjunction with the DataOut class. Examples of how this can be done can be found in step-29 and step-33.
What happens if the particle in the box has internal degrees of freedom? For example, if the particle were a spin- particle? In that case, we may want to start solving a vector-valued problem instead.
Our implementation of the deal.II library here uses the PETScWrappers and SLEPcWrappers and is suitable for running on serial machine architecture. However, for larger grids and with a larger number of degrees-of-freedom, we may want to run our application on parallel architectures. A parallel implementation of the above code can be particuarily useful here since the generalized eigenspectrum problem is somewhat more expensive to solve than the standard problems considered most of the earlier tutorials. Fortunately, modifying the above program to be MPI complient is a relatively straightforward procedure. A sketch of how this can be done can be found in step-17.
(If you are looking at a locally installed deal.II version, then the program can be found at /home /hazelsct /repositories /deal.ii /examples /step-36 /step-36.cc . Otherwise, this is only the path on some remote server.)
/ * Subversion Id: @ref step_36 "step-36".cc 19132 2009-07-30 12:03:40Z young * / / * Author: Toby D. Young, Polish Academy of Sciences, * / / * Wolfgang Bangerth, Texas A&M University * / / * Subversion Id: @ref step_36 "step-36".cc 19132 2009-07-30 12:03:40Z young * / / * * / / * Copyright (C) 2009 by the deal.II authors * / / * * / / * This file is subject to QPL and may not be distributed * / / * without copyright and license information. Please refer * / / * to the file deal.II/doc/license.html for the text and * / / * further information on this license. * /
#include <base/logstream.h> #include <base/quadrature_lib.h> #include <base/function.h> #include <base/function_parser.h> #include <base/parameter_handler.h> #include <base/utilities.h> #include <grid/tria.h> #include <grid/grid_generator.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <dofs/dof_handler.h> #include <dofs/dof_accessor.h> #include <dofs/dof_tools.h> #include <fe/fe_q.h> #include <fe/fe_values.h> #include <numerics/vectors.h> #include <numerics/matrices.h> #include <numerics/data_out.h> #include <lac/full_matrix.h> #include <lac/petsc_sparse_matrix.h> #include <lac/petsc_vector.h> #include <lac/slepc_solver.h> #include <fstream> #include <iostream> using namespace dealii;
template <int dim> class EigenvalueProblem { public: EigenvalueProblem (const std::string &prm_file); void run (); private: void make_grid_and_dofs (); void assemble_system (); void solve (); void output_results () const; Triangulation<dim> triangulation; FE_Q<dim> fe; DoFHandler<dim> dof_handler; PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; std::vector<PETScWrappers::Vector> eigenfunctions; std::vector<double> eigenvalues; ParameterHandler parameters; ConstraintMatrix constraints; };
template <int dim> EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file) : fe (1), dof_handler (triangulation) { parameters.declare_entry ("Global mesh refinement steps", "5", Patterns::Integer (0, 20), "The number of times the 1-cell coarse mesh should " "be refined globally for our computations."); parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "5", Patterns::Integer (0, 100), "The number of eigenvalues/eigenfunctions " "to be computed."); parameters.declare_entry ("Potential", "0", Patterns::Anything(), "A functional description of the potential."); parameters.read_input (prm_file); }
template <int dim> void EigenvalueProblem<dim>::make_grid_and_dofs () { GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (parameters.get_integer ("Global mesh refinement steps")); dof_handler.distribute_dofs (fe); DoFTools::make_zero_boundary_constraints (dof_handler, constraints); constraints.close (); stiffness_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); mass_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); eigenfunctions .resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions")); for (unsigned int i=0; i<eigenfunctions.size (); ++i) eigenfunctions[i].reinit (dof_handler.n_dofs ()); eigenvalues.resize (eigenfunctions.size ()); }
template <int dim> void EigenvalueProblem<dim>::assemble_system () { QGauss<dim> quadrature_formula(2); FEValues<dim> fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell); std::vector<unsigned int> local_dof_indices (dofs_per_cell); FunctionParser<dim> potential; potential.initialize (FunctionParser<dim>::default_variable_names (), parameters.get ("Potential"), typename FunctionParser<dim>::ConstMap()); std::vector<double> potential_values (n_q_points); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active (), endc = dof_handler.end (); for (; cell!=endc; ++cell) { fe_values.reinit (cell); cell_stiffness_matrix = 0; cell_mass_matrix = 0; potential.value_list (fe_values.get_quadrature_points(), potential_values); for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) { cell_stiffness_matrix (i, j) += (fe_values.shape_grad (i, q_point) * fe_values.shape_grad (j, q_point) + potential_values[q_point] * fe_values.shape_value (i, q_point) * fe_values.shape_value (j, q_point) ) * fe_values.JxW (q_point); cell_mass_matrix (i, j) += (fe_values.shape_value (i, q_point) * fe_values.shape_value (j, q_point) ) * fe_values.JxW (q_point); } cell->get_dof_indices (local_dof_indices); constraints .distribute_local_to_global (cell_stiffness_matrix, local_dof_indices, stiffness_matrix); constraints .distribute_local_to_global (cell_mass_matrix, local_dof_indices, mass_matrix); } stiffness_matrix.compress (); mass_matrix.compress (); }
template <int dim> void EigenvalueProblem<dim>::solve () { SolverControl solver_control (dof_handler.n_dofs(), 1e-9); SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control); eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); eigensolver.solve (stiffness_matrix, mass_matrix, eigenvalues, eigenfunctions, eigenfunctions.size()); for (unsigned int i=0; i<eigenfunctions.size(); ++i) eigenfunctions[i] /= eigenfunctions[i].linfty_norm (); }
template <int dim> void EigenvalueProblem<dim>::output_results () const { DataOut<dim> data_out; data_out.attach_dof_handler (dof_handler); for (unsigned int i=0; i<eigenfunctions.size(); ++i) data_out.add_data_vector (eigenfunctions[i], std::string("eigenfunction_") + Utilities::int_to_string(i)); Vector<double> projected_potential (dof_handler.n_dofs()); { FunctionParser<dim> potential; potential.initialize (FunctionParser<dim>::default_variable_names (), parameters.get ("Potential"), typename FunctionParser<dim>::ConstMap()); VectorTools::interpolate (dof_handler, potential, projected_potential); } data_out.add_data_vector (projected_potential, "interpolated_potential"); data_out.build_patches (); std::ofstream output ("eigenvectors.vtk"); data_out.write_vtk (output); }
template <int dim> void EigenvalueProblem<dim>::run () { make_grid_and_dofs (); std::cout << " Number of active cells: " << triangulation.n_active_cells () << std::endl << " Number of degrees of freedom: " << dof_handler.n_dofs () << std::endl << std::endl; assemble_system (); solve (); output_results (); for (unsigned int i=0; i<eigenvalues.size(); ++i) std::cout << " Eigenvalue " << i << " : " << eigenvalues[i] << std::endl; }
int main (int argc, char **argv) { try { SlepcInitialize (&argc, &argv, 0, 0); { deallog.depth_console (0); EigenvalueProblem<2> problem ("step-36.prm"); problem.run (); } SlepcFinalize (); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } std::cout << std::endl << "Job done." << std::endl; return 0; }