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 241 of file optimization_ex1.C.

References AssembleOptimization::A_matrix, libMesh::DofMap::add_dirichlet_boundary(), libMesh::System::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::System::get_matrix(), libMesh::ImplicitSystem::get_system_matrix(), libMesh::System::get_vector(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::LOCAL_VARIABLE_ORDER, 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::ExodusII_IO::write_equation_systems().

{
  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");
  infile.parse_command_line(argc, argv);

  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.
  ZeroFunction<> zf;

  // Most DirichletBoundary users will want to supply a "locally
  // indexed" functor
  DirichletBoundary dirichlet_bc({0,1,2,3}, {u_var}, 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.get_system_matrix().close();

  // We can get division-by-zero from ILU within the TAO solver, but
  // they can recover from it, so if we have FPEs turned on we should
  // turn them off.
  FPEDisabler disable_fpes;

  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;
}