Problem using `pow` in a `Sacado::Fad::SFad` template

[Tim Molter]

In my memristor model, I'm including a window function like this:

template <typename ScalarT>

void X_var_F( const ScalarT & V1, const ScalarT & V2, const ScalarT & X, const ScalarT & UV, const ScalarT & RON, const ScalarT & ROFF, const ScalarT & D, const ScalarT & PCOEFF, ScalarT & fval ){

  fval = UV*RON/D/D/(RON * X  + ROFF * (1 - X))*(V1-V2)*(1.0 - pow( (2.0*X - 1.0), 2.0*PCOEFF.val() ));

}

,where the last major term is the window. If I don't use `PCOEFF.val()`, but instead `PCOEFF.val()`, then `fval` = 

fval = 0 [ nan nan nan ]

, which stops the sim in it's tracks. Otherwise, with .val(), I get:

fval = 0 [ 2.5e-11 -2.5e-11 0 ]

, which is correct. In the TEAM model, the window functions all use the non- .val(). Those TEAM memristor simulations are also not working for me, failing at the DCOP step, and I suspect this is why. 

Any idea what's going on? 

[xyce-users]

It is difficult to answer this question without more detail.  For example, what is the value of PCOEFF and what are its derivative values?

The Sacado "pow" function may not be well-behaved at certain values of its arguments, and this may be system-dependent.  (Once upon a time, it was dangerous to call with zero as the first argument, for example, though a quick check of the code shows that they now try to guard against that).  When you use PCOEF.val(), you are assuring that the expression will never have the right derivatives if PCOEFF depends on probes, as the exponent will be treated as a double constant with no derivative information, so it is not correct to do so.

The Sacado pow function is not only attempting to evaluate pow(f,g) but also the derivatives with respect to all independent variables.  Doing so involves taking the log of f and the inverse of f.  If f is negative, this could lead to NANs in the derivatives.  You will need to check the values of these arguments before calling the pow function, or structure your code so that these can never be zero (difficult to do because the code may sample non-physical values of variables during its attempt to converge).  Device models must be written to be well-behaved at all input values, not just physically valid input values.

[Tim Molter]

I'll try to show what I'm doing. First my template methods, then my usage of them in my `Master:updateState` method. It's more or less elements of the TEAM and Yakopcic model combined. Also `P` (formerly `PCOEFF`) is just a constant in the model. Here's what I see in the console for the first iteration:

  P = 2 [ 0 0 0 ]

  X = 0.001 [ 0 0 1 ]

  fval = 0.00797603 [ nan nan nan ]

  ri.xVarFContribution = 0

  ri.dxFEqdVpos = nan

  ri.dxFEqdVneg = nan

  ri.dxFEqdx = nan

template <typename ScalarT>

void JogelkarWindowFunction( const ScalarT & X, const ScalarT & P, ScalarT & fval )

{

fval = 1.0 - pow( (2.0*X - 1.0), (2*P) );

if (DEBUG_DEVICE){

Xyce::dout()  << "  P = " <<  P << std::endl;

Xyce::dout()  << "  X = " <<  X << std::endl;

Xyce::dout()  << "  fval = " <<  fval << std::endl;

}

}

// linear current voltage model, Reff

template <typename ScalarT>

void Reff( const ScalarT & X, const ScalarT & RON, const ScalarT & ROFF,  ScalarT & fval )

{

fval = RON * X  + ROFF * (1 - X);

}

template <typename ScalarT>

void X_var_F( const ScalarT & V1, const ScalarT & V2, const ScalarT & X, const ScalarT & UV, const ScalarT & RON, const ScalarT & ROFF, const ScalarT & D, const ScalarT & P, ScalarT & fval ){

    ScalarT Rval;

    Reff( X, RON, ROFF, Rval );

    ScalarT Windowval;

    JogelkarWindowFunction(X, P, Windowval );

  fval = UV*RON/D/D/Rval*(V1-V2)*Windowval;

}

Master:updateState

        {

          // evaluate the state variable equation

          //Sacado::Fad::SFad<double,2> varDeltaV( 2, 0, (v_pos-v_neg) );

          Sacado::Fad::SFad<double,3> varV1( 3, 0, v_pos );

          Sacado::Fad::SFad<double,3> varV2( 3, 1, v_neg );

          Sacado::Fad::SFad<double,3> varX( 3, 2, x );

          Sacado::Fad::SFad<double,3> paramUV( ri.model_.uv_ );

          Sacado::Fad::SFad<double,3> paramD( ri.model_.D_ );

          Sacado::Fad::SFad<double,3> paramRON( ri.model_.Ron_ );

          Sacado::Fad::SFad<double,3> paramROFF( ri.model_.Roff_ );

          Sacado::Fad::SFad<double,3> paramP( ri.model_.p_ );

          Sacado::Fad::SFad<double,3> resultFad;

          X_var_F( varV1, varV2, varX, paramUV, paramRON, paramROFF, paramD, paramP, resultFad );

          ri.xVarFContribution = resultFad.val();

          ri.dxFEqdVpos = resultFad.dx(0);

          ri.dxFEqdVneg = resultFad.dx(1);

          ri.dxFEqdx = resultFad.dx(2);

if (DEBUG_DEVICE){

Xyce::dout()  << "  ri.xVarFContribution = " <<  ri.xVarFContribution << std::endl;

Xyce::dout()  << "  ri.dxFEqdVpos = " <<  ri.dxFEqdVpos << std::endl;

Xyce::dout()  << "  ri.dxFEqdVneg = " <<  ri.dxFEqdVneg << std::endl;

Xyce::dout()  << "  ri.dxFEqdx = " <<  ri.dxFEqdx << std::endl;

}

        }

[Tim Molter]

Here's a better formatted version...

I'll try to show what I'm doing. First my template methods, then my usage of them in my `Master:updateState` method. It's more or less elements of the TEAM and Yakopcic model combined. Also `P` (formerly `PCOEFF`) is just a constant in the model. Here's what I see in the console for the first iteration:

Output

P = 2 [ 0 0 0 ]
X = 0.001 [ 0 0 1 ]
fval = 0.00797603 [ nan nan nan ]
ri.xVarFContribution = 0
ri.dxFEqdVpos = nan
ri.dxFEqdVneg = nan
ri.dxFEqdx = nan

Templates

template <typename ScalarT>
void JogelkarWindowFunction( const ScalarT & X, const ScalarT & P, ScalarT & fval )
{
  fval = 1.0 - pow( (2.0*X - 1.0), (2*P) );
  if (DEBUG_DEVICE){
    Xyce::dout()  << "  P = " <<  P << std::endl;
    Xyce::dout()  << "  X = " <<  X << std::endl;
    Xyce::dout()  << "  fval = " <<  fval << std::endl;
  }
}


// linear current voltage model, Reff
template <typename ScalarT>

void Reff( const ScalarT & X, const ScalarT & RON, const ScalarT & ROFF,  ScalarT & fval )
{
  fval = RON * X  + ROFF * (1 - X);
}




template <typename ScalarT>

void X_var_F( const ScalarT & V1, const ScalarT & V2, const ScalarT & X, const ScalarT & UV, const ScalarT & RON, const ScalarT & ROFF, const ScalarT & D, const ScalarT & P, ScalarT & fval ){


  ScalarT Rval;
  Reff( X, RON, ROFF, Rval );
  ScalarT Windowval;
  JogelkarWindowFunction(X, P, Windowval );


  fval = UV*RON/D/D/Rval*(V1-V2)*Windowval;
}

Master:updateState

{


  // evaluate the state variable equation

[xyce-users]

Ok, now I understand.

To solve this, your best bet is not to use Sacado FAD types for anything except variables that depend on the solution --- Sacado is generally smart enough to know how to work with a mixture of doubles and SFads, so it is not necessary to make all expressions and templated functions use only FAD types.  And in fact, if you call pow(a,b) where a is a FAD type and b is a double, it does something different than if a is FAD and b is FAD with all zero derivatives.

So, what I'm suggesting is basically that you not try to mimic what was done in N_DEV_MemristorTEAM.C, which for some reason unnecessarily copies constants into FAD variables.  I have not tested these out, but I suggest you replace your templated functions like this:

template <typename ScalarT>
ScalarT JogelkarWindowFunction( const ScalarT & X, double P)
{
  ScalarT fval;

  fval = 1.0 - pow( (2.0*X - 1.0), (2*P) );
  if (DEBUG_DEVICE){
    Xyce::dout()  << "  P = " <<  P << std::endl;
    Xyce::dout()  << "  X = " <<  X << std::endl;
    Xyce::dout()  << "  fval = " <<  fval << std::endl;
  }

  return fval;

}
 
 
// linear current voltage model, Reff
template <typename ScalarT>

ScalarT Reff( const ScalarT & X, double RON, double ROFF)
{
  return(RON * X  + ROFF * (1 - X));
}
 


 
template <typename ScalarT>
ScalarT X_var_F( const ScalarT & V1, const ScalarT & V2, const ScalarT & X,
double UV, double RON, double ROFF, double D, double P){
 
 
  ScalarT Rval=Reff( X, RON, ROFF );
  ScalarT Windowval=JogelkarWindowFunction(X, P);


  return(UV*RON/D/D/Rval*(V1-V2)*Windowval);
}

Then, in your updateState, instead of what you're doing, do this:

  // evaluate the state variable equation
  Sacado::Fad::SFad<double,3> varV1( 3, 0, v_pos );
  Sacado::Fad::SFad<double,3> varV2( 3, 1, v_neg );
  Sacado::Fad::SFad<double,3> varX( 3, 2, x );

  Sacado::Fad::SFad<double,3> resultFad;


  resultFad = X_var_F( varV1, varV2, varX, ri.model_.uv_, ri.model_.Ron_, ri.model_.Roff_, ri.model_.D_, ri.model_.p_);

Because the model parameters have no dependence on any solution variables, it is completely unnecessary (and inefficient) to maintain them as FAD types, which track derivative information.    All the derivatives are set to zero when you declare the FAD variable as you have been doing, and clearly *some* Sacado functions (like pow) are tricky that way.

This should solve all your NAN issues for the functions that were having trouble with pow. 

[Tim Molter]

Worked like a charm. Much appreciated!