Derivative term in Reynolds Stress Equation A5 always evaluates to 0

[thclark]

In equation A5 of Perry and Marusic 1995 Part 1., the derivative dWc[1, Pi] / dPi always evaluates to 0.

This is because taking the derivative of equation 9 (for Wc) with respect to Pi gives a derivative in (1-eta), thus setting eta=1 forces the derivative to 0.

Equation A5 must therefore always simplify to N = Wc[1,Pi]. It is possible this slipped the net into the paper if that term is there as an appendage from a mathematica derivation, for example, but could represent an error which we're not accounting for.

So in stress.h, the derivative appears to be computed unnecessarily. Here is how to compute the derivative, using central differencing...

    // N (from eqn A5) using central differencing
    const double d_pi = 0.01 * pi_coles;
    double wc_minus = coles_wake(1.0, pi_coles - d_pi, true);
    double wc_plus =  coles_wake(1.0, pi_coles + d_pi, true);
    double n = coles_wake(1.0, pi_coles) + pi_coles * (wc_plus - wc_minus) / (2.0 * d_pi);
    std::cout << "Gradient (numeric diff) " << (wc_plus - wc_minus) / (2.0 * d_pi) << std::endl;

or Eigen's autodiff (recommended)

    // N (from eqn 5) using Eigen's AutoDiffScalar to determine dWc/dPi
    typedef Eigen::AutoDiffScalar<Eigen::VectorXd> ADScalar;
    ADScalar ads_one;
    ads_one.value() = 1.0;
    ads_one.derivatives() = Eigen::VectorXd::Unit(2,0);
    ADScalar ads_pi_coles;
    ads_pi_coles.value() = pi_coles;
    ads_pi_coles.derivatives() = Eigen::VectorXd::Unit(2,1);
    ADScalar ads_wc = coles_wake(ads_one, ads_pi_coles, false);
    double dwc_dpi_coles = ads_wc.derivatives()[0];
    std::cout << "Gradient (autodiff) " << dwc_dpi_coles << std::endl;
    double n = coles_wake(1.0, pi_coles) + pi_coles * dwc_dpi_coles;

In the meantime, we write this out - there's no point computing a derivative that's always 0. In stress.h:

    double n = coles_wake(1.0, pi_coles);

Contact Ivan or someone who knows better for help!