optimization_ex1.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
// <h1>Optimization Example 1 - Optimization of a quadratic objective function</h1>
// \author David Knezevic
// \date 2015
//
// In this example we demonstrate how to use OptimizationSystem to
// solve the optimization problem:
//   min_U 0.5*U^T A U - U^T F
// We enforce Dirichlet constraints on U so that the minimization is
// well-posed.  But note that we do not use OptimizationSystem's
// interface for imposing constraints in this case, so we can use an
// unconstrained solver (e.g. TAO's "Newton's method with line search"
// for unconstrained optimization is the default choice).
 
// C++ include files that we need
#include <iostream>
 
// libMesh includes
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/elem.h"
#include "libmesh/fpe_disabler.h"
#include "libmesh/boundary_info.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/getpot.h"
#include "libmesh/zero_function.h"
#include "libmesh/dirichlet_boundaries.h"
#include "libmesh/optimization_system.h"
#include "libmesh/optimization_solver.h"
 
// Bring in everything from the libMesh namespace
using namespace libMesh;
 
 
class AssembleOptimization :
  public OptimizationSystem::ComputeObjective,
  public OptimizationSystem::ComputeGradient,
  public OptimizationSystem::ComputeHessian
{
private:
 
  OptimizationSystem & _sys;
 
public:
 
  AssembleOptimization(OptimizationSystem & sys_in);
 
  void assemble_A_and_F();
 
  virtual Number objective (const NumericVector<Number> & soln,
                            OptimizationSystem & /*sys*/);
 
  virtual void gradient (const NumericVector<Number> & soln,
                         NumericVector<Number> & grad_f,
                         OptimizationSystem & /*sys*/);
 
  virtual void hessian (const NumericVector<Number> & soln,
                        SparseMatrix<Number> & H_f,
                        OptimizationSystem & /*sys*/);
 
  SparseMatrix<Number> * A_matrix;
 
  NumericVector<Number> * F_vector;
};
 
AssembleOptimization::AssembleOptimization(OptimizationSystem & sys_in) :
  _sys(sys_in)
{}
 
void AssembleOptimization::assemble_A_and_F()
{
  A_matrix->zero();
  F_vector->zero();
 
  const MeshBase & mesh = _sys.get_mesh();
 
  const unsigned int dim = mesh.mesh_dimension();
  const unsigned int u_var = _sys.variable_number ("u");
 
  const DofMap & dof_map = _sys.get_dof_map();
  FEType fe_type = dof_map.variable_type(u_var);
  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
  QGauss qrule (dim, fe_type.default_quadrature_order());
  fe->attach_quadrature_rule (&qrule);
 
  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<std::vector<Real>> & phi = fe->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
 
  std::vector<dof_id_type> dof_indices;
 
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
 
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      dof_map.dof_indices (elem, dof_indices);
 
      const unsigned int n_dofs = dof_indices.size();
 
      fe->reinit (elem);
 
      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);
 
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          for (unsigned int dof_i=0; dof_i<n_dofs; dof_i++)
            {
              for (unsigned int dof_j=0; dof_j<n_dofs; dof_j++)
                {
                  Ke(dof_i, dof_j) += JxW[qp] * (dphi[dof_j][qp]* dphi[dof_i][qp]);
                }
              Fe(dof_i) += JxW[qp] * phi[dof_i][qp];
            }
        }
 
#ifdef LIBMESH_ENABLE_CONSTRAINTS
      // This will zero off-diagonal entries of Ke corresponding to
      // Dirichlet dofs.
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
 
      // We want the diagonal of constrained dofs to be zero too
      for (unsigned int local_dof_index=0; local_dof_index<n_dofs; local_dof_index++)
        {
          dof_id_type global_dof_index = dof_indices[local_dof_index];
          if (dof_map.is_constrained_dof(global_dof_index))
            {
              Ke(local_dof_index, local_dof_index) = 0.;
            }
        }
#endif // LIBMESH_ENABLE_CONSTRAINTS
 
      A_matrix->add_matrix (Ke, dof_indices);
      F_vector->add_vector (Fe, dof_indices);
    }
 
  A_matrix->close();
  F_vector->close();
}
 
Number AssembleOptimization::objective (const NumericVector<Number> & soln,
                                        OptimizationSystem & /*sys*/)
{
  std::unique_ptr<NumericVector<Number>> AxU = soln.zero_clone();
 
  A_matrix->vector_mult(*AxU, soln);
 
  Number
    UTxAxU = AxU->dot(soln),
    UTxF   = F_vector->dot(soln);
 
  return 0.5 * UTxAxU - UTxF;
}
 
void AssembleOptimization::gradient (const NumericVector<Number> & soln,
                                     NumericVector<Number> & grad_f,
                                     OptimizationSystem & /*sys*/)
{
  grad_f.zero();
 
  // Since we've enforced constraints on soln, A and F,
  // this automatically sets grad_f to zero for constrained
  // dofs.
  A_matrix->vector_mult(grad_f, soln);
  grad_f.add(-1, *F_vector);
}
 
 
void AssembleOptimization::hessian (const NumericVector<Number> & /*soln*/,
                                    SparseMatrix<Number> & H_f,
                                    OptimizationSystem & /*sys*/)
{
  LOG_SCOPE("hessian()", "AssembleOptimization");
  H_f.zero();
  H_f.add(1., *A_matrix);
}
 
 
int main (int argc, char ** argv)
{
  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;
}