From 3816a45609f4a8ef208ffcaadc53a46230009be6 Mon Sep 17 00:00:00 2001 From: "Toby D. Young" Date: Mon, 6 Jul 2009 07:14:29 +0000 Subject: [PATCH] Fixed a bug in the matrix.init () git-svn-id: https://svn.dealii.org/trunk@19014 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-36/doc/intro.dox | 17 +- deal.II/examples/step-36/step-36.cc | 226 ++++++++++++++++--------- deal.II/examples/step-36/step-36.prm | 2 +- 3 files changed, 165 insertions(+), 80 deletions(-) diff --git a/deal.II/examples/step-36/doc/intro.dox b/deal.II/examples/step-36/doc/intro.dox index 5aa971f5c6..e6abfb85b0 100644 --- a/deal.II/examples/step-36/doc/intro.dox +++ b/deal.II/examples/step-36/doc/intro.dox @@ -7,8 +7,21 @@ Bangerth.

Introduction

+Everything here is built on top of classes provided by deal.II that +wrap around the linear algebra implementation of the SLEPc +library; which wraps itself around the PETSc 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$. diff --git a/deal.II/examples/step-36/step-36.cc b/deal.II/examples/step-36/step-36.cc index 008ee2bb21..ba24fd7c45 100644 --- a/deal.II/examples/step-36/step-36.cc +++ b/deal.II/examples/step-36/step-36.cc @@ -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 #include #include @@ -19,10 +20,6 @@ #include #include #include -#include -#include -#include -#include #include #include #include @@ -37,21 +34,27 @@ #include #include - // We need the interfaces for solvers - // that SLEPc provides: + // PETSc appears here because SLEPc + // depends on him: +#include +#include + + // and then we need to actually import + // the interfaces for solvers that + // SLEPc provides: #include - // and then some standard C++: + // and some standard C++: #include #include - // 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 // EigenvalueProblem // class template} @@ -68,9 +71,9 @@ private: void solve (); void output_results () const; - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; + Triangulation triangulation; + FE_Q fe; + DoFHandler 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 // EigenvalueProblem // class} - // @sect4{EigenvalueProblem::EigenvalueProblem} + // @sect3{EigenvalueProblem::EigenvalueProblem} template EigenvalueProblem::EigenvalueProblem (const std::string &prm_file) @@ -93,6 +96,10 @@ EigenvalueProblem::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::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::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 void EigenvalueProblem::assemble_system () { QGauss quadrature_formula(2); FEValues 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::assemble_system () "x,y,z"), parameters.get ("Potential"), typename FunctionParser::ConstMap()); + + // Here we call the actual + // quadrature-point values of the + // potential defined in the prm file + // which we read in earlier std::vector potential_values (n_q_points); // Initialize objects denoting zero @@ -159,13 +185,12 @@ void EigenvalueProblem::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::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::assemble_system () potential.value_list (fe_values.get_quadrature_points(), potential_values); - for (unsigned int q_point=0; q_point::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 void EigenvalueProblem::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 void EigenvalueProblem::output_results () const { @@ -283,32 +323,54 @@ void EigenvalueProblem::output_results () const std::ofstream output ("eigenvectors.vtk"); data_out.write_vtk (output); + + for (unsigned int i=0; i void EigenvalueProblem::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 main function} + // @sect9{The main 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; } diff --git a/deal.II/examples/step-36/step-36.prm b/deal.II/examples/step-36/step-36.prm index 11209580aa..42efdad73d 100644 --- a/deal.II/examples/step-36/step-36.prm +++ b/deal.II/examples/step-36/step-36.prm @@ -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) -- 2.39.5