gross-pitaevskii/example_02.hpp

Go to the documentation of this file.

// @HEADER
// ************************************************************************
//
//               Rapid Optimization Library (ROL) Package
//                 Copyright (2014) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact lead developers:
//              Drew Kouri   (dpkouri@sandia.gov) and
//              Denis Ridzal (dridzal@sandia.gov)
//
// ************************************************************************
// @HEADER

#include <iostream>

#include "ROL_Stream.hpp"
#include "Teuchos_GlobalMPISession.hpp"

#include "ROL_ParameterList.hpp"
#include "ROL_StdVector.hpp"
#include "ROL_Objective.hpp"
#include "ROL_Constraint.hpp"
#include "ROL_Algorithm.hpp"
#include "ROL_CompositeStep.hpp"
#include "ROL_ConstraintStatusTest.hpp"

#include "numerics/FiniteDifference.hpp"


using namespace ROL;

template <class Real, class Element=Real>
class OptStdVector;  // Optimization space.

template <class Real, class Element=Real>
class OptDualStdVector;  // Dual optimization space.

template <class Real, class Element=Real>
class ConStdVector;  // Constraint space.

template <class Real, class Element=Real>
class ConDualStdVector;  // Dual constraint space.

// Vector space definitions:


// Optimization space.
template <class Real, class Element>
class OptStdVector : public Vector<Real> {

  typedef std::vector<Element>       vector;
  typedef typename vector::size_type uint;

private:
ROL::Ptr<std::vector<Element> >  std_vec_;
mutable ROL::Ptr<OptDualStdVector<Real> >  dual_vec_;

ROL::Ptr<FiniteDifference<Real> > fd_;


public:

OptStdVector(const ROL::Ptr<std::vector<Element> > & std_vec, ROL::Ptr<FiniteDifference<Real> >fd) : 
    std_vec_(std_vec), dual_vec_(ROL::nullPtr), fd_(fd) {}

void plus( const Vector<Real> &x ) {
  const OptStdVector &ex = dynamic_cast<const OptStdVector&>(x);
  ROL::Ptr<const vector> xvalptr = ex.getVector();
  uint dimension  = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] += (*xvalptr)[i];
  }
}

void scale( const Real alpha ) {
  uint dimension = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] *= alpha;
  }
}


Real dot( const Vector<Real> &x ) const {
  Real val = 0;
  const OptStdVector<Real, Element> & ex = dynamic_cast<const OptStdVector&>(x);
  ROL::Ptr<const vector> xvalptr = ex.getVector();
   
  ROL::Ptr<vector> kxvalptr = ROL::makePtr<vector>(std_vec_->size(), 0.0);

  fd_->apply(xvalptr,kxvalptr);

  uint dimension  = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    val += (*std_vec_)[i]*(*kxvalptr)[i];
  }
  return val;
}

Real norm() const {
  Real val = 0;
  val = std::sqrt( dot(*this) );
  return val;
}

ROL::Ptr<Vector<Real> > clone() const {
  return ROL::makePtr<OptStdVector>( ROL::makePtr<vector>(std_vec_->size()), fd_ );
}

ROL::Ptr<const vector> getVector() const {
  return std_vec_;
}

ROL::Ptr<vector> getVector()  {
  return std_vec_;
}

ROL::Ptr<Vector<Real> > basis( const int i ) const {
  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
  ROL::Ptr<OptStdVector> e = ROL::makePtr<OptStdVector>( e_ptr, fd_ );
  (*e_ptr)[i]= 1.0;
  return e;
}

int dimension() const {return static_cast<int>(std_vec_->size());}


const Vector<Real> & dual() const {
  ROL::Ptr<vector> dual_vecp = ROL::makePtr<vector>(*std_vec_);
  dual_vec_ = ROL::makePtr<OptDualStdVector<Real>>( dual_vecp, fd_ );
  fd_->apply(dual_vecp); 
  return *dual_vec_;
}

}; // class OptStdVector


// Dual optimization space.
template <class Real, class Element>
class OptDualStdVector : public Vector<Real> {

  typedef std::vector<Element>       vector;
  typedef typename vector::size_type uint;

private:
ROL::Ptr<std::vector<Element> >  std_vec_;
mutable ROL::Ptr<OptStdVector<Real> >  dual_vec_;
ROL::Ptr<FiniteDifference<Real> > fd_;

public:

OptDualStdVector(const ROL::Ptr<std::vector<Element> > & std_vec, ROL::Ptr<FiniteDifference<Real> >fd) : 
    std_vec_(std_vec), dual_vec_(ROL::nullPtr), fd_(fd) {}

void plus( const Vector<Real> &x ) {
  const OptDualStdVector &ex = dynamic_cast<const OptDualStdVector&>(x);
  ROL::Ptr<const vector> xvalptr = ex.getVector();
  uint dimension  = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] += (*xvalptr)[i];
  }
}

void scale( const Real alpha ) {
  uint dimension = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] *= alpha;
  }
}

Real dot( const Vector<Real> &x ) const {
  Real val = 0;
  const OptDualStdVector<Real, Element> & ex = dynamic_cast<const OptDualStdVector<Real, Element>&>(x);
  ROL::Ptr<const vector> kxvalptr = ex.getVector();
  ROL::Ptr<vector> xvalptr = ROL::makePtr<vector>(std_vec_->size(), 0.0);
  fd_->solve(kxvalptr,xvalptr);
  
  uint dimension  = std_vec_->size();
  for (unsigned i=0; i<dimension; i++) {
    val += (*std_vec_)[i]*(*xvalptr)[i];
  }
  return val;
}

Real norm() const {
  Real val = 0;
  val = std::sqrt( dot(*this) );
  return val;
}

ROL::Ptr<Vector<Real> > clone() const {
  return ROL::makePtr<OptDualStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()), fd_ );
}

ROL::Ptr<const std::vector<Element> > getVector() const {
  return std_vec_;
}

ROL::Ptr<std::vector<Element> > getVector()  {
  return std_vec_;
}

ROL::Ptr<Vector<Real> > basis( const int i ) const {
  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(), 0.0 );
  ROL::Ptr<OptDualStdVector> e = ROL::makePtr<OptDualStdVector>( e_ptr,fd_ );
  (*e_ptr)[i] = 1.0;
  return e;
}

int dimension() const {return static_cast<int>(std_vec_->size());}

const Vector<Real> & dual() const {
    ROL::Ptr<vector> dual_vecp = ROL::makePtr<vector>(*std_vec_); 
    dual_vec_ = ROL::makePtr<OptStdVector<Real>>( dual_vecp, fd_ );
    
    fd_->solve(dual_vecp);
    return *dual_vec_;
}

}; // class OptDualStdVector




// Constraint space.
template <class Real, class Element>
class ConStdVector : public Vector<Real> {

  typedef std::vector<Element> vector;
  typedef typename vector::size_type uint;

private:
ROL::Ptr<std::vector<Element> >  std_vec_;
mutable ROL::Ptr<ConDualStdVector<Real> >  dual_vec_;

public:

ConStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}

void plus( const Vector<Real> &x ) {
  const ConStdVector &ex = dynamic_cast<const ConStdVector&>(x);
  ROL::Ptr<const vector> xvalptr = ex.getVector();
  uint dimension  = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] += (*xvalptr)[i];
  }
}

void scale( const Real alpha ) {
  uint dimension = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] *= alpha;
  }
}

Real dot( const Vector<Real> &x ) const {
  Real val = 0;
  const ConStdVector<Real, Element> & ex = dynamic_cast<const ConStdVector<Real, Element>&>(x);
  ROL::Ptr<const vector> xvalptr = ex.getVector();

  uint dimension  = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    val += (*std_vec_)[i]*(*xvalptr)[i];
  }
  return val;
}

Real norm() const {
  Real val = 0;
  val = std::sqrt( dot(*this) );
  return val;
}

ROL::Ptr<Vector<Real> > clone() const {
  return ROL::makePtr<ConStdVector>( ROL::makePtr<vector>(std_vec_->size()));
}

ROL::Ptr<const std::vector<Element> > getVector() const {
  return std_vec_;
}

ROL::Ptr<std::vector<Element> > getVector() {
  return std_vec_;
}

ROL::Ptr<Vector<Real> > basis( const int i ) const {
  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
  ROL::Ptr<ConStdVector> e = ROL::makePtr<ConStdVector>( e_ptr);
  (*e_ptr)[i] = 1.0;
  return e;
}

int dimension() const {return static_cast<int>(std_vec_->size());}

const Vector<Real> & dual() const {
  dual_vec_ = ROL::makePtr<ConDualStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
  return *dual_vec_;
}

}; // class ConStdVector


// Dual constraint space.
template <class Real, class Element>
class ConDualStdVector : public Vector<Real> {

  typedef std::vector<Element>       vector;
  typedef typename vector::size_type uint;

private:

ROL::Ptr<std::vector<Element> >        std_vec_;
mutable ROL::Ptr<ConStdVector<Real> >  dual_vec_;

public:

ConDualStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}

void plus( const Vector<Real> &x ) {
  const ConDualStdVector &ex = dynamic_cast<const ConDualStdVector&>(x);
  ROL::Ptr<const vector> xvalptr = ex.getVector();
  uint dimension  = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] += (*xvalptr)[i];
  }
}

void scale( const Real alpha ) {
  uint dimension = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    (*std_vec_)[i] *= alpha;
  }
}

Real dot( const Vector<Real> &x ) const {
  Real val = 0;
  const ConDualStdVector<Real, Element> & ex = dynamic_cast<const ConDualStdVector<Real, Element>&>(x);
  ROL::Ptr<const vector> xvalptr = ex.getVector();
  uint dimension  = std_vec_->size();
  for (uint i=0; i<dimension; i++) {
    val += (*std_vec_)[i]*(*xvalptr)[i];
  }
  return val;
}

Real norm() const {
  Real val = 0;
  val = std::sqrt( dot(*this) );
  return val;
}

ROL::Ptr<Vector<Real> > clone() const {
  return ROL::makePtr<ConDualStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()));
}

ROL::Ptr<const std::vector<Element> > getVector() const {
  return std_vec_;
}

ROL::Ptr<std::vector<Element> > getVector() {
  return std_vec_;
}

ROL::Ptr<Vector<Real> > basis( const int i ) const {
  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
  ROL::Ptr<ConDualStdVector> e = ROL::makePtr<ConDualStdVector>( e_ptr );
  (*e_ptr)[i] = 1.0;
  return e;
}

int dimension() const {return static_cast<int>(std_vec_->size());}

const Vector<Real> &  dual() const {
  dual_vec_ = ROL::makePtr<ConStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
  return *dual_vec_;
}

}; // class ConDualStdVector

/*** End of declaration of four vector spaces. ***/



template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real> >
class Objective_GrossPitaevskii : public Objective<Real> {

    typedef std::vector<Real> vector;
    typedef typename vector::size_type uint;

    private:

        Real g_;    

        uint  nx_;     

        Real dx_;     
        
        ROL::Ptr<const std::vector<Real> > Vp_;    

        ROL::Ptr<FiniteDifference<Real> > fd_;


        void applyK(const Vector<Real> &v, Vector<Real> &kv) {

              

            // Pointer to direction vector 
            ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();

            // Pointer to action of Hessian on direction vector 
            ROL::Ptr<vector> kvp = dynamic_cast<XDual&>(kv).getVector();

            Real dx2 = dx_*dx_;

            (*kvp)[0] = (2.0*(*vp)[0]-(*vp)[1])/dx2;
  
            for(uint i=1;i<nx_-1;++i) {
                (*kvp)[i] = (2.0*(*vp)[i]-(*vp)[i-1]-(*vp)[i+1])/dx2;
            } 

            (*kvp)[nx_-1] = (2.0*(*vp)[nx_-1]-(*vp)[nx_-2])/dx2;

        } 

    public: 

        Objective_GrossPitaevskii(const Real &g, const Vector<Real> &V, ROL::Ptr<FiniteDifference<Real> > fd) : g_(g),  
            Vp_((dynamic_cast<const StdVector<Real>&>(V)).getVector()), fd_(fd)  {

            nx_ = Vp_->size(); 
            dx_ = (1.0/(1.0+nx_));
        }
           

    Real value( const Vector<Real> &psi, Real &tol ) {

            

        // Pointer to opt vector 
        ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();

        // Pointer to K applied to opt vector 
        ROL::Ptr<vector> kpsip = ROL::makePtr<vector>(nx_, 0.0);
        XDual kpsi(kpsip,fd_);

        Real J = 0;

        applyK(psi,kpsi);

        for(uint i=0;i<nx_;++i) {
            J += (*psip)[i]*(*kpsip)[i] + (*Vp_)[i]*pow((*psip)[i],2) + g_*pow((*psip)[i],4);
        } 
      
        // Rescale for numerical integration by trapezoidal rule
        J *= 0.5*dx_;

        return J;
    }


    void gradient( Vector<Real> &g, const Vector<Real> &psi, Real &tol ) {

            

        // Pointer to opt vector 
        ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();

        // Pointer to gradient vector 
        ROL::Ptr<vector> gp = dynamic_cast<XDual&>(g).getVector();

        // Pointer to K applied to opt vector 
        ROL::Ptr<vector> kpsip = ROL::makePtr<vector>(nx_, 0.0);
        XDual kpsi(kpsip,fd_);

        applyK(psi,kpsi);

        for(uint i=0;i<nx_;++i) {
            (*gp)[i] = ((*kpsip)[i] + (*Vp_)[i]*(*psip)[i] + 2.0*g_*pow((*psip)[i],3))*dx_;
        } 
      
    }




    void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol ) {

          

        // Pointer to opt vector 
        ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();

        // Pointer to direction vector 
        ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();

        // Pointer to action of Hessian on direction vector 
        ROL::Ptr<vector> hvp = dynamic_cast<XDual&>(hv).getVector();

        applyK(v,hv);
 
        for(uint i=0;i<nx_;++i) {
            (*hvp)[i] *= dx_;
            (*hvp)[i] += ( (*Vp_)[i] + 6.0*g_*pow((*psip)[i],2) )*(*vp)[i]*dx_;
        } 

   }

};


template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real>, class CPrim=StdVector<Real>, class CDual=StdVector<Real> >
class Normalization_Constraint : public Constraint<Real> {

    typedef std::vector<Real> vector;
    typedef typename vector::size_type uint;

    private:     
    uint nx_;
    Real dx_;
    ROL::Ptr<FiniteDifference<Real> > fd_;
    bool exactsolve_; 

    public:
    Normalization_Constraint(int n, Real dx, ROL::Ptr<FiniteDifference<Real> > fd, bool exactsolve) : 
        nx_(n), dx_(dx), fd_(fd), exactsolve_(exactsolve) {}          


    void value(Vector<Real> &c, const Vector<Real> &psi, Real &tol){

          

        // Pointer to constraint vector (only one element)
        ROL::Ptr<vector> cp = dynamic_cast<CPrim&>(c).getVector();

        // Pointer to opt vector 
        ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();

        (*cp)[0] = -1.0;
        for(uint i=0;i<nx_;++i) {
            (*cp)[0] += dx_*pow((*psip)[i],2);
        } 
    }


    void applyJacobian(Vector<Real> &jv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol){

          

        // Pointer to action of Jacobian of constraint on direction vector (yields scalar)
        ROL::Ptr<vector> jvp = dynamic_cast<CPrim&>(jv).getVector();

        // Pointer to direction vector     
        ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();

        // Pointer to opt vector 
        ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();
     
        (*jvp)[0] = 0;
        for(uint i=0;i<nx_;++i) {
            (*jvp)[0] += 2.0*dx_*(*psip)[i]*(*vp)[i];
        }
    }


    void applyAdjointJacobian(Vector<Real> &ajv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol){

          

        // Pointer to action of adjoint of Jacobian of constraint on direction vector (yields vector)
        ROL::Ptr<vector> ajvp = dynamic_cast<XDual&>(ajv).getVector();

        // Pointer to direction vector     
        ROL::Ptr<const vector> vp = dynamic_cast<const CDual&>(v).getVector();
 
        // Pointer to opt vector 
        ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();

        for(uint i=0;i<nx_;++i) {
            (*ajvp)[i] = 2.0*dx_*(*psip)[i]*(*vp)[0];
        }
    }


    void applyAdjointHessian(Vector<Real> &ahuv, const Vector<Real> &u, const Vector<Real> &v, 
                             const Vector<Real> &psi, Real &tol){


          

        // The pointer to action of constraint Hessian in u,v inner product
        ROL::Ptr<vector> ahuvp = dynamic_cast<XDual&>(ahuv).getVector();

        // Pointer to direction vector u     
        ROL::Ptr<const vector> up = dynamic_cast<const CDual&>(u).getVector();

        // Pointer to direction vector v     
        ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();

        // Pointer to opt vector 
        ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();
  
        for(uint i=0;i<nx_;++i) {
            (*ahuvp)[i] = 2.0*dx_*(*vp)[i]*(*up)[0];        
        }  
    }
     
    std::vector<Real> solveAugmentedSystem(Vector<Real> &v1, Vector<Real> &v2, const Vector<Real> &b1, 
                                           const Vector<Real> &b2, const Vector<Real> &psi, Real &tol) {

            

        if(exactsolve_) {
            ROL::Ptr<vector> v1p = dynamic_cast<XPrim&>(v1).getVector();    
            ROL::Ptr<vector> v2p = dynamic_cast<CDual&>(v2).getVector();
            ROL::Ptr<const vector> b1p = dynamic_cast<const XDual&>(b1).getVector();
            ROL::Ptr<const vector> b2p = dynamic_cast<const CPrim&>(b2).getVector();
            ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();
        
            ROL::Ptr<vector> jacp = ROL::makePtr<vector>(nx_, 0.0);
            ROL::Ptr<vector> b1dp = ROL::makePtr<vector>(nx_, 0.0);

            for(uint i=0;i<nx_;++i) {
                (*jacp)[i] = (*psip)[i];
                (*b1dp)[i] = (*b1p)[i];
            }
         
            // The Jacobian of the constraint is \f$c'(\psi)=2dx\psi\f$
            XDual jac(jacp,fd_);
            jac.scale(2.0*dx_);

            // A Dual-in-name-only version of b1, so we can compute the desired inner products involving inv(K) 
            XDual b1d(b1dp,fd_);
        
            // \f$ (c'K^{-1}*c'^\ast)^{-1} \f$ 
            Real d = 1.0/jac.dot(jac);
            Real p = jac.dot(b1d);

            (*v2p)[0] = d*(p-(*b2p)[0]);
     
            v1.set(jac.dual());
            v1.scale(-(*v2p)[0]);
            v1.plus(b1d.dual());  

            return std::vector<Real>(0);
        }     
        else{
            return Constraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,psi,tol);
        }
    }
};