From 8b24976f080ab920810e5d18a02dcc0e2bfda6a5 Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 14 Apr 2004 19:40:21 +0000 Subject: [PATCH] This yields reasonable good results, so let's check it in for the moment. git-svn-id: https://svn.dealii.org/trunk@8999 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-15/step-15.cc | 196 ++++++++++++++-------------- 1 file changed, 100 insertions(+), 96 deletions(-) diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index 2f4e11d66e..00b33b7b93 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -61,15 +61,15 @@ class MinimizationProblem MinimizationProblem (); ~MinimizationProblem (); void run (); + void output_results (const unsigned int cycle) const; private: void setup_system (); - void assemble_step (const bool p); + void assemble_step (); double line_search (const Vector & update) const; void do_step (); void initialize (); void refine_grid (); - void output_results (const unsigned int cycle) const; static double energy (const DoFHandler &dof_handler, const Vector &function); @@ -105,13 +105,19 @@ class InitializationValues : public Function<1> double InitializationValues::value (const Point<1> &p, const unsigned int) const { - return std::pow(p(0), 1.); + const double base = std::pow(p(0), 1./3.); + const double random = 2.*rand()/RAND_MAX-1; + if (base+.1*random < 0 ) + return 0; + else + return base+.1*random; } template -MinimizationProblem::MinimizationProblem () : +MinimizationProblem::MinimizationProblem () + : fe (1), dof_handler (triangulation) {} @@ -157,10 +163,9 @@ double gradient_power (const Tensor<1,dim> &v, template -void MinimizationProblem::assemble_step (const bool p) +void MinimizationProblem::assemble_step () { - if (p) - matrix.reinit (sparsity_pattern); + matrix.reinit (sparsity_pattern); residual.reinit (dof_handler.n_dofs()); QGauss3 quadrature_formula; @@ -207,45 +212,13 @@ void MinimizationProblem::assemble_step (const bool p) for (unsigned int i=0; idiameter() * cell->diameter() + + fe_values.shape_value(i,q_point) * + fe_values.shape_value(j,q_point)) * + fe_values.JxW(q_point); cell_rhs(i) += -((6. * x_minus_u3 * gradient_power (local_solution_grads[q_point], @@ -261,22 +234,21 @@ void MinimizationProblem::assemble_step (const bool p) ) * fe_values.JxW(q_point)); - }; - }; + } + } cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i double MinimizationProblem::line_search (const Vector &update) const { - double alpha = 0.01; - double optimal_energy = energy (dof_handler, present_solution); + double alpha = 0.; Vector tmp (present_solution.size()); - for (double a=.01; a<=10; a*=1.5) + for (unsigned int step=0; step<5; ++step) { tmp = present_solution; - tmp.add (a, update); - const double e = energy(dof_handler, tmp); + tmp.add (alpha, update); + const double f_s = energy (dof_handler, tmp); + + const double dalpha = (alpha != 0 ? alpha/100 : 0.01); + + tmp = present_solution; + tmp.add (alpha+dalpha, update); + const double f_s_plus = energy (dof_handler, tmp); + + tmp = present_solution; + tmp.add (alpha-dalpha, update); + const double f_s_minus = energy (dof_handler, tmp); - std::cout << "XX" << a << ' ' << e << std::endl; - if (e < optimal_energy) + const double f_s_prime = (f_s_plus-f_s_minus) / (2*dalpha); + const double f_s_doubleprime = ((f_s_plus-2*f_s+f_s_minus) / + (dalpha*dalpha)); + + if (std::fabs(f_s_prime) < 1e-7*std::fabs(f_s)) + break; + + if (std::fabs(f_s_doubleprime) < 1e-7*std::fabs(f_s_prime)) + break; + + double step_length = -f_s_prime / f_s_doubleprime; + for (unsigned int i=0; i<3; ++i) { - optimal_energy = e; - alpha = a; + tmp = present_solution; + tmp.add (alpha+step_length, update); + const double e = energy (dof_handler, tmp); + + if (e >= f_s) + step_length /= 2; + else + break; } + alpha += step_length; } - std::cout << " Step length : " << alpha << ' ' << optimal_energy << std::endl; - return alpha; } @@ -333,35 +329,46 @@ MinimizationProblem::line_search (const Vector &update) const template void MinimizationProblem::do_step () { - assemble_step (true); + assemble_step (); Vector update (present_solution.size()); { - SolverControl solver_control (1000, - 1e-3*residual.l2_norm()); + SolverControl solver_control (residual.size(), + 1e-2*residual.l2_norm()); PrimitiveVectorMemory<> vector_memory; SolverCG<> solver (solver_control, vector_memory); PreconditionSSOR<> preconditioner; - preconditioner.initialize(matrix, 1.2); + preconditioner.initialize(matrix); solver.solve (matrix, update, residual, preconditioner); hanging_node_constraints.distribute (update); } - - present_solution.add (line_search (update), update); + + const double step_length = line_search (update); + present_solution.add (step_length, update); } -template -void MinimizationProblem::initialize () +template <> +void MinimizationProblem<1>::initialize () { dof_handler.distribute_dofs (fe); present_solution.reinit (dof_handler.n_dofs()); VectorTools::interpolate (dof_handler, InitializationValues(), present_solution); + DoFHandler<1>::cell_iterator cell; + cell = dof_handler.begin(0); + while (cell->has_children()) + cell = cell->child(0); + present_solution(cell->vertex_dof_index(0,0)) = 0; + + cell = dof_handler.begin(0); + while (cell->has_children()) + cell = cell->child(1); + present_solution(cell->vertex_dof_index(1,0)) = 1; } @@ -441,7 +448,7 @@ MinimizationProblem::energy (const DoFHandler &dof_handler, gradient_power (local_solution_grads[q_point], 6) * fe_values.JxW (q_point)); - }; + } return energy; } @@ -483,36 +490,27 @@ void MinimizationProblem::run () triangulation.refine_global (4); initialize (); - for (unsigned int refinement_cycle=0; refinement_cycle<5; - ++refinement_cycle) + double last_energy = energy (dof_handler, present_solution); + + while (true) { - std::cout << "Cycle " << refinement_cycle << ':' << std::endl; - - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl; - setup_system (); unsigned int iteration=0; for (; iteration<5; ++iteration) - { - do_step (); + do_step (); + + const double this_energy = energy (dof_handler, present_solution); + std::cout << " Energy: " << this_energy << std::endl; + + if ((last_energy-this_energy) < 1e-5*last_energy) + break; + + last_energy = this_energy; - if (residual.l2_norm() < 1.e-4) - break; - }; - output_results (refinement_cycle); - - std::cout << " Iterations : " - << iteration - << std::endl; - std::cout << " Energy : " - << energy (dof_handler, present_solution) - << std::endl; - refine_grid (); - }; + } + std::cout << std::endl; } @@ -522,8 +520,14 @@ int main () { deallog.depth_console (0); - MinimizationProblem<1> minimization_problem_1d; - minimization_problem_1d.run (); + for (unsigned int realization=0; realization<100; ++realization) + { + std::cout << "Realization " << realization << ":" << std::endl; + + MinimizationProblem<1> minimization_problem_1d; + minimization_problem_1d.run (); + minimization_problem_1d.output_results (realization); + } } catch (std::exception &exc) { @@ -547,6 +551,6 @@ int main () << "----------------------------------------------------" << std::endl; return 1; - }; + } return 0; } -- 2.39.5