ROL_DiodeCircuit.hpp

Go to the documentation of this file.

#ifndef ROL_DIODECIRCUIT_HPP
#define ROL_DIODECIRCUIT_HPP

#include "ROL_Objective.hpp"
#include "ROL_StdVector.hpp"
#include "ROL_ScaledStdVector.hpp"
#include "ROL_BoundConstraint.hpp"

#include <iostream>
#include <fstream>
#include <string>

namespace ROL {
namespace ZOO {

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

    typedef std::vector<Real>            vector;
    typedef Vector<Real>                 V;
    typedef StdVector<Real>              STDV;
    typedef PrimalScaledStdVector<Real>  PSV;
    typedef DualScaledStdVector<Real>    DSV; 
    typedef typename vector::size_type   uint;

  private:
    Real Vth_; 
    Teuchos::RCP<std::vector<Real> > Imeas_;
    Teuchos::RCP<std::vector<Real> > Vsrc_; 
    bool lambertw_; 
    Real noise_; 
    bool use_adjoint_;
    int use_hessvec_;
  
    Teuchos::RCP<const vector> getVector( const V& x ) {
      using Teuchos::dyn_cast;  using Teuchos::getConst;
      try { 
        return dyn_cast<const STDV>(getConst(x)).getVector();
      }
      catch (std::exception &e) {
        try { 
          return dyn_cast<const PSV>(getConst(x)).getVector();
        }
        catch (std::exception &e) {
          return dyn_cast<const DSV>(getConst(x)).getVector();
        }
      }
    }

    Teuchos::RCP<vector> getVector( V& x ) {
      using Teuchos::dyn_cast;
      try {
        return dyn_cast<STDV>(x).getVector(); 
      }
      catch (std::exception &e) {
        try {
          return dyn_cast<PSV>(x).getVector(); 
        }
        catch (std::exception &e) {
          return dyn_cast<DSV>(x).getVector(); 
        }
      }
    }

  public:

    Objective_DiodeCircuit(Real Vth, Real Vsrc_min, Real Vsrc_max, Real Vsrc_step, Real true_Is, Real true_Rs, bool lambertw, Real noise, bool use_adjoint, int use_hessvec){
      lambertw_ = lambertw; 
      Vth_ = Vth;      
      use_adjoint_ = use_adjoint;
      use_hessvec_ = use_hessvec;
      int n = (Vsrc_max-Vsrc_min)/Vsrc_step + 1;
      Vsrc_ = Teuchos::rcp(new std::vector<Real>(n,0.0));
      Imeas_ = Teuchos::rcp(new std::vector<Real>(n,0.0));
      std::ofstream output ("Measurements.dat");
      Real left = 0.0; Real right = 1.0;
      if(lambertw_){
        // Using Lambert-W function
        std::cout << "Generating data using Lambert-W function." << std::endl;
        for(int i=0;i<n;i++){
          (*Vsrc_)[i] = Vsrc_min+i*Vsrc_step;
          (*Imeas_)[i] = lambertWCurrent(true_Is,true_Rs,(*Vsrc_)[i]);
          if(noise>0.0){
            (*Imeas_)[i] += noise*pow(10,int(log10((*Imeas_)[i])))*(( (Real)rand() / (Real)RAND_MAX ) * (right - left) + left);
          }
          // Write generated data into file
          if(output.is_open()){
            output << std::setprecision(8) << std::scientific << (*Vsrc_)[i] << "  " << (*Imeas_)[i] << "\n";
          }
        }       
      }
      else{
        // Using Newton's method
        std::cout << "Generating data using Newton's method." << std::endl;
        for(int i=0;i<n;i++){
          (*Vsrc_)[i] = Vsrc_min+i*Vsrc_step;
          Real I0 = 1.e-12; // initial guess for Newton
          (*Imeas_)[i] = Newton(I0,Vsrc_min+i*Vsrc_step,true_Is,true_Rs);       
          if(noise>0.0){
            (*Imeas_)[i] += noise*pow(10,int(log10((*Imeas_)[i])))*(( (Real)rand() / (Real)RAND_MAX ) * (right - left) + left);
          }
          // Write generated data into file
          if(output.is_open()){
            output << std::setprecision(8) << std::scientific << (*Vsrc_)[i] << "  " << (*Imeas_)[i] << "\n";
          }
        }
      }
        
      output.close();

    }

    Objective_DiodeCircuit(Real Vth, std::ifstream& input_file, bool lambertw, Real noise, bool use_adjoint, int use_hessvec){
      lambertw_ = lambertw; 
      Vth_ = Vth;     
      use_adjoint_ = use_adjoint;
      use_hessvec_ = use_hessvec;
      std::string line;
      int dim = 0;
      for(int k=0;std::getline(input_file,line);++k){dim=k;} // count number of lines
      input_file.clear(); // reset to beginning of file
      input_file.seekg(0,std::ios::beg); 
      Vsrc_ = Teuchos::rcp(new std::vector<Real>(dim,0.0));
      Imeas_ = Teuchos::rcp(new std::vector<Real>(dim,0.0));
      double Vsrc, Imeas;
      std::cout << "Using input file to generate data." << "\n";
      for(int i=0;i<dim;i++){
        input_file >> Vsrc;
        input_file >> Imeas;
        (*Vsrc_)[i] = Vsrc;
        (*Imeas_)[i] = Imeas;
      }
      input_file.close();

    }

    void set_method(bool lambertw){
      lambertw_ = lambertw;
    }

    Real diode(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return I-Is*(exp((Vsrc-I*Rs)/Vth_)-1);
    }
    
    Real diodeI(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return 1+Is*exp((Vsrc-I*Rs)/Vth_)*(Rs/Vth_);
    }

    Real diodeIs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return 1-exp((Vsrc-I*Rs)/Vth_);
    }

    Real diodeRs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return Is*exp((Vsrc-I*Rs)/Vth_)*(I/Vth_);
    }
    
    Real diodeII(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return -Is*exp((Vsrc-I*Rs)/Vth_)*(Rs/Vth_)*(Rs/Vth_);
    }

    Real diodeIIs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return exp((Vsrc-I*Rs)/Vth_)*(Rs/Vth_);
    }

    Real diodeIRs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return (Is/Vth_)*exp((Vsrc-I*Rs)/Vth_)*(1-(I*Rs)/Vth_);
    }

    Real diodeIsIs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return 0;
    }
    
    Real diodeIsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return exp((Vsrc-I*Rs)/Vth_)*(I/Vth_);
    }

    Real diodeRsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      return -Is*exp((Vsrc-I*Rs)/Vth_)*(I/Vth_)*(I/Vth_);
    }

    Real Newton(const Real I, const Real Vsrc, const Real Is, const Real Rs){
      double EPS = 1.e-16;
      double TOL = 1.e-13;
      int MAXIT = 200;
      Real IN = I;
      Real fval  = diode(IN,Vsrc,Is,Rs);
      Real dfval = 0.0;
      Real IN_tmp   = 0.0;
      Real fval_tmp = 0.0;
      Real alpha = 1.0;
      for ( int i = 0; i < MAXIT; i++ ) {
        if ( std::fabs(fval) < TOL ) {
          // std::cout << "converged with |fval| = " << std::fabs(fval) << " and TOL = " << TOL << "\n";
          break;
        }
        dfval = diodeI(IN,Vsrc,Is,Rs);
        if( std::fabs(dfval) < EPS ){
          std::cout << "denominator is too small" << std::endl;
          break;
        }
        
        alpha    = 1.0;
        IN_tmp   = IN - alpha*fval/dfval;
        fval_tmp = diode(IN_tmp,Vsrc,Is,Rs);
        while ( std::fabs(fval_tmp) >= (1.0-1.e-4*alpha)*std::fabs(fval) ) {
          alpha   /= 2.0;
          IN_tmp   = IN - alpha*fval/dfval;
          fval_tmp = diode(IN_tmp,Vsrc,Is,Rs);
          if ( alpha < std::sqrt(EPS) ) { 
            // std::cout << "Step tolerance met\n";
            break;
          }
        }
        IN   = IN_tmp;
        fval = fval_tmp;
        // if ( i == MAXIT-1){
        //   std::cout << "did not converge  " << std::fabs(fval) << "\n";
        // }
      }
      return IN;
    }    
    
    void lambertw(Real x, Real &w, int &ierr, Real &xi){
      int i=0, maxit = 10;
      const Real turnpt = -exp(-1.), c1 = 1.5, c2 = .75;
      Real r, r2, r3, s, mach_eps, relerr = 1., diff;
      mach_eps = 2.e-15;   // float:2e-7
      ierr = 0;
      
      if( x > c1){
        w = c2*log(x);
        xi = log( x/ w) - w;
      }
      else{
        if( x >= 0.0){
          w = x;
          if( x == 0. ) return;
          if( x < (1-c2) ) w = x*(1.-x + c1*x*x);
          xi = - w;
        }
        else{
          if( x >= turnpt){
                if( x > -0.2 ){
                  w = x*(1.0-x + c1*x*x);
                  xi = log(1.0-x + c1*x*x) - w;
                }
                else{
                  diff = x-turnpt;
                  if( diff < 0.0 ) diff = -diff;
                  w = -1 + sqrt(2.0*exp(1.))*sqrt(x-turnpt);
                  if( diff == 0.0 ) return;
                  xi = log( x/ w) - w;
                }
          }
          else{
            ierr = 1; // x is not in the domain.
            w = -1.0;
            return;
          }
        }
      }
      
      while( relerr > mach_eps  && i<maxit){
        r = xi/(w+1.0);   //singularity at w=-1
        r2 = r*r;
        r3 = r2*r;
        s  = 6.*(w+1.0)*(w+1.0);
        w = w * (  1.0 + r + r2/(2.0*( w+1.0)) - (2. * w -1.0)*r3/s  );
        if( w * x < 0.0 ) w = -w;
        xi = log( x/ w) - w;
        
        if( x>1.0 ){
          relerr =  xi / w;
        }
        else{
          relerr =  xi;
        }
        if(relerr < 0.0 ) relerr = -relerr;
        ++i;
      }
      if( i == maxit ) ierr = 2;
    }
    
    Real lambertWCurrent(Real Is, Real Rs, Real Vsrc){
      Real arg1 = (Vsrc + Is*Rs)/Vth_;
      Real evd = exp(arg1);
      Real lambWArg = Is*Rs*evd/Vth_;
      Real lambWReturn;
      int ierr;
      Real lambWError;
      lambertw(lambWArg, lambWReturn, ierr, lambWError);
      if(ierr == 1){std::cout << "LambertW error: argument is not in the domain" <<  std::endl; return -1.0;}
      if(ierr == 2){std::cout << "LambertW error: BUG!" <<  std::endl;}
      Real Id = -Is+Vth_*(lambWReturn)/Rs;
      //Real Gd = lambWReturn / ((1 + lambWReturn)*RS);     
      return Id;
    }


    void solve_circuit(Vector<Real> &I, const Vector<Real> &S){
      using Teuchos::RCP;
      RCP<vector> Ip = getVector(I);
      RCP<const vector> Sp = getVector(S);

      uint n = Ip->size();
      
      if(lambertw_){
        // Using Lambert-W function
        Real lambval;
        for(uint i=0;i<n;i++){
          lambval = lambertWCurrent((*Sp)[0],(*Sp)[1],(*Vsrc_)[i]);
          (*Ip)[i] = lambval;
        }
      }
      else{     
        // Using Newton's method      
        Real I0 = 1.e-12; // Initial guess for Newton
        for(uint i=0;i<n;i++){
          (*Ip)[i] = Newton(I0,(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
        }
      }
    }

    
    Real value(const Vector<Real> &S, Real &tol){
      using Teuchos::RCP;  using Teuchos::rcp;
      RCP<const vector> Sp = getVector(S);
      uint n = Imeas_->size();
      STDV I( rcp( new vector(n,0.0) ) );
      RCP<vector> Ip = getVector(I);

      // Solve state equation
      solve_circuit(I,S);
      Real val = 0;
      
      for(uint i=0;i<n;i++){
        val += ((*Ip)[i]-(*Imeas_)[i])*((*Ip)[i]-(*Imeas_)[i]);
      }
      return val/2.0;
    }
    
    void solve_adjoint(Vector<Real> &lambda, const Vector<Real> &I, const Vector<Real> &S){
      
      using Teuchos::RCP;
      RCP<vector> lambdap = getVector(lambda);
      RCP<const vector> Ip = getVector(I);
      RCP<const vector> Sp = getVector(S);
      
      uint n = Ip->size();
      for(uint i=0;i<n;i++){
        (*lambdap)[i] = ((*Imeas_)[i]-(*Ip)[i])/diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
      }
    }

    void solve_sensitivity_Is(Vector<Real> &sens, const Vector<Real> &I, const Vector<Real> &S){
 
      using Teuchos::RCP;
      RCP<vector> sensp = getVector(sens);
      RCP<const vector> Ip = getVector(I);
      RCP<const vector> Sp = getVector(S);     
      
      uint n = Ip->size();
      for(uint i=0;i<n;i++){
        (*sensp)[i] = -diodeIs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])/diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
      }
    }
    
    void solve_sensitivity_Rs(Vector<Real> &sens, const Vector<Real> &I, const Vector<Real> &S){
           
      using Teuchos::RCP;
      RCP<vector> sensp = getVector(sens);
      RCP<const vector> Ip = getVector(I);
      RCP<const vector> Sp = getVector(S);
      
      uint n = Ip->size();
      for(uint i=0;i<n;i++){
        (*sensp)[i] = -diodeRs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])/diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
      }
    }
    
    void gradient(Vector<Real> &g, const Vector<Real> &S, Real &tol){

      using Teuchos::RCP;  using Teuchos::rcp;
      RCP<vector> gp = getVector(g);
      RCP<const vector> Sp = getVector(S);
      
      uint n = Imeas_->size();
      
      STDV I( rcp( new vector(n,0.0) ) );
      RCP<vector> Ip = getVector(I);
      
      // Solve state equation      
      solve_circuit(I,S);
      
      if(use_adjoint_){      
        // Compute the gradient of the reduced objective function using adjoint computation
        STDV lambda( rcp( new vector(n,0.0) ) );
        RCP<vector> lambdap = getVector(lambda);
        
        // Solve adjoint equation
        solve_adjoint(lambda,I,S);
      
        // Compute gradient
        (*gp)[0] = 0.0; (*gp)[1] = 0.0;
        for(uint i=0;i<n;i++){
          (*gp)[0] += diodeIs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])*(*lambdap)[i];
          (*gp)[1] += diodeRs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])*(*lambdap)[i];            
        }      
      }
      else{
        // Compute the gradient of the reduced objective function using sensitivities
        STDV sensIs( rcp( new vector(n,0.0) ) );
        STDV sensRs( rcp( new vector(n,0.0) ) );
        // Solve sensitivity equations
        solve_sensitivity_Is(sensIs,I,S);
        solve_sensitivity_Rs(sensRs,I,S);
        
        RCP<vector> sensIsp = getVector(sensIs);
        RCP<vector> sensRsp = getVector(sensRs);

        // Write sensitivities into file
        std::ofstream output ("Sensitivities.dat");
        for(uint k=0;k<n;k++){
          if(output.is_open()){
            output << std::scientific << (*sensIsp)[k] << " " << (*sensRsp)[k] << "\n";
          }     
        }
        output.close();
        // Compute gradient
        (*gp)[0] = 0.0; (*gp)[1] = 0.0;
        for(uint i=0;i<n;i++){
          (*gp)[0] += ((*Ip)[i]-(*Imeas_)[i])*(*sensIsp)[i];
          (*gp)[1] += ((*Ip)[i]-(*Imeas_)[i])*(*sensRsp)[i];    
        }      
      }
    }
    

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

      using Teuchos::RCP;  using Teuchos::rcp;    

      if(use_hessvec_==0){
        Objective<Real>::hessVec(hv, v, S, tol);
      }
      else if(use_hessvec_==1){
        RCP<vector> hvp = getVector(hv);
        RCP<const vector> vp = getVector(v);
        RCP<const vector> Sp = getVector(S);
        
        uint n = Imeas_->size();
        
        STDV I( rcp( new vector(n,0.0) ) );
        RCP<vector> Ip = getVector(I);
        
        // Solve state equation      
        solve_circuit(I,S);
        
        STDV lambda( rcp( new vector(n,0.0) ) );
        RCP<vector> lambdap = getVector(lambda);
        
        // Solve adjoint equation
        solve_adjoint(lambda,I,S);
        
        STDV w( rcp( new vector(n,0.0) ) );
        RCP<vector> wp = getVector(w);
        
        // Solve state sensitivity equation
        for(uint i=0;i<n;i++){
          (*wp)[i] = ( (*vp)[0] * diodeIs( (*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1] ) + (*vp)[1] * diodeRs( (*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1] ) ) / diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
        }
        
        STDV p( rcp( new vector(n,0.0) ) );
        RCP<vector> pp = getVector(p);
        
        // Solve for p
        for(uint j=0;j<n;j++){
          (*pp)[j] = ( (*wp)[j] + (*lambdap)[j] * diodeII( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] ) * (*wp)[j] - (*lambdap)[j] * diodeIIs( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] ) * (*vp)[0] - (*lambdap)[j] * diodeIRs( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] ) * (*vp)[1] ) / diodeI( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] );
        }
        
        // Assemble Hessian-vector product
        (*hvp)[0] = 0.0;(*hvp)[1] = 0.0;
        for(uint k=0;k<n;k++){
          (*hvp)[0] += diodeIs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) * (*pp)[k] - (*lambdap)[k] * (*wp)[k] * diodeIIs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) + (*lambdap)[k] * (*vp)[0] * diodeIsIs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) + (*lambdap)[k] * (*vp)[1] * diodeIsRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] );
          (*hvp)[1] += diodeRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) * (*pp)[k] - (*lambdap)[k] * (*wp)[k] * diodeIRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) + (*lambdap)[k] * (*vp)[0] * diodeIsRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) + (*lambdap)[k] * (*vp)[1] * diodeRsRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] );
        }
      }
      else if(use_hessvec_==2){
        //Gauss-Newton approximation
        RCP<vector> hvp = getVector(hv);
        RCP<const vector> vp = getVector(v);
        RCP<const vector> Sp = getVector(S);
        
        uint n = Imeas_->size();

        STDV I( rcp( new vector(n,0.0) ) );
        RCP<vector> Ip = getVector(I);

        // Solve state equation                                                                                
        solve_circuit(I,S);

        // Compute sensitivities
        STDV sensIs( rcp( new vector(n,0.0) ) );
        STDV sensRs( rcp( new vector(n,0.0) ) );

        // Solve sensitivity equations                                                                          
        solve_sensitivity_Is(sensIs,I,S);
        solve_sensitivity_Rs(sensRs,I,S);
        RCP<vector> sensIsp = getVector(sensIs);
        RCP<vector> sensRsp = getVector(sensRs);
        
        // Compute approximate Hessian
        Real H11 = 0.0; Real H12 = 0.0; Real H22 = 0.0;
        for(uint k=0;k<n;k++){
          H11 += (*sensIsp)[k]*(*sensIsp)[k];
          H12 += (*sensIsp)[k]*(*sensRsp)[k];
          H22 += (*sensRsp)[k]*(*sensRsp)[k];
        }
        
        // Compute approximate Hessian-times-vector
        (*hvp)[0] = H11*(*vp)[0] + H12*(*vp)[1];
        (*hvp)[1] = H12*(*vp)[0] + H22*(*vp)[1];
      }
      else{
        ROL::Objective<Real>::hessVec( hv, v, S, tol ); // Use parent class function    
      }
    }

    void generate_plot(Real Is_lo, Real Is_up, Real Is_step, Real Rs_lo, Real Rs_up, Real Rs_step){
      Teuchos::RCP<std::vector<double> > S_rcp = Teuchos::rcp(new std::vector<double>(2,0.0) );
      StdVector<double> S(S_rcp);
      std::ofstream output ("Objective.dat");

      Real Is = 0.0;
      Real Rs = 0.0;
      Real val = 0.0;
      Real tol = 1.e-16;
      int n = (Is_up-Is_lo)/Is_step + 1;
      int m = (Rs_up-Rs_lo)/Rs_step + 1;
      for(int i=0;i<n;i++){
        Is = Is_lo + i*Is_step;
        for(int j=0;j<m;j++){
          Rs = Rs_lo + j*Rs_step;
          (*S_rcp)[0] = Is;
          (*S_rcp)[1] = Rs;
          val = value(S,tol);
          if(output.is_open()){
            output << std::scientific << Is << " " << Rs << " " << val << "\n";
          }
        }
      }
      output.close();
    }


  };  // class Objective_DiodeCircuit 


  // template<class Real>
  // void getDiodeCircuit( Teuchos::RCP<Objective<Real> > &obj, Vector<Real> &x0, Vector<Real> &x ) {
  //   // Cast Initial Guess and Solution Vectors                                     
  //   Teuchos::RCP<std::vector<Real> > x0p =
  //     Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalScaledStdVector<Real> >(x0)).getVector());
  //   Teuchos::RCP<std::vector<Real> > xp =
  //     Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalScaledStdVector<Real> >(x)).getVector());

  //   int n = xp->size();

  //   // Resize Vectors                                                                                              
  //   n = 2;
  //   x0p->resize(n);
  //   xp->resize(n);

  //   // Instantiate Objective Function                                                                              
  //   obj = Teuchos::rcp( new Objective_DiodeCircuit<Real> (0.02585,0.0,1.0,1.e-2));
  //   //ROL::Objective_DiodeCircuit<Real> obj(0.02585,0.0,1.0,1.e-2);

  //   // Get Initial Guess
  //   (*x0p)[0] = 1.e-13;
  //   (*x0p)[1] = 0.2;
    
  //   // Get Solution
  //   (*xp)[0] = 1.e-12;
  //   (*xp)[1] = 0.25;
    
  // }


} //end namespace ZOO
} //end namespace ROL

#endif