optimization_ex2.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 396 of file optimization_ex2.C.

References AssembleOptimization::A_matrix, 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(), AssembleOptimization::F_vector, libMesh::System::get_matrix(), libMesh::ImplicitSystem::get_system_matrix(), libMesh::System::get_vector(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), mesh, libMesh::MeshTools::n_elem(), libMesh::OptimizationSystem::optimization_solver, libMesh::out, 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 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");

  // TAO is giving us problems in parallel?
  if (init.comm().size() != 1)
    {
      libMesh::out << "This example can currently only be run in serial." << std::endl;
      return 77;
    }

  GetPot infile("optimization_ex2.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.add_variable("u",
                      Utility::string_to_enum<Order>   (approx_order),
                      Utility::string_to_enum<FEFamily>(fe_family));

  system.optimization_solver->objective_object = &assemble_opt;
  system.optimization_solver->gradient_object = &assemble_opt;
  system.optimization_solver->hessian_object = &assemble_opt;
  system.optimization_solver->equality_constraints_object = &assemble_opt;
  system.optimization_solver->equality_constraints_jacobian_object = &assemble_opt;
  system.optimization_solver->inequality_constraints_object = &assemble_opt;
  system.optimization_solver->inequality_constraints_jacobian_object = &assemble_opt;
  system.optimization_solver->lower_and_upper_bounds_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");


  equation_systems.init();
  equation_systems.print_info();

  assemble_opt.assemble_A_and_F();

  {
    std::vector<std::set<numeric_index_type>> constraint_jac_sparsity;
    std::set<numeric_index_type> sparsity_row;
    sparsity_row.insert(17);
    constraint_jac_sparsity.push_back(sparsity_row);
    sparsity_row.clear();

    sparsity_row.insert(23);
    constraint_jac_sparsity.push_back(sparsity_row);
    sparsity_row.clear();

    sparsity_row.insert(98);
    sparsity_row.insert(185);
    constraint_jac_sparsity.push_back(sparsity_row);

    system.initialize_equality_constraints_storage(constraint_jac_sparsity);
  }

  {
    std::vector<std::set<numeric_index_type>> constraint_jac_sparsity;
    std::set<numeric_index_type> sparsity_row;
    sparsity_row.insert(200);
    sparsity_row.insert(201);
    constraint_jac_sparsity.push_back(sparsity_row);

    system.initialize_inequality_constraints_storage(constraint_jac_sparsity);
  }

  // We need to close the matrix so that we can use it to store the
  // Hessian during the solve.
  system.get_system_matrix().close();
  system.solve();

  // Print convergence information
  system.optimization_solver->print_converged_reason();

  // Write the solution to file if the optimization solver converged
  if (system.optimization_solver->get_converged_reason() > 0)
    {
      std::stringstream filename;
      ExodusII_IO (mesh).write_equation_systems("optimization_soln.exo",
                                                equation_systems);
    }

  return 0;
}