From: bangerth Date: Fri, 13 Apr 2012 09:30:41 +0000 (+0000) Subject: Implement the framework for a line search but only return 0.1. Getting a proper line... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4ef3ba64fd6d2a178e1e820d69aa191d12f9155e;p=dealii-svn.git Implement the framework for a line search but only return 0.1. Getting a proper line search to converge turns out to be difficult and will be left as an exercise. git-svn-id: https://svn.dealii.org/trunk@25410 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index 10318f62f4..3a111251fb 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -88,6 +88,9 @@ class Step15 void assemble_system (); void solve (); void refine_grid (); + void set_boundary_values (); + double compute_residual (const double alpha) const; + double determine_step_length() const; Triangulation triangulation; @@ -105,7 +108,6 @@ class Step15 - double res; unsigned int refinement; // As described in the Introduction, the first Newton iteration @@ -231,7 +233,7 @@ void Step15::assemble_system () system_rhs = 0; FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | + update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; @@ -285,8 +287,7 @@ void Step15::assemble_system () * fe_values.JxW(q_point)); } - cell_rhs(i) -=0.1 * - (fe_values.shape_grad(i, q_point) * coeff + cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff * gradients[q_point] * fe_values.JxW(q_point)); } } @@ -306,30 +307,11 @@ void Step15::assemble_system () hanging_node_constraints.condense (system_rhs); std::map boundary_values; - // As described above, there is a different boundary condition - // in the first Newton step than in the later ones. This is - // implemented with the help of the bool first_step, which - // will later be false for all times. Starting with the zero- - // function in the first step, we have to set the boundary - // condition $\delta u^{0}=g$ on $\partial \Omega $: - if(first_step) - { - VectorTools::interpolate_boundary_values (dof_handler, - 0, - BoundaryValues(), - boundary_values); - } - // In later steps, the Newton update has to have homogeneous - // boundary conditions, in order for the solution to have the - // right ones. - - else{ VectorTools::interpolate_boundary_values (dof_handler, 0, ZeroFunction(), boundary_values); - } MatrixTools::apply_boundary_values (boundary_values, system_matrix, @@ -337,6 +319,83 @@ void Step15::assemble_system () system_rhs); } + +template +double Step15::compute_residual (const double alpha) const +{ + const QGauss quadrature_formula(3); + + Vector residual (dof_handler.n_dofs()); + + Vector linearization_point (dof_handler.n_dofs()); + linearization_point = present_solution; + linearization_point.add (alpha, newton_update); + + FEValues fe_values (fe, quadrature_formula, + 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(); + + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_rhs = 0; + + fe_values.reinit (cell); + + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { + + // To setup up the linear system, the gradient of the old solution + // in the quadrature points is needed. For this purpose there is + // is a function, which will write these gradients in a vector, + // where every component of the vector is a vector itself: + + std::vector > gradients(n_q_points); + fe_values.get_function_gradients(linearization_point, gradients); + + // Having the gradients of the old solution in the quadrature + // points, we are able to compute the coefficients $a_{n}$ + // in these points. + + const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]); + + // The assembly of the system then is the same as always, except + // of the damping parameter of the Newton method, which we set on + // 0.1 in this case. + + for (unsigned int i = 0; i < dofs_per_cell; ++i) { + cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff + * gradients[q_point] * fe_values.JxW(q_point)); + } + } + + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + for (std::map::const_iterator p = boundary_values.begin(); + p != boundary_values.end(); ++p) + residual(p->first) = 0; + + return residual.l2_norm(); +} + // @sect4{Step15::solve} // The solve function is the same as always, we just have to @@ -346,8 +405,7 @@ void Step15::assemble_system () template void Step15::solve () { - res=system_rhs.l2_norm(); - SolverControl solver_control (1000, res*1e-6); + SolverControl solver_control (1000, system_rhs.l2_norm()*1e-6); SolverMinRes<> solver (solver_control); PreconditionSSOR<> preconditioner; @@ -359,10 +417,17 @@ void Step15::solve () hanging_node_constraints.distribute (newton_update); // In this step, the old solution is updated to the new one: - - present_solution += newton_update; + const double alpha = determine_step_length(); + std::cout << " step length alpha=" << alpha << std::endl; + present_solution.add (alpha, newton_update); } + +template +double Step15::determine_step_length() const +{ + return 0.1; +} // @sect4{Step15::refine_grid} // The first part of this function is the same as in step 6. @@ -429,21 +494,7 @@ void Step15::refine_grid () solution_transfer.interpolate(present_solution,tmp); present_solution=tmp; - // Having refined the mesh, there might be new nodal points on - // the boundary. These have just interpolated values, but - // not the right boundary values. This is fixed up, by - // setting all boundary nodals explicit to the right value: - - std::map boundary_values2; - VectorTools::interpolate_boundary_values(dof_handler, - 0, - BoundaryValues(), - boundary_values2); - for (std::map::const_iterator - p = boundary_values2.begin(); - p != boundary_values2.end(); - ++p) - present_solution(p->first) = p->second; + set_boundary_values (); // On the new mesh, there are different hanging nodes, which shall // be enlisted in a matrix like before. To ensure there are no @@ -461,6 +512,22 @@ void Step15::refine_grid () hanging_node_constraints.distribute(present_solution); } + +template +void Step15::set_boundary_values () +{ + // Having refined the mesh, there might be new nodal points on + // the boundary. These have just interpolated values, but + // not the right boundary values. This is fixed up, by + // setting all boundary nodals explicit to the right value: + + std::map boundary_values2; + VectorTools::interpolate_boundary_values(dof_handler, 0, + BoundaryValues(), boundary_values2); + for (std::map::const_iterator p = + boundary_values2.begin(); p != boundary_values2.end(); ++p) + present_solution(p->first) = p->second; +} // @sect4{Step15::run} // In the run function, the first grid is build. Also in this @@ -489,7 +556,8 @@ void Step15::run () // iteration scheme. Later the Newton method will continue until the // residual is less than $10^{-3}$. - while(first_step || (res>1e-3)) + double previous_res = 0; + while(first_step || (previous_res>1e-3)) { // In the first step, we compute the solution on the two times globally @@ -511,9 +579,13 @@ void Step15::run () setup_system(); + if (first_step) + set_boundary_values (); + // On every mesh there are done five Newton steps, in order to get a // better solution, before the mesh gets too fine and the computations // take more time. + std::cout<<"initial residual:"<::run () // have to be computed. assemble_system (); + previous_res = system_rhs.l2_norm(); + solve (); first_step=false; - std::cout<<"residual:"<