gross-pitaevskii/example_01.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @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"


using namespace ROL;

template<class Real>
class Objective_GrossPitaevskii : public Objective<Real> {

  typedef std::vector<Real>  vector;
  typedef Vector<Real>       V;
  typedef StdVector<Real>    SV;

  typedef typename vector::size_type uint;

private:


  Real g_;    

  uint  nx_;     

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

  ROL::Ptr<const vector> getVector( const V& x ) {
    
    return dynamic_cast<const SV&>(x).getVector();
  }

  ROL::Ptr<vector> getVector( V& x ) {
    
    return dynamic_cast<SV&>(x).getVector();
  }




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

    using namespace Teuchos;

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

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

    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) : g_(g),  
    Vp_(getVector(V))  {
    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 = getVector(psi);

    // Pointer to K applied to opt vector 
    ROL::Ptr<V> kpsi = psi.clone();
    ROL::Ptr<vector> kpsip = getVector(*kpsi);

    Real J = 0;

    this->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 = getVector(psi);

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

    // Pointer to K applied to opt vector 
    ROL::Ptr<V> kpsi = psi.clone();
    ROL::Ptr<vector> kpsip = getVector(*kpsi);

    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 = getVector(psi);

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

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

    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 Normalization_Constraint : public Constraint<Real> {

  typedef std::vector<Real>  vector;
  typedef Vector<Real>       V;
  typedef StdVector<Real>    SV;

  typedef typename vector::size_type uint;


private:     
  uint nx_;
  Real dx_;

  ROL::Ptr<const vector> getVector( const V& x ) {
    
    return dynamic_cast<const SV&>(x).getVector();
  }

  ROL::Ptr<vector> getVector( V& x ) {
    
    return dynamic_cast<SV&>(x).getVector();
  }

public:
  Normalization_Constraint(int n, Real dx) : nx_(n), dx_(dx) {}          


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

    

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

    // Pointer to optimization vector     
    ROL::Ptr<const vector> psip = getVector(psi);

    (*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 = getVector(jv);

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

    // Pointer to optimization vector     
    ROL::Ptr<const vector> psip = getVector(psi);
      
        (*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 = getVector(ajv);

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

    // Pointer to optimization vector     
    ROL::Ptr<const vector> psip = getVector(psi);

    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 = getVector(ahuv);

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

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

    // Pointer to optimization vector     
    ROL::Ptr<const vector> psip = getVector(psi);

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