I'm interested in having an interface to optimization software from libMesh.
I've attached some code that provides: - A generic optimization solver interface (OptimizationSolver), modeled on NonlinearSolver. - OptimizationSystem, modeled on NonlinearImplicitSystem. - TaoOptimizationSolver, which is a subclass of OptimizationSolver for the specific case of TAO. Note that TAO is included in PETSc now (since version 3.5 I believe). I've implemented a simple test problem in tao_test.C that solves: min 0.5 * U^T A U - U^T. You can run, say, "./example-opt -tao_monitor -tao_view -tao_type nls" and you should get the same solution as if you solved AU = F. If there's broad interest in this, I can clean up the code and turn it into a PR. Some ideas for future additions: - Add support for inequality and equality constraints (inequality constraints are my main motivation for considering optimization problems). I looked into this a bit with TAO and it seems like we have to provide a Jacobian for the constraints, which is a rectangular matrix. - Add interfaces for other optimization suites, e.g. Ipopt or NLOpt David
// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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
// C++ includes
// Local Includes
#include "optimization_solver.h"
#include "tao_optimization_solver.h"
namespace libMesh
{
template <typename T>
inline
OptimizationSolver<T>::OptimizationSolver (sys_type& s) :
ParallelObject(s),
objective_object(NULL),
gradient_object(NULL),
hessian_object(NULL),
equality_constraints(NULL),
_system(s),
_is_initialized (false)
{
}
template <typename T>
inline
OptimizationSolver<T>::~OptimizationSolver ()
{
this->clear ();
}
#if defined(LIBMESH_HAVE_PETSC)
template <typename T>
UniquePtr<OptimizationSolver<T> >
OptimizationSolver<T>::build(sys_type& s, const SolverPackage solver_package)
{
// Build the appropriate solver
switch (solver_package)
{
#ifdef LIBMESH_HAVE_PETSC
case PETSC_SOLVERS:
return UniquePtr<OptimizationSolver<T> >(new TaoOptimizationSolver<T>(s));
#endif // LIBMESH_HAVE_PETSC
default:
libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
}
libmesh_error_msg("We'll never get here!");
return UniquePtr<OptimizationSolver<T> >();
}
#else // LIBMESH_HAVE_PETSC
template <typename T>
UniquePtr<OptimizationSolver<T> >
OptimizationSolver<T>::build(sys_type&, const SolverPackage)
{
libmesh_not_implemented_msg("ERROR: libMesh was compiled without optimization solver support");
}
#endif
//------------------------------------------------------------------
// Explicit instantiations
template class OptimizationSolver<Number>;
} // namespace libMesh
// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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
#ifndef LIBMESH_OPTIMIZATION_SOLVER_H
#define LIBMESH_OPTIMIZATION_SOLVER_H
// Local includes
#include "libmesh/libmesh_common.h"
#include "libmesh/enum_solver_package.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/nonlinear_implicit_system.h"
#include "libmesh/libmesh.h"
#include "libmesh/parallel_object.h"
#include "libmesh/auto_ptr.h"
#include "optimization_system.h"
// C++ includes
#include <cstddef>
namespace libMesh
{
// forward declarations
template <typename T> class SparseMatrix;
template <typename T> class NumericVector;
template <typename T> class Preconditioner;
/**
* This class provides a uniform interface for optimization solvers. This base
* class is overloaded to provide nonlinear solvers from different packages
* like PETSC.
*
* @author David Knezevic, 2015
*/
template <typename T>
class OptimizationSolver : public ReferenceCountedObject<NonlinearSolver<T> >,
public ParallelObject
{
public:
/**
* The type of system
*/
typedef OptimizationSystem sys_type;
/**
* Constructor. Initializes Solver data structures
*/
explicit
OptimizationSolver (sys_type& s);
/**
* Destructor.
*/
virtual ~OptimizationSolver ();
/**
* Builds a \p NonlinearSolver using the nonlinear solver package specified by
* \p solver_package
*/
static UniquePtr<OptimizationSolver<T> > build(sys_type& s,
const SolverPackage solver_package = libMesh::default_solver_package());
/**
* @returns true if the data structures are
* initialized, false otherwise.
*/
bool initialized () const { return _is_initialized; }
/**
* Release all memory and clear data structures.
*/
virtual void clear () {}
/**
* Initialize data structures if not done so already.
*/
virtual void init () = 0;
/**
* Solves the optimization problem.
*/
virtual void solve () = 0;
/**
* Prints a useful message about why the latest nonlinear solve
* con(di)verged.
*/
virtual void print_converged_reason() { libmesh_not_implemented(); }
/**
* Object that computes the objective function f(X)
* at the input iterate X.
*/
OptimizationSystem::ComputeObjective *objective_object;
/**
* Object that computes the gradient grad_f(X) of the objective function
* at the input iterate X.
*/
OptimizationSystem::ComputeGradient *gradient_object;
/**
* Object that computes the Hessian H_f(X) of the objective function
* at the input iterate X.
*/
OptimizationSystem::ComputeHessian *hessian_object;
/**
* Object that computes the equality constraints vector C_eq(X).
* This will lead to the constraints C_eq(X) = 0 being imposed.
*/
OptimizationSystem::ComputeEqualityConstraints *equality_constraints;
/**
* @returns a constant reference to the system we are using to
* define the optimization problem.
*/
const sys_type & system () const { return _system; }
/**
* @returns a writeable reference to the system we are using to
* define the optimization problem.
*/
sys_type & system () { return _system; }
protected:
/**
* A reference to the system we are solving.
*/
sys_type& _system;
/**
* Flag indicating if the data structures have been initialized.
*/
bool _is_initialized;
};
} // namespace libMesh
#endif // LIBMESH_OPTIMIZATION_SOLVER_H
// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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
// C++ includes
// Local includes
#include "libmesh/equation_systems.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "optimization_solver.h"
#include "optimization_system.h"
namespace libMesh
{
// ------------------------------------------------------------
// OptimizationSystem implementation
OptimizationSystem::OptimizationSystem (EquationSystems& es,
const std::string& name_in,
const unsigned int number_in) :
Parent(es, name_in, number_in),
optimization_solver(OptimizationSolver<Number>::build(*this)),
C_eq(NumericVector<Number>::build(this->comm()))
{
}
OptimizationSystem::~OptimizationSystem ()
{
// Clear data
this->clear();
}
void OptimizationSystem::clear ()
{
// clear the optimization solver
optimization_solver->clear();
// clear the parent data
Parent::clear();
}
void OptimizationSystem::reinit ()
{
optimization_solver->clear();
Parent::reinit();
}
void OptimizationSystem::initialize_equality_constraints_storage(
unsigned int n_eq_constraints)
{
C_eq->init(n_eq_constraints, false, SERIAL);
}
void OptimizationSystem::solve ()
{
// Log how long the optimization solve takes.
START_LOG("solve()", "OptimizationSystem");
optimization_solver->init();
// Solve the optimization problem.
optimization_solver->solve ();
// Stop logging the nonlinear solve
STOP_LOG("solve()", "OptimizationSystem");
// Update the system after the solve
this->update();
}
} // namespace libMesh
// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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
#ifndef LIBMESH_OPTIMIZATION_SYSTEM_H
#define LIBMESH_OPTIMIZATION_SYSTEM_H
// Local Includes
#include "libmesh/implicit_system.h"
// C++ includes
namespace libMesh
{
// Forward declarations
template<typename T> class OptimizationSolver;
/**
* This System subclass enables us to assemble an objective function,
* gradient, Hessian and bounds for optimization problems.
*/
// ------------------------------------------------------------
// OptimizationSystem class definition
class OptimizationSystem : public ImplicitSystem
{
public:
/**
* Constructor. Optionally initializes required
* data structures.
*/
OptimizationSystem (EquationSystems& es,
const std::string& name,
const unsigned int number);
/**
* Destructor.
*/
virtual ~OptimizationSystem ();
/**
* The type of system.
*/
typedef OptimizationSystem sys_type;
/**
* The type of the parent.
*/
typedef ImplicitSystem Parent;
/**
* Abstract base class to be used to calculate the objective
* function for optimization.
*/
class ComputeObjective
{
public:
virtual ~ComputeObjective () {}
/**
* This function will be called to compute the objective function
* to be minimized, and must be implemented by the user in a
* derived class. @return the value of the objective function at
* the iterate \p X.
*/
virtual Number objective (const NumericVector<Number>& X,
sys_type& S) = 0;
};
/**
* Abstract base class to be used to calculate the gradient of
* an objective function.
*/
class ComputeGradient
{
public:
virtual ~ComputeGradient () {}
/**
* This function will be called to compute the gradient of the
* objective function, and must be implemented by the user in
* a derived class. Set \p grad_f to be the gradient at the
* iterate \p X.
*/
virtual void gradient (const NumericVector<Number>& X,
NumericVector<Number>& grad_f,
sys_type& S) = 0;
};
/**
* Abstract base class to be used to calculate the Hessian
* of an objective function.
*/
class ComputeHessian
{
public:
virtual ~ComputeHessian () {}
/**
* This function will be called to compute the Hessian of
* the objective function, and must be implemented by the
* user in a derived class. Set \p H_f to be the gradient
* at the iterate \p X.
*/
virtual void hessian (const NumericVector<Number>& X,
SparseMatrix<Number>& H_f,
sys_type& S) = 0;
};
/**
* Abstract base class to be used to calculate the equality constraints.
*/
class ComputeEqualityConstraints
{
public:
virtual ~ComputeEqualityConstraints () {}
/**
* This function will be called to evaluate the equality constraints
* vector C_eq(X). This will impose the constraints C_eq(X) = 0.
*/
virtual void equality_constraints (const NumericVector<Number>& X,
NumericVector<Number>& C_eq,
sys_type& S) = 0;
};
/**
* Initialize the vector and matrix that store the equality constraint
* and corresponding Jacobian.
*/
void initialize_equality_constraint_storage(
unsigned int n_eq_constraints,
std::vector<unsigned int> n_dofs_per_constraint);
/**
* @returns a clever pointer to the system.
*/
sys_type & system () { return *this; }
/**
* Clear all the data structures associated with
* the system.
*/
virtual void clear ();
/**
* Reinitializes the member data fields associated with
* the system, so that, e.g., \p assemble() may be used.
*/
virtual void reinit ();
/**
* Assembles & solves the nonlinear system R(x) = 0.
*/
virtual void solve ();
/**
* Initialize storage for the \p n_eq_constraints
* equality constraints.
*/
void initialize_equality_constraints_storage(
unsigned int n_eq_constraints);
/**
* @returns \p "Optimization". Helps in identifying
* the system type in an equation system file.
*/
virtual std::string system_type () const { return "Optimization"; }
/**
* The \p OptimizationSolver that is used for performing the optimization.
*/
UniquePtr<OptimizationSolver<Number> > optimization_solver;
/**
* The vector that stores equality constraints.
*/
UniquePtr<NumericVector<Number> > C_eq;
};
} // namespace libMesh
#endif // LIBMESH_OPTIMIZATION_SYSTEM_H
// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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
#include "libmesh/libmesh_common.h"
#ifdef LIBMESH_HAVE_PETSC
// C++ includes
// Local Includes
#include "libmesh/petsc_vector.h"
#include "libmesh/petsc_matrix.h"
#include "libmesh/dof_map.h"
#include "tao_optimization_solver.h"
namespace libMesh
{
//--------------------------------------------------------------------
// Functions with C linkage to pass to Tao. Tao will call these
// methods as needed.
//
// Since they must have C linkage they have no knowledge of a namespace.
// Give them an obscure name to avoid namespace pollution.
extern "C"
{
//---------------------------------------------------------------
// This function is called by Tao to evaluate objective function at x
PetscErrorCode
__libmesh_tao_objective (Tao tao, Vec x, PetscReal* objective, void *ctx)
{
START_LOG("objective()", "TaoOptimizationSolver");
PetscErrorCode ierr = 0;
libmesh_assert(x);
libmesh_assert(objective);
libmesh_assert(ctx);
// ctx should be a pointer to the solver (it was passed in as void*)
TaoOptimizationSolver<Number>* solver =
static_cast<TaoOptimizationSolver<Number>*> (ctx);
OptimizationSystem &sys = solver->system();
// We'll use current_local_solution below, so let's ensure that it's consistent
// with the vector x that was passed in.
PetscVector<Number>& X_sys = *cast_ptr<PetscVector<Number>*>(sys.solution.get());
PetscVector<Number> X(x, sys.comm());
// Perform a swap so that sys.solution points to X
X.swap(X_sys);
// Impose constraints on X
sys.get_dof_map().enforce_constraints_exactly(sys);
// Update sys.current_local_solution based on X
sys.update();
// Swap back
X.swap(X_sys);
if (solver->objective_object != NULL)
{
(*objective) =
solver->objective_object->objective(
*(sys.current_local_solution), sys);
}
else
{
libmesh_error_msg("Objective function not defined in __libmesh_tao_objective");
}
STOP_LOG("objective()", "TaoOptimizationSolver");
return ierr;
}
//---------------------------------------------------------------
// This function is called by Tao to evaluate the gradient at x
PetscErrorCode
__libmesh_tao_gradient(Tao tao, Vec x, Vec g, void *ctx)
{
START_LOG("gradient()", "TaoOptimizationSolver");
PetscErrorCode ierr = 0;
libmesh_assert(x);
libmesh_assert(g);
libmesh_assert(ctx);
// ctx should be a pointer to the solver (it was passed in as void*)
TaoOptimizationSolver<Number>* solver =
static_cast<TaoOptimizationSolver<Number>*> (ctx);
OptimizationSystem &sys = solver->system();
// We'll use current_local_solution below, so let's ensure that it's consistent
// with the vector x that was passed in.
PetscVector<Number>& X_sys = *cast_ptr<PetscVector<Number>*>(sys.solution.get());
PetscVector<Number> X(x, sys.comm());
// Perform a swap so that sys.solution points to X
X.swap(X_sys);
// Impose constraints on X
sys.get_dof_map().enforce_constraints_exactly(sys);
// Update sys.current_local_solution based on X
sys.update();
// Swap back
X.swap(X_sys);
// We'll also pass the gradient in to the assembly routine
// so let's make a PETSc vector for that too.
PetscVector<Number> gradient(g, sys.comm());
// Clear the gradient prior to assembly
gradient.zero();
if (solver->gradient_object != NULL)
{
solver->gradient_object->gradient(
*(sys.current_local_solution), gradient, sys);
}
else
{
libmesh_error_msg("Gradient function not defined in __libmesh_tao_gradient");
}
gradient.close();
STOP_LOG("gradient()", "TaoOptimizationSolver");
return ierr;
}
//---------------------------------------------------------------
// This function is called by Tao to evaluate the Hessian at x
PetscErrorCode
__libmesh_tao_hessian(Tao tao, Vec x, Mat h, Mat pc, void *ctx)
{
START_LOG("hessian()", "TaoOptimizationSolver");
PetscErrorCode ierr = 0;
libmesh_assert(x);
libmesh_assert(hessian);
libmesh_assert(pc);
libmesh_assert(ctx);
// ctx should be a pointer to the solver (it was passed in as void*)
TaoOptimizationSolver<Number>* solver =
static_cast<TaoOptimizationSolver<Number>*> (ctx);
OptimizationSystem &sys = solver->system();
// We'll use current_local_solution below, so let's ensure that it's consistent
// with the vector x that was passed in.
PetscVector<Number>& X_sys = *cast_ptr<PetscVector<Number>*>(sys.solution.get());
PetscVector<Number> X(x, sys.comm());
// Perform a swap so that sys.solution points to X
X.swap(X_sys);
// Impose constraints on X
sys.get_dof_map().enforce_constraints_exactly(sys);
// Update sys.current_local_solution based on X
sys.update();
// Swap back
X.swap(X_sys);
// Let's also set up matrices that will be
PetscMatrix<Number> PC(pc, sys.comm());
PetscMatrix<Number> hessian(h, sys.comm());
PC.attach_dof_map(sys.get_dof_map());
hessian.attach_dof_map(sys.get_dof_map());
if (solver->hessian_object != NULL)
{
// Following PetscNonlinearSolver by passing in PC. It's not clear
// why we pass in PC and not hessian though?
solver->hessian_object->hessian(
*(sys.current_local_solution), PC, sys);
}
else
{
libmesh_error_msg("Hessian function not defined in __libmesh_tao_hessian");
}
PC.close();
hessian.close();
STOP_LOG("hessian()", "TaoOptimizationSolver");
return ierr;
}
//---------------------------------------------------------------
// This function is called by Tao to evaluate the equality constraints at x
PetscErrorCode
__libmesh_tao_equality_constraints(Tao tao, Vec x, Vec ce, void *ctx)
{
START_LOG("equality_constraints()", "TaoOptimizationSolver");
PetscErrorCode ierr = 0;
libmesh_assert(x);
libmesh_assert(ce);
libmesh_assert(ctx);
// ctx should be a pointer to the solver (it was passed in as void*)
TaoOptimizationSolver<Number>* solver =
static_cast<TaoOptimizationSolver<Number>*> (ctx);
OptimizationSystem &sys = solver->system();
// We'll use current_local_solution below, so let's ensure that it's consistent
// with the vector x that was passed in.
PetscVector<Number>& X_sys = *cast_ptr<PetscVector<Number>*>(sys.solution.get());
PetscVector<Number> X(x, sys.comm());
// Perform a swap so that sys.solution points to X
X.swap(X_sys);
// Impose constraints on X
sys.get_dof_map().enforce_constraints_exactly(sys);
// Update sys.current_local_solution based on X
sys.update();
// Swap back
X.swap(X_sys);
// We'll also pass the gradient in to the assembly routine
// so let's make a PETSc vector for that too.
PetscVector<Number> eq_constraints(ce, sys.comm());
// Clear the gradient prior to assembly
eq_constraints.zero();
if (solver->equality_constraints != NULL)
{
solver->equality_constraints->equality_constraints(
*(sys.current_local_solution), eq_constraints, sys);
}
else
{
libmesh_error_msg("Constraints function not defined in __libmesh_tao_equality_constraints");
}
eq_constraints.close();
STOP_LOG("equality_constraints()", "TaoOptimizationSolver");
return ierr;
}
} // end extern "C"
//---------------------------------------------------------------------
//---------------------------------------------------------------------
// TaoOptimizationSolver<> methods
template <typename T>
TaoOptimizationSolver<T>::TaoOptimizationSolver (OptimizationSystem& system_in)
:
OptimizationSolver<T>(system_in),
_reason(TAO_CONVERGED_FATOL/*==0*/) // Arbitrary initial value...
{
}
template <typename T>
TaoOptimizationSolver<T>::~TaoOptimizationSolver ()
{
this->clear ();
}
template <typename T>
void TaoOptimizationSolver<T>::clear ()
{
if (this->initialized())
{
this->_is_initialized = false;
PetscErrorCode ierr=0;
ierr = TaoDestroy(&_tao);
LIBMESH_CHKERRABORT(ierr);
}
}
template <typename T>
void TaoOptimizationSolver<T>::init ()
{
// Initialize the data structures if not done so already.
if (!this->initialized())
{
this->_is_initialized = true;
PetscErrorCode ierr=0;
ierr = TaoCreate(this->comm().get(),&_tao);
LIBMESH_CHKERRABORT(ierr);
}
}
template <typename T>
void TaoOptimizationSolver<T>::solve ()
{
START_LOG("solve()", "TaoOptimizationSolver");
this->init ();
this->system().solution->zero();
PetscMatrix<T>* hessian = cast_ptr<PetscMatrix<T>*>(this->system().matrix);
PetscVector<T>* gradient = cast_ptr<PetscVector<T>*>(this->system().rhs);
PetscVector<T>* x = cast_ptr<PetscVector<T>*>(this->system().solution.get());
PetscVector<T>* ce = cast_ptr<PetscVector<T>*>(this->system().C_eq.get());
// Set the starting guess to zero.
x->zero();
PetscErrorCode ierr = 0;
// Set solution vec and an initial guess
ierr = TaoSetInitialVector(_tao, x->vec());
LIBMESH_CHKERRABORT(ierr);
// We have to have an objective function
libmesh_assert( this->objective_object );
// Set routines for objective, gradient, hessian evaluation
ierr = TaoSetObjectiveRoutine(_tao, __libmesh_tao_objective, this);
LIBMESH_CHKERRABORT(ierr);
if ( this->gradient_object )
{
ierr = TaoSetGradientRoutine(_tao, __libmesh_tao_gradient, this);
LIBMESH_CHKERRABORT(ierr);
}
if ( this->hessian_object )
{
ierr = TaoSetHessianRoutine(_tao, hessian->mat(), hessian->mat(), __libmesh_tao_hessian, this);
LIBMESH_CHKERRABORT(ierr);
}
// Optionally set equality constraints
if ( this->equality_constraints )
{
ierr = TaoSetEqualityConstraintsRoutine(_tao, ce->vec(), __libmesh_tao_equality_constraints, this);
LIBMESH_CHKERRABORT(ierr);
}
// Check for Tao command line options
ierr = TaoSetFromOptions(_tao);
LIBMESH_CHKERRABORT(ierr);
// Perform the optimization
ierr = TaoSolve(_tao);
LIBMESH_CHKERRABORT(ierr);
// Store the convergence/divergence reason
ierr = TaoGetConvergedReason(_tao, &_reason);
LIBMESH_CHKERRABORT(ierr);
// Print termination information
libMesh::out << "Converged reason: " << _reason << std::endl;
if (_reason <= 0)
{
libMesh::out << "Tao failed to converge." << std::endl;
}
else
{
libMesh::out << "Tao converged." << std::endl;
}
STOP_LOG("solve()", "TaoOptimizationSolver");
}
//template <typename T>
//void TaoOptimizationSolver<T>::print_converged_reason()
//{
//
// libMesh::out << "Tao optimization solver convergence/divergence reason: "
// << TaoConvergedReasons[this->get_converged_reason()] << std::endl;
//}
template <typename T>
TaoConvergedReason TaoOptimizationSolver<T>::get_converged_reason()
{
PetscErrorCode ierr=0;
if (this->initialized())
{
ierr = TaoGetConvergedReason(_tao, &_reason);
LIBMESH_CHKERRABORT(ierr);
}
return _reason;
}
//------------------------------------------------------------------
// Explicit instantiations
template class TaoOptimizationSolver<Number>;
} // namespace libMesh
#endif // #ifdef LIBMESH_HAVE_PETSC
// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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
#ifndef LIBMESH_TAO_OPTIMIZATION_SOLVER_H
#define LIBMESH_TAO_OPTIMIZATION_SOLVER_H
#include "libmesh/libmesh_config.h"
// Petsc include files.
#ifdef LIBMESH_HAVE_PETSC
// Local includes
#include "libmesh/petsc_macro.h"
// Include header for the Tao optimization library
EXTERN_C_FOR_PETSC_BEGIN
# include <petsctao.h>
EXTERN_C_FOR_PETSC_END
#include "optimization_solver.h"
// C++ includes
namespace libMesh
{
// Allow users access to these functions in case they want to reuse them. Note that users shouldn't
// need access to these most of the time as they are used internally by this object.
extern "C"
{
PetscErrorCode __libmesh_tao_objective (Tao tao, Vec x, PetscReal* objective, void *ctx);
PetscErrorCode __libmesh_tao_gradient(Tao tao, Vec x, Vec g, void *ctx);
PetscErrorCode __libmesh_tao_hessian(Tao tao, Vec x, Mat h, Mat pc, void *ctx);
PetscErrorCode __libmesh_tao_equality_constraints(Tao tao, Vec x, Vec ce, void *ctx);
}
/**
* This class provides an interface to the Tao optimization solvers.
*
* @author David Knezevic, 2015
*/
template <typename T>
class TaoOptimizationSolver : public OptimizationSolver<T>
{
public:
/**
* The type of system that we use in conjunction with this solver.
*/
typedef OptimizationSystem sys_type;
/**
* Constructor. Initializes Tao data structures.
*/
explicit
TaoOptimizationSolver (sys_type& system);
/**
* Destructor.
*/
~TaoOptimizationSolver ();
/**
* Release all memory and clear data structures.
*/
virtual void clear ();
/**
* Initialize data structures if not done so already.
*/
virtual void init ();
/**
* Returns the raw PETSc Tao context pointer.
*/
Tao tao() { this->init(); return _tao; }
/**
* Call the Tao solver.
*/
virtual void solve ();
// /**
// * Prints a useful message about why the latest nonlinear solve
// * con(di)verged.
// */
// virtual void print_converged_reason();
/**
* Returns the currently-available (or most recently obtained, if the Tao object has
* been destroyed) convergence reason. Refer to Tao docs for the meaning of different
* TaoConvergedReason.
*/
TaoConvergedReason get_converged_reason();
protected:
/**
* Optimization solver context
*/
Tao _tao;
/**
* Store the reason for Tao convergence/divergence for use even after _tao
* has been cleared. Note that print_converged_reason() will always *try* to
* get the current reason with TaoGetConvergedReason(), but if the Tao object
* has already been cleared, it will fall back on this stored value. Note that
* this value is therefore necessarily *not* cleared by the clear() function.
*/
TaoConvergedReason _reason;
private:
friend PetscErrorCode __libmesh_tao_objective (Tao tao, Vec x, PetscReal* objective, void *ctx);
friend PetscErrorCode __libmesh_tao_gradient(Tao tao, Vec x, Vec g, void *ctx);
friend PetscErrorCode __libmesh_tao_hessian(Tao tao, Vec x, Mat h, Mat pc, void *ctx);
friend PetscErrorCode __libmesh_tao_equality_constraints(Tao tao, Vec x, Vec ce, void *ctx);
};
} // namespace libMesh
#endif // #ifdef LIBMESH_HAVE_PETSC
#endif // LIBMESH_TAO_OPTIMIZATION_SOLVER_H
// 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/boundary_info.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/getpot.h"
#include "libmesh/zero_function.h"
#include "libmesh/dirichlet_boundaries.h"
#include "optimization_system.h"
#include "optimization_solver.h"
// Bring in everything from the libMesh namespace
using namespace libMesh;
/**
* This class encapsulate all functionality required for assembling
* and solving a linear elastic model with contact.
*/
class AssembleOptimization :
public OptimizationSystem::ComputeObjective,
public OptimizationSystem::ComputeGradient,
public OptimizationSystem::ComputeHessian
{
private:
/**
* Keep a reference to the OptimizationSystem.
*/
OptimizationSystem& _sys;
public:
/**
* Constructor.
*/
AssembleOptimization(
OptimizationSystem &sys_in);
/**
* The optimization problem we consider here is:
* min_U 0.5*U^T A U - U^T F.
* In this method, we assemble A and F.
*/
void assemble_A_and_F();
/**
* Evaluate the objective function.
*/
virtual Number objective (
const NumericVector<Number>& soln,
OptimizationSystem& /*sys*/);
/**
* Evaluate the gradient.
*/
virtual void gradient (
const NumericVector<Number>& soln,
NumericVector<Number>& grad_f,
OptimizationSystem& /*sys*/);
/**
* Evaluate the Hessian.
*/
virtual void hessian (
const NumericVector<Number>& soln,
SparseMatrix<Number>& H_f,
OptimizationSystem& /*sys*/);
// /**
// * Evaluate the vector C_eq(x) that defines equality
// * constraints according to C_eq(x) = 0.
// */
// virtual void equality_constraints (
// const NumericVector<Number>& soln,
// NumericVector<Number>& C_eq,
// OptimizationSystem& /*sys*/);
/**
* Sparse matrix for storing the matrix A. We use
* this to facilitate computation of objective, gradient
* and hessian.
*/
SparseMatrix<Number>* A_matrix;
/**
* Vector for storing F. We use this to facilitate
* computation of objective, gradient and hessian.
*/
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);
UniquePtr<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;
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
const Elem* elem = *el;
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];
}
}
// 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.;
}
}
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*/)
{
UniquePtr< NumericVector<Number> > AxU = soln.zero_clone();
A_matrix->vector_mult(*AxU, soln);
Real UTxAxU = AxU->dot(soln);
Real 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 constaints 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*/)
{
H_f.zero();
H_f.add(1., *A_matrix);
}
//void AssembleOptimization::equality_constraints (
// const NumericVector<Number>& soln,
// NumericVector<Number>& C_eq,
// OptimizationSystem& /*sys*/)
//{
// dof_id_type count = 0;
// for(dof_id_type dof_index=0; dof_index<_sys.n_dofs(); dof_index++)
// {
// if(_sys.get_dof_map().is_constrained_dof(dof_index))
// {
// C_eq.set(count, soln(dof_index));
// count++;
// }
// }
//}
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
GetPot infile("tao_test.in");
const std::string approx_order = infile("approx_order", "FIRST");
const std::string fe_family = infile("fe_family", "LAGRANGE");
Mesh mesh(init.comm());
MeshTools::Generation::build_square (mesh,
10,
10,
-1., 1.,
-1., 1.,
QUAD9);
mesh.print_info();
EquationSystems equation_systems (mesh);
OptimizationSystem& system =
equation_systems.add_system<OptimizationSystem> ("Optimization");
AssembleOptimization assemble_opt(system);
unsigned int u_var = 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;
// Seems like -tao_type ipm is the only solver that can handle
// equality constraints.
// system.optimization_solver->equality_constraints = &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");
// Apply Dirichlet constraints. We can then use DofMap::is_constrained_dof
// to figure out which dofs we want to apply equality constraints to.
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;
DirichletBoundary dirichlet_bc(boundary_ids,
variables,
&zf);
system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
equation_systems.init();
equation_systems.print_info();
assemble_opt.assemble_A_and_F();
// Let's assume that we impose an equality constraint
// for every constrained dof in the system
system.initialize_equality_constraints_storage(system.n_constrained_dofs());
// We need to close the matrix so that we can use it to store the
// Hessian during the solve.
system.matrix->close();
system.solve();
std::stringstream filename;
ExodusII_IO (mesh).write_equation_systems(
"optimization_soln.exo",
equation_systems);
return 0;
}
tao_test.in
Description: Binary data
------------------------------------------------------------------------------ BPM Camp - Free Virtual Workshop May 6th at 10am PDT/1PM EDT Develop your own process in accordance with the BPMN 2 standard Learn Process modeling best practices with Bonita BPM through live exercises http://www.bonitasoft.com/be-part-of-it/events/bpm-camp-virtual- event?utm_ source=Sourceforge_BPM_Camp_5_6_15&utm_medium=email&utm_campaign=VA_SF
_______________________________________________ Libmesh-devel mailing list Libmesh-devel@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-devel
