<commit_before>/* -*-  mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */
/* -------------------------------------------------------------------------
ATS

License: see $ATS_DIR/COPYRIGHT
Author: <NAME>

EWCModelBase provides some of the functionality of EWCModel for inverse
evaluating.

------------------------------------------------------------------------- */

#include "ewc_model_base.hh"

#define DEBUG_FLAG 0

namespace Amanzi {

// ----------------------------------------------------------------------
// Lightweight wrapper to forward-evaluate the model.
// ----------------------------------------------------------------------
int EWCModelBase::Evaluate(double T, double p,
        double& energy, double& wc) {
  AmanziGeometry::Point res(2);
  int ierr = EvaluateEnergyAndWaterContent_(T,p,res);
  energy = res[0];
  wc = res[1];
  return ierr;
};


/* ----------------------------------------------------------------------
Solves a given energy and water content (at a given, fixed porosity), for
temperature and pressure.

Note this cannot really work for a saturated cell, as d_wc / d{T,p} = 0
<
Error codes:

  1 = Singular Jacobian at some point in the evaluation.  Often this is
      because the cell is saturated or becomes saturated and below
      freezing.
  2 = Iteration did not converge in max_steps (hard-coded to be 100 for
      now).
---------------------------------------------------------------------- */
int EWCModelBase::InverseEvaluate(double energy, double wc,
        double& T, double& p, bool verbose) {

  // -- scaling for the norms
  double wc_scale = 10.;
  double e_scale = 10000.;
  double T_corr_cap = 2.;
  double p_corr_cap = 200000.;
  double tol = 1.e-6;
  double max_steps = 100;
  double stepnum = 0;

  // get the initial residual
  AmanziGeometry::Point res(2);
  WhetStone::Tensor jac(2,2);
  int ierr = EvaluateEnergyAndWaterContentAndJacobian_(T,p,res,jac);
  if (ierr) {
    std::cout << "Error in evaluation: " << ierr << std::endl;
    return ierr + 10;
  }

  if (verbose) {
    std::cout << "Inverse Evaluating, e=" << energy << ", wc=" << wc << std::endl;
    std::cout << "   guess T,p (res) = " << T << ", " << p << " (" << res[0] << ", " << res[1] << ")" << std::endl;
  }

  AmanziGeometry::Point f(2);
  f[0] = energy;
  f[1] = wc;
  res = res - f;

  // check convergence
  AmanziGeometry::Point scaled_res(res);
  scaled_res[0] = res[0] / e_scale; scaled_res[1] = res[1] / wc_scale;
  double norm = AmanziGeometry::norm(scaled_res);

  bool converged = norm < tol;

  // workspace
  AmanziGeometry::Point x(2);
  AmanziGeometry::Point x_tmp(2);
  x[0] = T; x[1] = p;
  x_tmp[0] = T; x_tmp[1] = p;

  while (!converged) {
    // calculate the update size
    double detJ = jac.determinant();
    AmanziGeometry::Point correction;

    if (std::abs(detJ) < 1.e-20) {
      std::cout << " Zero determinant of Jacobian:" << std::endl;
      std::cout << "   [" << jac(0,0) << "," << jac(0,1) << "]" << std::endl;
      std::cout << "   [" << jac(1,0) << "," << jac(1,1) << "]" << std::endl;
      std::cout << "  at T,p = " << x_tmp[0] << ", " << x_tmp[1] << std::endl;
      std::cout << "  with res(e,wc) = " << res[0] << ", " << res[1] << std::endl;
      return 1;
    }


    jac.inverse();
    correction = jac * res;

    // cap the correction
    double scale = 1.;
    if (std::abs(correction[0]) > T_corr_cap) {
      scale = T_corr_cap / std::abs(correction[0]);
    }
    if (std::abs(correction[1]) > p_corr_cap) {
      double pscale = p_corr_cap / std::abs(correction[1]);
      scale = std::min(scale,pscale);
    }
    correction *= scale;

    // perform the update
    x_tmp = x - correction;
    ierr = EvaluateEnergyAndWaterContentAndJacobian_(x_tmp[0],x_tmp[1],res,jac);
    if (ierr) {
      std::cout << "Error in evaluation: " << ierr << std::endl;
      return ierr + 10;
    }
    res = res - f;

    // check convergence and damping
    scaled_res[0] = res[0] / e_scale; scaled_res[1] = res[1] / wc_scale;
    double norm_new = AmanziGeometry::norm(scaled_res);

    if (verbose) {
      std::cout << "  Iter: " << stepnum;
      std::cout << " corrected T,p (res) [norm] = " << x_tmp[0] << ", " << x_tmp[1] << " (" << res[0] << ", " << res[1] << ") ["
                << norm_new << "]" << std::endl;
    }

    double damp = 1.;
    bool backtracking_required = false;
    while (norm_new > norm) {
      backtracking_required = true;

      // backtrack
      damp *= 0.5;
      x_tmp = x - (damp * correction);

      // evaluate the damped value
      ierr = EvaluateEnergyAndWaterContent_(x_tmp[0],x_tmp[1],res);
      if (ierr) {
        std::cout << "Error in evaluation: " << ierr << std::endl;
        return ierr + 10;
      }
      res = res - f;

      // check the new residual
      scaled_res[0] = res[0] / e_scale; scaled_res[1] = res[1] / wc_scale;
      norm_new = AmanziGeometry::norm(scaled_res);

      if (verbose) {
        std::cout << "    Damping: " << stepnum;
        std::cout << " corrected T,p (res) [norm] = " << x_tmp[0] << ", " << x_tmp[1] << " (" << res[0] << ", " << res[1] << ") ["
                  << norm_new << "]" << std::endl;
      }

    }

    if (backtracking_required) {
      // must recalculate the Jacobian at the new value
      ierr = EvaluateEnergyAndWaterContentAndJacobian_(x_tmp[0],x_tmp[1],res,jac);
      if (ierr) {
        std::cout << "Error in evaluation: " << ierr << std::endl;
        return ierr + 10;
      }
      res = res - f;
    }

    // iterate
    x = x_tmp;
    norm = norm_new;

    AmanziGeometry::Point scaled_correction = damp * correction;
    scaled_correction[1] = scaled_correction[1] / 100000.;
    converged = norm < tol || AmanziGeometry::norm(scaled_correction) < 1.e-5;

    stepnum++;
    if (stepnum > max_steps && !converged) {
      std::cout << " Nonconverged after " << max_steps << " steps with norm (tol) "
                << norm << " (" << tol << ")" << std::endl;
      return 2;
    }
  }<commit_msg>updates convergence test for allowing to continue with smaller corrections<commit_after>/* -*-  mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */
/* -------------------------------------------------------------------------
ATS

License: see $ATS_DIR/COPYRIGHT
Author: <NAME>

EWCModelBase provides some of the functionality of EWCModel for inverse
evaluating.

------------------------------------------------------------------------- */

#include "ewc_model_base.hh"

#define DEBUG_FLAG 0

namespace Amanzi {

// ----------------------------------------------------------------------
// Lightweight wrapper to forward-evaluate the model.
// ----------------------------------------------------------------------
int EWCModelBase::Evaluate(double T, double p,
        double& energy, double& wc) {
  AmanziGeometry::Point res(2);
  int ierr = EvaluateEnergyAndWaterContent_(T,p,res);
  energy = res[0];
  wc = res[1];
  return ierr;
};


/* ----------------------------------------------------------------------
Solves a given energy and water content (at a given, fixed porosity), for
temperature and pressure.

Note this cannot really work for a saturated cell, as d_wc / d{T,p} = 0
<
Error codes:

  1 = Singular Jacobian at some point in the evaluation.  Often this is
      because the cell is saturated or becomes saturated and below
      freezing.
  2 = Iteration did not converge in max_steps (hard-coded to be 100 for
      now).
---------------------------------------------------------------------- */
int EWCModelBase::InverseEvaluate(double energy, double wc,
        double& T, double& p, bool verbose) {

  // -- scaling for the norms
  double wc_scale = 10.;
  double e_scale = 10000.;
  double T_corr_cap = 2.;
  double p_corr_cap = 200000.;
  double tol = 1.e-6;
  double max_steps = 100;
  double stepnum = 0;

  // get the initial residual
  AmanziGeometry::Point res(2);
  WhetStone::Tensor jac(2,2);
  int ierr = EvaluateEnergyAndWaterContentAndJacobian_(T,p,res,jac);
  if (ierr) {
    std::cout << "Error in evaluation: " << ierr << std::endl;
    return ierr + 10;
  }

  if (verbose) {
    std::cout << "Inverse Evaluating, e=" << energy << ", wc=" << wc << std::endl;
    std::cout << "   guess T,p (res) = " << T << ", " << p << " (" << res[0] << ", " << res[1] << ")" << std::endl;
  }

  AmanziGeometry::Point f(2);
  f[0] = energy;
  f[1] = wc;
  res = res - f;

  // check convergence
  AmanziGeometry::Point scaled_res(res);
  scaled_res[0] = res[0] / e_scale; scaled_res[1] = res[1] / wc_scale;
  double norm = AmanziGeometry::norm(scaled_res);

  bool converged = norm < tol;

  // workspace
  AmanziGeometry::Point x(2);
  AmanziGeometry::Point x_tmp(2);
  x[0] = T; x[1] = p;
  x_tmp[0] = T; x_tmp[1] = p;

  while (!converged) {
    // calculate the update size
    double detJ = jac.determinant();
    AmanziGeometry::Point correction;

    if (std::abs(detJ) < 1.e-20) {
      std::cout << " Zero determinant of Jacobian:" << std::endl;
      std::cout << "   [" << jac(0,0) << "," << jac(0,1) << "]" << std::endl;
      std::cout << "   [" << jac(1,0) << "," << jac(1,1) << "]" << std::endl;
      std::cout << "  at T,p = " << x_tmp[0] << ", " << x_tmp[1] << std::endl;
      std::cout << "  with res(e,wc) = " << res[0] << ", " << res[1] << std::endl;
      return 1;
    }


    jac.inverse();
    correction = jac * res;

    // cap the correction
    double scale = 1.;
    if (std::abs(correction[0]) > T_corr_cap) {
      scale = T_corr_cap / std::abs(correction[0]);
    }
    if (std::abs(correction[1]) > p_corr_cap) {
      double pscale = p_corr_cap / std::abs(correction[1]);
      scale = std::min(scale,pscale);
    }
    correction *= scale;

    // perform the update
    x_tmp = x - correction;
    ierr = EvaluateEnergyAndWaterContentAndJacobian_(x_tmp[0],x_tmp[1],res,jac);
    if (ierr) {
      std::cout << "Error in evaluation: " << ierr << std::endl;
      return ierr + 10;
    }
    res = res - f;

    // check convergence and damping
    scaled_res[0] = res[0] / e_scale; scaled_res[1] = res[1] / wc_scale;
    double norm_new = AmanziGeometry::norm(scaled_res);

    if (verbose) {
      std::cout << "  Iter: " << stepnum;
      std::cout << " corrected T,p (res) [norm] = " << x_tmp[0] << ", " << x_tmp[1] << " (" << res[0] << ", " << res[1] << ") ["
                << norm_new << "]" << std::endl;
    }

    double damp = 1.;
    bool backtracking_required = false;
    while (norm_new > norm) {
      backtracking_required = true;

      // backtrack
      damp *= 0.5;
      x_tmp = x - (damp * correction);

      // evaluate the damped value
      ierr = EvaluateEnergyAndWaterContent_(x_tmp[0],x_tmp[1],res);
      if (ierr) {
        std::cout << "Error in evaluation: " << ierr << std::endl;
        return ierr + 10;
      }
      res = res - f;

      // check the new residual
      scaled_res[0] = res[0] / e_scale; scaled_res[1] = res[1] / wc_scale;
      norm_new = AmanziGeometry::norm(scaled_res);

      if (verbose) {
        std::cout << "    Damping: " << stepnum;
        std::cout << " corrected T,p (res) [norm] = " << x_tmp[0] << ", " << x_tmp[1] << " (" << res[0] << ", " << res[1] << ") ["
                  << norm_new << "]" << std::endl;
      }

    }

    if (backtracking_required) {
      // must recalculate the Jacobian at the new value
      ierr = EvaluateEnergyAndWaterContentAndJacobian_(x_tmp[0],x_tmp[1],res,jac);
      if (ierr) {
        std::cout << "Error in evaluation: " << ierr << std::endl;
        return ierr + 10;
      }
      res = res - f;
    }

    // iterate
    x = x_tmp;
    norm = norm_new;

    AmanziGeometry::Point scaled_correction = damp * correction;
    scaled_correction[1] = scaled_correction[1] / 100000.;
    converged = norm < tol || AmanziGeometry::norm(scaled_correction) < 1.e-10;

    stepnum++;
    if (stepnum > max_steps && !converged) {
      std::cout << " Nonconverged after " << max_steps << " steps with norm (tol) "
                << norm << " (" << tol << ")" << std::endl;
      return 2;
    }
  }
{'fraction_non_alphanumeric': 0.22879814845428997, 'fraction_numerical': 0.018763431972226813, 'mean_word_length': 2.095957011258956, 'pattern_counts': {'":': 0, '<': 363, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3427427', 'n_tokens_mistral': 4288, 'n_tokens_neox': 4025, 'n_words': 985}