]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fixed a bug in the matrix.init ()
authorToby D. Young <tyoung@ippt.pan.pl>
Mon, 6 Jul 2009 07:14:29 +0000 (07:14 +0000)
committerToby D. Young <tyoung@ippt.pan.pl>
Mon, 6 Jul 2009 07:14:29 +0000 (07:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@19014 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-36/doc/intro.dox
deal.II/examples/step-36/step-36.cc
deal.II/examples/step-36/step-36.prm

index 5aa971f5c6d6a768502b36670182964d49575afa..e6abfb85b0053cfdf7c696c6c443c78d59508d9d 100644 (file)
@@ -7,8 +7,21 @@ Bangerth.  </i>
 <a name="Intro"></a>
 <h1>Introduction</h1>
 
+Everything here is built on top of classes provided by deal.II that
+wrap around the linear algebra implementation of the <a
+href="http://www.grycap.upv.es/slepc/" target="_top">SLEPc</a>
+library; which wraps itself around the <a
+href="http://www.mcs.anl.gov/petsc/" target="_top">PETSc</a> library.
+
+Start with the wave equation of Schr\"odinger:
 @f[
-       ( \frac{1}{2} \partial_x \partial_y + f(x,y) ) \Phi
+       ( \frac{1}{2} \partial_x \partial_y + V(x,y) ) \Phi
        = 
        \varepsilon \Psi \quad,
-@f]
\ No newline at end of file
+@f]
+Then transform using the finite element {\it ansatz}
+@f[
+       \phi = \tilde{\phi} \quad,
+@f]
+where $\tilde{\phi_i}$ is known as the reduced wavefunction of the
+state $i$.
index 008ee2bb216e6d9f6a5ea18b2069534662e8f0b7..ba24fd7c452abaad854ab2b1deb1e7e13a815b40 100644 (file)
@@ -10,8 +10,9 @@
 /*    to the file deal.II/doc/license.html for the  text  and     */
 /*    further information on this license.                        */
 
-                                // @sect3{Include files}
-
+                                // @sect3{Include files} We bundle the
+                                // "usual" deal.II include files as we
+                                // did in step-4:
 #include <base/logstream.h>
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/parameter_handler.h>
 #include <base/utilities.h>
 #include <lac/full_matrix.h>
-#include <lac/petsc_sparse_matrix.h>
-#include <lac/petsc_vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
 #include <lac/compressed_simple_sparsity_pattern.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 
-                                // We need the interfaces for solvers
-                                // that SLEPc provides:
+                                // PETSc appears here because SLEPc
+                                // depends on him:
+#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>
 
-                                // and then some standard C++:
+                                // and some standard C++:
 #include <fstream>
 #include <iostream>
 
-                               // and the final step, as in previous
-                               // programs, is to import all the
-                               // deal.II class and function names
-                               // into the global namespace:
+                               // and the finally, as in previous
+                               // programs, we import all the deal.II
+                               // class and function names into the
+                               // global namespace:
 using namespace dealii;
 
-                                // @sect3{The
+                                // @sect1{The
                                 // <code>EigenvalueProblem</code>
                                 // class template}
 
@@ -68,9 +71,9 @@ private:
   void solve ();
   void output_results () const;
   
-  Triangulation<dim>   triangulation;
-  FE_Q<dim>            fe;
-  DoFHandler<dim>      dof_handler;
+  Triangulation<dim> triangulation;
+  FE_Q<dim>          fe;
+  DoFHandler<dim>    dof_handler;
   
                                 // These are data types pertaining to
                                 // the generalized problem.
@@ -81,11 +84,11 @@ private:
   ParameterHandler parameters;
 };
 
-                                // @sect3{Implementation of the
+                                // @sect2{Implementation of the
                                 // <code>EigenvalueProblem</code>
                                 // class}
 
-                                // @sect4{EigenvalueProblem::EigenvalueProblem}
+                                // @sect3{EigenvalueProblem::EigenvalueProblem}
 
 template <int dim>
 EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
@@ -93,6 +96,10 @@ EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
   fe (1),
   dof_handler (triangulation)
 {
+
+                                // Declare some of the needed data
+                                // from a file; you can always change
+                                // this!
   parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "10",
                            Patterns::Integer (0, 100),
                            "The number of eigenvalues/eigenfunctions "
@@ -101,6 +108,8 @@ EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
                            Patterns::Anything(),
                            "A functional description of the potential.");
   
+                                // Entries are declared, so now we
+                                // read them in...
   parameters.read_input (prm_file);
 }
 
@@ -113,29 +122,41 @@ void EigenvalueProblem<dim>::make_grid_and_dofs ()
   triangulation.refine_global (5);
   dof_handler.distribute_dofs (fe);
 
-  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
-                                      dof_handler.n_dofs());
-  DoFTools::make_sparsity_pattern (dof_handler, csp);
-  stiffness_matrix.reinit (csp);
-  mass_matrix.reinit (csp);
+  //   CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+  //                                  dof_handler.n_dofs());
+  //   DoFTools::make_sparsity_pattern (dof_handler, csp);
+  //   stiffness_matrix.reinit (csp);
+  //   mass_matrix.reinit (csp);
 
+  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());
+  
+                                // Set the collective eigenfunction
+                                // block to be as big as we wanted!
   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());
+  for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+    eigenfunctions[i].reinit (dof_handler.n_dofs ());
+
+                                // and do the same for the eigenvalue
+                                // output, which had better be the
+                                // same size as the eigenfunction
+                                // block
   eigenvalues
-    .resize (eigenfunctions.size());
+    .resize (eigenfunctions.size ());
 }
 
 
-                                // @sect4{EigenvalueProblem::assemble_system}
+                                // @sect5{EigenvalueProblem::assemble_system}
 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_values | update_gradients |
                            update_quadrature_points | update_JxW_values);
   
   const unsigned int dofs_per_cell = fe.dofs_per_cell;
@@ -152,6 +173,11 @@ void EigenvalueProblem<dim>::assemble_system ()
                         "x,y,z"),
                        parameters.get ("Potential"),
                        typename FunctionParser<dim>::ConstMap());
+
+                                // Here we call the actual
+                                // quadrature-point values of the
+                                // potential defined in the prm file
+                                // which we read in earlier
   std::vector<double> potential_values (n_q_points);
   
                                 // Initialize objects denoting zero
@@ -159,13 +185,12 @@ void EigenvalueProblem<dim>::assemble_system ()
                                 // present grid.
   ConstraintMatrix constraints;
   constraints.clear ();
-  DoFTools::make_zero_boundary_constraints (dof_handler,
-                                           constraints);
+  DoFTools::make_zero_boundary_constraints (dof_handler, constraints);
   constraints.close ();
   
   typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
+    cell = dof_handler.begin_active (),
+    endc = dof_handler.end ();
   for (; cell!=endc; ++cell)
     {
       fe_values.reinit (cell);
@@ -175,24 +200,28 @@ void EigenvalueProblem<dim>::assemble_system ()
       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)
+      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)
-               += ((0.5                               *
-                    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_stiffness_matrix (i, j)
+               += (fe_values.shape_grad (i, q_point) *
+                   0.5                               *
+                   fe_values.shape_grad (j, q_point) 
+                   + 
+                   fe_values.shape_value (i, q_point) *
+                   //               potential_values[q_point] *
+                   0 * // infinite potential well
+                   fe_values.shape_value (j, q_point)
+                   ) * fe_values.JxW (q_point);
              
-             cell_mass_matrix(i,j)
+             cell_mass_matrix (i, j)
                += (fe_values.shape_value (i, q_point) *
-                   fe_values.shape_value (j, q_point) *
-                   fe_values.JxW (q_point));
+                   fe_values.shape_value (j, q_point) 
+                   ) * fe_values.JxW (q_point);
            }
 
                                // Now we have the local system we
@@ -216,46 +245,57 @@ void EigenvalueProblem<dim>::assemble_system ()
   stiffness_matrix.compress();
   mass_matrix.compress();
 
-                                  // make sure that the diagonal entries of
-                                  // constrained degrees of freedom are
-                                  // non-zero to ensure that the matrix is
-                                  // actually invertible
-  for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-    if (constraints.is_constrained(i))
-      {
-       stiffness_matrix.set (i, i, 1);
-       mass_matrix.set (i, i, 1);
-      }
+                               // make sure that the diagonal entries
+                               // of constrained degrees of freedom
+                               // are non-zero to ensure that the
+                               // matrix is actually invertible
+  //   for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+  //     if (constraints.is_constrained(i))
+  //       {
+  //   stiffness_matrix.set (i, i, 1);
+  //   mass_matrix.set (i, i, 1);
+  //       }
   
-  stiffness_matrix.compress();
-  mass_matrix.compress();  
+                                // finally set the matrices in an
+                                // assembled state so SLEPc likes:
+  //   stiffness_matrix.compress ();
+  //   mass_matrix.compress ();  
 }
 
 
-                                // @sect4{EigenvalueProblem::solve}
+                                // @sect6{EigenvalueProblem::solve}
+                                // 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
 template <int dim>
 void EigenvalueProblem<dim>::solve () 
 {
-                                // assign the accuracy to which we
-                                // would like to solve the system,
+                                // We start by assigning the accuracy
+                                // to which we would like to solve the
+                                // system,
   SolverControl solver_control (1000, 1e-6);
-  SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
 
-                                // choose which part of the spectrum to
-                                // solve for
+                                // and assign our solver of
+                                // choice. Here we want to use the
+                                // Krylov-Schur solver, which is
+                                // pretty darn fast and robust:
+  SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
 
+                                // Lets assign the solver which part
+                                // of the spectrum we want to solve
   eigensolver.set_which_eigenpairs (EPS_SMALLEST_MAGNITUDE);
 
-                                // then actually solve the system,
-
+                                // Finally, we actually solve the
+                                // generalized eigenproblem:
   eigensolver.solve (stiffness_matrix, mass_matrix, 
                     eigenvalues, eigenfunctions, 
                     eigenfunctions.size());
 }
 
 
-                                // @sect4{EigenvalueProblem::output_results}
-
+                                // @sect7{EigenvalueProblem::output_results}
 template <int dim>
 void EigenvalueProblem<dim>::output_results () const
 {
@@ -283,32 +323,54 @@ void EigenvalueProblem<dim>::output_results () const
 
   std::ofstream output ("eigenvectors.vtk");
   data_out.write_vtk (output);
+
+  for (unsigned int i=0; i<eigenvalues.size(); ++i)
+    std::cout << std::endl 
+             << "      eigenvalue " << i 
+             << " : " << eigenvalues[i];
+
 }
 
 
-                                 // @sect4{EigenvalueProblem::run}
+                                 // @sect8{EigenvalueProblem::run}
 
                                  // This is the function which has the
-                                // top-level control over everything. It is
-                                // the same as for the previous example.
+                                // top-level control over
+                                // everything. It is very similar as
+                                // for the previous examples.
 template <int dim>
 void EigenvalueProblem<dim>::run () 
 {
   std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
   
-  make_grid_and_dofs();
+  make_grid_and_dofs ();
+
+                                 // While we are here, lets count the
+                                 // number of active cells and degrees
+                                 // of freedom like we always do.
+  std::cout << "   Number of active cells:       "
+           << triangulation.n_active_cells()
+           << std::endl
+           << "   Number of degrees of freedom: "
+           << dof_handler.n_dofs();
+  
   assemble_system ();
   solve ();
   output_results ();
 }
 
 
-                                 // @sect3{The <code>main</code> function}
+                                 // @sect9{The <code>main</code> function}
 int main (int argc, char **argv) 
 {
   try
     {
-      SlepcInitialize(&argc,&argv,0,0);
+
+                                 // Here is another difference from
+                                 // other steps: We initialize the
+                                 // SLEPc work space which inherently
+                                 // initializes the PETSc work space
+      SlepcInitialize (&argc,&argv,0,0);
 
       {
        deallog.depth_console (0);
@@ -317,8 +379,12 @@ int main (int argc, char **argv)
        problem.run ();
       }
 
-      SlepcFinalize();
+                                 // and then unitialize the SLEPc
+                                 // work space when the job is done:
+      SlepcFinalize ();
     }
+
+                                 // or panic if something goes wrong.
   catch (std::exception &exc)
     {
       std::cerr << std::endl << std::endl
@@ -344,5 +410,11 @@ int main (int argc, char **argv)
       return 1;
     }
   
+                                 // or show that we happy and didn't
+                                 // crap out on the calculation...
+  std::cout << std::endl 
+           << "Completed" 
+           << std::endl;
+
   return 0;
 }
index 11209580aa416840acb0204ded52843b8a4cc922..42efdad73dcabc3a579be38da888ab1e0e947c5e 100644 (file)
@@ -1,2 +1,2 @@
-set Number of eigenvalues/eigenfunctions = 5
+set Number of eigenvalues/eigenfunctions = 6
 set Potential = if (x^2 + y^2 < 0.75^2, if (x*y > 0, -10, -5), 0)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.