From: Wolfgang Bangerth Date: Mon, 6 Jul 2009 22:18:58 +0000 (+0000) Subject: Document half of the program. X-Git-Tag: v8.0.0~7520 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7cd70477902b5ca20963b9d2fef2d7cb4ae7afab;p=dealii.git Document half of the program. git-svn-id: https://svn.dealii.org/trunk@19039 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-36/step-36.cc b/deal.II/examples/step-36/step-36.cc index a6de265997..472adcd01a 100644 --- a/deal.II/examples/step-36/step-36.cc +++ b/deal.II/examples/step-36/step-36.cc @@ -10,9 +10,15 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ - // @sect3{Include files} We bundle the - // "usual" deal.II include files as we - // did in step-4: + // @sect3{Include files} + + // 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 #include #include @@ -34,72 +40,93 @@ #include #include - // PETSc appears here because SLEPc - // depends on this library: + // PETSc appears here because SLEPc + // depends on this library: #include #include - // and then we need to actually import - // the interfaces for solvers that - // SLEPc provides: + // And then we need to actually import + // the interfaces for solvers that + // SLEPc provides: #include - // and some standard C++: + // We also need some standard C++: #include #include - // and the finally, as in previous - // programs, we import all the deal.II - // class and function names into the - // global namespace: + // Finally, as in previous programs, we + // import all the deal.II class and function + // names into the global namespace: using namespace dealii; - // @sect1{The - // EigenvalueProblem - // class template} + // @sect3{The EigenvalueProblem class template} + // Following is the class declaration for the + // main class template. It looks pretty much + // exactly like what has already been shown + // in step-4: template class EigenvalueProblem { -public: - EigenvalueProblem (const std::string &prm_file); - void run (); + public: + EigenvalueProblem (const std::string &prm_file); + void run (); -private: - void make_grid_and_dofs (); - void assemble_system (); - void solve (); - void output_results () const; + private: + void make_grid_and_dofs (); + void assemble_system (); + void solve (); + void output_results () const; - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; - - // These are data types pertaining to - // the generalized problem. - PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; - std::vector eigenfunctions; - std::vector eigenvalues; - - ParameterHandler parameters; + Triangulation triangulation; + FE_Q fe; + DoFHandler 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 those for the eigenfunctions we + // want to compute, along with the + // corresponding eigenvectors: + PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; + std::vector eigenfunctions; + std::vector 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 $U_i=0$. + ConstraintMatrix constraints; }; - // @sect2{Implementation of the - // EigenvalueProblem - // class} + // @sect3{Implementation of the EigenvalueProblem class} - // @sect3{EigenvalueProblem::EigenvalueProblem} + // @sect4{EigenvalueProblem::EigenvalueProblem} + // First 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 EigenvalueProblem::EigenvalueProblem (const std::string &prm_file) - : - fe (1), - dof_handler (triangulation) + : + fe (1), + dof_handler (triangulation) { - - // Declare some of the needed data - // from a file; you can always change - // this! parameters.declare_entry ("Global mesh refinement steps", "5", Patterns::Integer (0, 20), "The number number of times the 1-cell coarse mesh should " @@ -112,13 +139,64 @@ 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); } - // @sect3{EigenvalueProblem::make_grid_and_dofs} + // @sect4{EigenvalueProblem::make_grid_and_dofs} + + // The next function creates a mesh on the + // domain $[-1,1]^d$, 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 $u|_{\partial\Omega}=0$. + // + // 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: + // @code + // 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); + // @code + // instead of the two reinit() + // calls for the stiffness and mass matrices. + // + // 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. + // + // Absent 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 void EigenvalueProblem::make_grid_and_dofs () { @@ -126,50 +204,50 @@ void EigenvalueProblem::make_grid_and_dofs () triangulation.refine_global (parameters.get_integer ("Global mesh refinement steps")); dof_handler.distribute_dofs (fe); - CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(), - dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, csp); - csp.compress (); - - // What is going on here? - - // This does not work! - // stiffness_matrix.reinit (csp); - // mass_matrix.reinit (csp); + DoFTools::make_zero_boundary_constraints (dof_handler, constraints); + constraints.close (); - // But this does... TODO: Fix it! - stiffness_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), + 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(), + mass_matrix.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); - // with this done we stream-out the - // sparsity pattern - std::ofstream out ("constrained_sparsity_pattern.gpl"); - csp.print_gnuplot (out); - - // The next step is to take care of - // the eigenspectrum. In this case, - // the outputs are eigenfunctions and - // eigenvalues. Set the collective - // eigenfunction block to be as big as - // we wanted! + // The next step is to take care of the + // eigenspectrum. In this case, the outputs + // are eigenfunctions and eigenvalues, so + // we set the size of the list of + // eigenfunctions and eigenvalues to be as + // large as asked for in the input file: eigenfunctions .resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions")); for (unsigned int i=0; i void EigenvalueProblem::assemble_system () { @@ -191,20 +269,9 @@ void EigenvalueProblem::assemble_system () potential.initialize (FunctionParser::default_variable_names (), 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 - // boundary constraints for the - // present grid. - ConstraintMatrix constraints; - constraints.clear (); - DoFTools::make_zero_boundary_constraints (dof_handler, constraints); - constraints.close (); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (), @@ -218,12 +285,9 @@ void EigenvalueProblem::assemble_system () potential.value_list (fe_values.get_quadrature_points(), potential_values); - for (unsigned int q_point=0; q_point::assemble_system () potential_values[q_point] * fe_values.shape_value (i, q_point) * fe_values.shape_value (j, q_point) - ) * fe_values.JxW (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); + ) * fe_values.JxW (q_point); } - // Now we have the local system we - // transfer it into the global objects - // and take care of zero boundary - // constraints, + // 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 @@ -256,55 +320,43 @@ void EigenvalueProblem::assemble_system () mass_matrix); } - // and finally, set the matrices and - // vectors in an assembled state. + // 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(); - - // 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 () { - // We start by assigning 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); - // and assign our solver of - // choice. Here we want to use the - // Krylov-Schur solver, which is - // pretty darn fast and robust: + // 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 + // Lets assign the solver which part + // of the spectrum we want to solve eigensolver.set_which_eigenpairs (EPS_SMALLEST_MAGNITUDE); - // Finally, we actually solve the - // generalized eigenproblem: + // Finally, we actually solve the + // generalized eigenproblem: eigensolver.solve (stiffness_matrix, mass_matrix, eigenvalues, eigenfunctions, eigenfunctions.size()); @@ -318,7 +370,7 @@ void EigenvalueProblem::solve () } - // @sect3{EigenvalueProblem::output_results} + // @sect4{EigenvalueProblem::output_results} template void EigenvalueProblem::output_results () const { @@ -331,7 +383,7 @@ void EigenvalueProblem::output_results () const std::string("eigenfunction_") + Utilities::int_to_string(i)); - // How does this work? + // How does this work? Vector projected_potential (dof_handler.n_dofs()); FunctionParser potential; potential.initialize (FunctionParser::default_variable_names (), @@ -347,7 +399,7 @@ void EigenvalueProblem::output_results () const } - // @sect3{EigenvalueProblem::run} + // @sect4{EigenvalueProblem::run} // This is the function which has the // top-level control over @@ -360,9 +412,9 @@ void EigenvalueProblem::run () make_grid_and_dofs (); - // While we are here, lets count the - // number of active cells and degrees - // of freedom like we always do. + // 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 @@ -386,10 +438,11 @@ int 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 + // 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); { @@ -399,12 +452,14 @@ int main (int argc, char **argv) problem.run (); } - // and then unitialize the SLEPc - // work space when the job is done: + // ...and then unitialize the SLEPc + // work space when the job is done: SlepcFinalize (); } - // or panic if something goes wrong. + // 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 @@ -430,10 +485,9 @@ int main (int argc, char **argv) return 1; } - // or show that we happy and didn't - // crap out on the calculation... + // ...or show that we are happy: std::cout << std::endl - << "Completed" + << "Done." << std::endl; return 0;