optimization_ex1.C File Reference

Go to the source code of this file.

Classes
class	AssembleOptimization
	This class encapsulate all functionality required for assembling the objective function, gradient, and hessian. More...

Functions
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 239 of file optimization_ex1.C.

{
  LibMeshInit init (argc, argv);
 
#ifndef LIBMESH_HAVE_PETSC_TAO
 
  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support");
 
#elif !defined(LIBMESH_ENABLE_GHOSTED)
 
  libmesh_example_requires(false, "--enable-ghosted");
 
#elif LIBMESH_USE_COMPLEX_NUMBERS
 
  // According to
  // http://www.mcs.anl.gov/research/projects/tao/documentation/installation.html
  // TAO & PETSc-complex are currently mutually exclusive
  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support & real-numbers only");
 
#endif
 
  // We use a 2D domain.
  libmesh_example_requires(LIBMESH_DIM > 1, "--disable-1D-only");
 
  // We use Dirichlet boundary conditions here
#ifndef LIBMESH_ENABLE_DIRICHLET
  libmesh_example_requires(false, "--enable-dirichlet");
#endif
 
  if (libMesh::on_command_line ("--use-eigen"))
    {
      libMesh::err << "This example requires an OptimizationSolver, and therefore does not "
                   << "support --use-eigen on the command line."
                   << std::endl;
      return 0;
    }
 
  GetPot infile("optimization_ex1.in");
  const std::string approx_order = infile("approx_order", "FIRST");
  const std::string fe_family = infile("fe_family", "LAGRANGE");
  const unsigned int n_elem = infile("n_elem", 10);
 
  Mesh mesh(init.comm());
  MeshTools::Generation::build_square (mesh,
                                       n_elem,
                                       n_elem,
                                       -1., 1.,
                                       -1., 1.,
                                       QUAD9);
 
  mesh.print_info();
 
  EquationSystems equation_systems (mesh);
 
  OptimizationSystem & system =
    equation_systems.add_system<OptimizationSystem> ("Optimization");
 
  // The default is to use PETSc/Tao solvers, but let the user change
  // the optimization solver package on the fly.
  {
    const std::string optimization_solver_type = infile("optimization_solver_type",
                                                        "PETSC_SOLVERS");
    SolverPackage sp = Utility::string_to_enum<SolverPackage>(optimization_solver_type);
    std::unique_ptr<OptimizationSolver<Number>> new_solver =
      OptimizationSolver<Number>::build(system, sp);
    system.optimization_solver.reset(new_solver.release());
  }
 
  // Set tolerances and maximum iteration counts directly on the optimization solver.
  system.optimization_solver->max_objective_function_evaluations = 128;
  system.optimization_solver->objective_function_relative_tolerance = 1.e-4;
  system.optimization_solver->verbose = true;
 
  AssembleOptimization assemble_opt(system);
 
  system.optimization_solver->objective_object = &assemble_opt;
  system.optimization_solver->gradient_object  = &assemble_opt;
  system.optimization_solver->hessian_object   = &assemble_opt;
 
  // system.matrix and system.rhs are used for the gradient and Hessian,
  // so in this case we add an extra matrix and vector to store A and F.
  // This makes it easy to write the code for evaluating the objective,
  // gradient, and hessian.
  system.add_matrix("A_matrix");
  system.add_vector("F_vector");
  assemble_opt.A_matrix = &system.get_matrix("A_matrix");
  assemble_opt.F_vector = &system.get_vector("F_vector");
 
#ifdef LIBMESH_ENABLE_DIRICHLET
  unsigned int u_var = system.add_variable("u",
                                           Utility::string_to_enum<Order>   (approx_order),
                                           Utility::string_to_enum<FEFamily>(fe_family));
 
  // Apply Dirichlet constraints. This will be used to apply constraints
  // to the objective function, gradient and Hessian.
  std::set<boundary_id_type> boundary_ids;
  boundary_ids.insert(0);
  boundary_ids.insert(1);
  boundary_ids.insert(2);
  boundary_ids.insert(3);
  std::vector<unsigned int> variables;
  variables.push_back(u_var);
  ZeroFunction<> zf;
 
  // Most DirichletBoundary users will want to supply a "locally
  // indexed" functor
  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
                                 LOCAL_VARIABLE_ORDER);
  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
#endif // LIBMESH_ENABLE_DIRICHLET
 
  equation_systems.init();
  equation_systems.print_info();
 
  assemble_opt.assemble_A_and_F();
 
  // We need to close the matrix so that we can use it to store the
  // Hessian during the solve.
  system.matrix->close();
  system.solve();
 
  // Print convergence information
  system.optimization_solver->print_converged_reason();
 
#ifdef LIBMESH_HAVE_EXODUS_API
  std::stringstream filename;
  ExodusII_IO (mesh).write_equation_systems("optimization_soln.exo",
                                            equation_systems);
#endif
 
  return 0;
}

References AssembleOptimization::A_matrix, libMesh::DofMap::add_dirichlet_boundary(), libMesh::ImplicitSystem::add_matrix(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::System::add_vector(), AssembleOptimization::assemble_A_and_F(), libMesh::MeshTools::Generation::build_square(), libMesh::SparseMatrix< T >::close(), libMesh::err, AssembleOptimization::F_vector, libMesh::System::get_dof_map(), libMesh::ImplicitSystem::get_matrix(), libMesh::System::get_vector(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::LOCAL_VARIABLE_ORDER, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshTools::n_elem(), libMesh::on_command_line(), libMesh::OptimizationSystem::optimization_solver, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::OptimizationSystem::solve(), and libMesh::MeshOutput< MT >::write_equation_systems().