gem_ipopt.c

Go to the documentation of this file.

/* src/gem_ipopt.c
 *
 * Copyright (C) 2011-2018 Dongdong Li
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundataion; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABLITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include "ipopt/IpStdCInterface.h"
#include <stdlib.h>
#include <assert.h>
#include <stdio.h>

extern int gem_ipopt ();

static Bool eval_f (Index n, Number* x, Bool new_x,
            Number* obj_value, UserDataPtr user_data);

static Bool eval_grad_f (Index n, Number* x, Bool new_x,
                 Number* grad_f, UserDataPtr user_data);

static Bool eval_g (Index n, Number* x, Bool new_x,
            Index m, Number* g, UserDataPtr user_data);

static Bool eval_jac_g (Index n, Number *x, Bool new_x,
                Index m, Index nele_jac,
                Index *iRow, Index *jCol, Number *values,
                UserDataPtr user_data);

static Bool eval_h (Index n, Number *x, Bool new_x, Number obj_factor,
            Index m, Number *lambda, Bool new_lambda,
            Index nele_hess, Index *iRow, Index *jCol,
            Number *values, UserDataPtr user_data);

static Bool intermediate_cb (Index alg_mod, Index iter_count, Number obj_value,
                     Number inf_pr, Number inf_du, Number mu, Number d_norm,
                     Number regularization_size, Number alpha_du,
                     Number alpha_pr, Index ls_trials, UserDataPtr user_data);

extern int gem_ipopt()
{
  Index n=-1;                          /* number of variables */
  Index m=-1;                          /* number of constraints */
  Number* x_L = NULL;                  /* lower bounds on x */
  Number* x_U = NULL;                  /* upper bounds on x */
  Number* g_L = NULL;                  /* lower bounds on g */
  Number* g_U = NULL;                  /* upper bounds on g */
  IpoptProblem nlp = NULL;             /* IpoptProblem */
  enum ApplicationReturnStatus status; /* Solve return code */
  Number* x = NULL;                    /* starting point and solution vector */
  Number* mult_g = NULL;               /* constraint multipliers
                      at the solution */
  Number* mult_x_L = NULL;             /* lower bound multipliers
                      at the solution */
  Number* mult_x_U = NULL;             /* upper bound multipliers
                      at the solution */
  Number obj;                          /* objective value */
  Index i;                             /* generic counter */

  /* Number of nonzeros in the Jacobian of the constraints */
  Index nele_jac = (total_comp_num + 1) * total_species_num;
  /* Number of nonzeros in the Hessian of the Lagrangian (lower or
     upper triangual part only) */
  Index nele_hess = total_species_num * total_species_num;
  /* indexing style for matrices */
  Index index_style = 0; /* C-style; start counting of rows and column
                indices at 0 */

  /* set the number of variables and allocate space for the bounds */
  n = total_species_num;
  x_L = (Number*)malloc(sizeof(Number) * n);
  x_U = (Number*)malloc(sizeof(Number) * n);

  /* set the values for the variable bounds */
  for (i = 0; i < n; i++) 
    {
      x_L[i] = 0.0;
      x_U[i] = 2.0e+19;
    }

  /* set the number of constraints and allocate space for the bounds */
  m = total_comp_num + 1; // mass balance + charge blance
  g_L = (Number*)malloc(sizeof(Number) * m);
  g_U = (Number*)malloc(sizeof(Number) * m);
  /* set the values of the constraint bounds */

  for (i = 0; i < m; i++)
    {
      g_L[i] = 0;
      g_U[i] = 0;
    }

  /* create the IpoptProblem */
  nlp = CreateIpoptProblem(n, x_L, x_U, m, g_L, g_U, nele_jac, nele_hess,
                           index_style, &eval_f, &eval_g, &eval_grad_f,
                           &eval_jac_g, &eval_h);

  /* We can free the memory now - the values for the bounds have been
     copied internally in CreateIpoptProblem */
  free(x_L);
  free(x_U);
  free(g_L);
  free(g_U);

  /* Set some options.  Note the following ones are only examples,
     they might not be suitable for your problem. */
  AddIpoptNumOption(nlp, "tol", 1e-10);
  AddIpoptStrOption(nlp, "mu_strategy", "adaptive");
  AddIpoptStrOption(nlp, "hessian_approximation", "limited-memory");
  AddIpoptIntOption(nlp, "print_level", 0);
  

  /* allocate space for the initial point and set the values */
  x = (Number*)malloc(sizeof(Number) * n);

  for (i = 0; i < n; i ++)
    {
      x[i] = 1.0;
    }

  /* allocate space to store the bound multipliers at the solution */
  mult_g = (Number*)malloc(sizeof(Number) * m);
  mult_x_L = (Number*)malloc(sizeof(Number) * n);
  mult_x_U = (Number*)malloc(sizeof(Number) * n);

  /* Set the callback method for intermediate user-control.  This is
   * not required, just gives you some intermediate control in case
   * you need it. */
  /* SetIntermediateCallback(nlp, intermediate_cb); */

  /* solve the problem */
  status = IpoptSolve(nlp, x, NULL, &obj, mult_g, mult_x_L, mult_x_U, NULL);

  if (status == Solve_Succeeded) 
    {
      printf("\nSUCCESSED IPOPT OPTIMIZATION.\n");
      successed = TRUE;
    }
  else 
    {
      printf("\nERROR OCCURRED DURING IPOPT OPTIMIZATION.\n");
      successed = FALSE;
    }

  if (status == Solve_Succeeded) 
    {
      n_eq = (double *)malloc (sizeof(double) * n);

      for (i=0; i < n; i++) 
    {
      n_eq[i] = x[i];
    }

      gibbs_eq = obj;

      successed = TRUE;
    }

  /* free allocated memory */
  FreeIpoptProblem(nlp);
  free(x);
  free(mult_g);
  free(mult_x_L);
  free(mult_x_U);

  return (int)status;
}

/* Function Implementations */
static Bool eval_f(Index n, Number* x, Bool new_x,
            Number* obj_value, UserDataPtr user_data)
{
  assert(n == total_species_num);

  *obj_value = 0;

  double t = system_T/1000.0;
  double T = system_T;
  double P = system_P;

  int i, j;
  int species_num_index = 0;

  //Total Gibbs energy of gas phases
  for (i = 0; i < gas_phase_num; i++)
    {
      double sumn = 0;
      double ngi[phases.gas[i].species_num];

      for (j = 0; j < phases.gas[i].species_num; j++)
    {
      ngi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.gas[i].species_num; j++)
    {
      sumn += ngi[j];
    }

      for (j = 0; j < phases.gas[i].species_num; j++)
    {
      double A = phases.gas[i].gas_species[j].thermo_coeff[0];
      double B = phases.gas[i].gas_species[j].thermo_coeff[1];
      double C = phases.gas[i].gas_species[j].thermo_coeff[2];
      double D = phases.gas[i].gas_species[j].thermo_coeff[3];
      double E = phases.gas[i].gas_species[j].thermo_coeff[4];
      double F = phases.gas[i].gas_species[j].thermo_coeff[5];
      double G = phases.gas[i].gas_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      // kJ
      *obj_value += ngi[j] * (G0 + R * T * log (ngi[j]/sumn) / 1000);

      species_num_index ++;
    }
    }

  //Total Gibbs energy of aqueous phases
  for (i = 0; i < aqueous_phase_num; i++)
    {
      double sumn = 0;
      double naqi[phases.aqueous[i].species_num];

      for (j = 0; j < phases.aqueous[i].species_num; j++)
    {
      naqi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.aqueous[i].species_num; j++)
    {
      sumn += naqi[j];
    }

      for (j = 0; j < phases.aqueous[i].species_num; j++)
    {
      double A = phases.aqueous[i].aqueous_species[j].thermo_coeff[0];
      double B = phases.aqueous[i].aqueous_species[j].thermo_coeff[1];
      double C = phases.aqueous[i].aqueous_species[j].thermo_coeff[2];
      double D = phases.aqueous[i].aqueous_species[j].thermo_coeff[3];
      double E = phases.aqueous[i].aqueous_species[j].thermo_coeff[4];
      double F = phases.aqueous[i].aqueous_species[j].thermo_coeff[5];
      double G = phases.aqueous[i].aqueous_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      double aj = psc_a (phases.aqueous[i], j, naqi, T, P);

      *obj_value += naqi[j] * (G0 + R * T * log (aj) / 1000);

      species_num_index ++;
    } 
    }

  //Total Gibbs energy of solid phases
  for (i = 0; i < solid_phase_num; i++)
    {
      double sumn = 0;
      double nsi[phases.solids[i].species_num];

      for (j = 0; j < phases.solids[i].species_num; j++)
    {
      nsi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.solids[i].species_num; j++)
    {
      double A = phases.solids[i].solid_species[j].thermo_coeff[0];
      double B = phases.solids[i].solid_species[j].thermo_coeff[1];
      double C = phases.solids[i].solid_species[j].thermo_coeff[2];
      double D = phases.solids[i].solid_species[j].thermo_coeff[3];
      double E = phases.solids[i].solid_species[j].thermo_coeff[4];
      double F = phases.solids[i].solid_species[j].thermo_coeff[5];
      double G = phases.solids[i].solid_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      *obj_value += nsi[j] * G0;

      species_num_index ++;
    }
    }

  //Total Gibbs energy of solid solution phases
  for (i = 0; i < solidsolution_phase_num; i++)
    {
      double sumn = 0;
      double nssi[phases.solidsolution[i].species_num];

      for (j = 0; j < phases.solidsolution[i].species_num; j++)
    {
      nssi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.solidsolution[i].species_num; j++)
    {
      sumn += nssi[j];
    }

      for (j = 0; j < phases.solidsolution[i].species_num; j++)
    {
      double A = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[0];
      double B = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[1];
      double C = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[2];
      double D = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[3];
      double E = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[4];
      double F = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[5];
      double G = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      double aj = rkg_a (phases.solidsolution[i], j, nssi, T, P);

      *obj_value += nssi[j] * (G0 + R * T * log (aj) / 1000);

      species_num_index ++;
    }
    }

  return TRUE;
}

static Bool eval_grad_f(Index n, Number* x, Bool new_x,
                 Number* grad_f, UserDataPtr user_data)
{
  assert(n == total_species_num);

  double t = system_T/1000.0;
  double T = system_T;
  double P = system_P;

  int species_num_index = 0;

  int i, j;

  // Gas phases
  for (i = 0; i < gas_phase_num; i ++)
    {
      double sumn = 0;
      double ngi[phases.gas[i].species_num];

      for (j = 0; j < phases.gas[i].species_num; j++)
    {
      ngi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.gas[i].species_num; j++)
    {
      sumn += ngi[j];
    }

      for (j = 0; j < phases.gas[i].species_num; j ++)
    {
      double A = phases.gas[i].gas_species[j].thermo_coeff[0];
      double B = phases.gas[i].gas_species[j].thermo_coeff[1];
      double C = phases.gas[i].gas_species[j].thermo_coeff[2];
      double D = phases.gas[i].gas_species[j].thermo_coeff[3];
      double E = phases.gas[i].gas_species[j].thermo_coeff[4];
      double F = phases.gas[i].gas_species[j].thermo_coeff[5];
      double G = phases.gas[i].gas_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      grad_f[species_num_index] = (G0 + R * T * log (ngi[j]/sumn) / 1000);

      species_num_index ++;
    }
    }

  // Aqueous phases
  for (i = 0; i < aqueous_phase_num; i ++)
    {
      double sumn = 0;
      double naqi[phases.aqueous[i].species_num];

      for (j = 0; j < phases.aqueous[i].species_num; j++)
    {
      naqi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.aqueous[i].species_num; j++)
    {
      sumn += naqi[j];
    }

      for (j = 0; j < phases.aqueous[i].species_num; j ++)
    {
      double A = phases.aqueous[i].aqueous_species[j].thermo_coeff[0];
      double B = phases.aqueous[i].aqueous_species[j].thermo_coeff[1];
      double C = phases.aqueous[i].aqueous_species[j].thermo_coeff[2];
      double D = phases.aqueous[i].aqueous_species[j].thermo_coeff[3];
      double E = phases.aqueous[i].aqueous_species[j].thermo_coeff[4];
      double F = phases.aqueous[i].aqueous_species[j].thermo_coeff[5];
      double G = phases.aqueous[i].aqueous_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      double aj = psc_a (phases.aqueous[i], j, naqi, T, P);

      grad_f[species_num_index] = (G0 + R * T * log (aj) / 1000);

      species_num_index ++;
    }
    }

  // Solid phases;
  for (i = 0; i < solid_phase_num; i ++)
    {
      double sumn = 0;
      double nsi[phases.solids[i].species_num];

      for (j = 0; j < phases.solids[i].species_num; j++)
    {
      nsi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.solids[i].species_num; j++)
    {
      sumn += nsi[j];
    }

      for (j = 0; j < phases.solids[i].species_num; j ++)
    {
      double A = phases.solids[i].solid_species[j].thermo_coeff[0];
      double B = phases.solids[i].solid_species[j].thermo_coeff[1];
      double C = phases.solids[i].solid_species[j].thermo_coeff[2];
      double D = phases.solids[i].solid_species[j].thermo_coeff[3];
      double E = phases.solids[i].solid_species[j].thermo_coeff[4];
      double F = phases.solids[i].solid_species[j].thermo_coeff[5];
      double G = phases.solids[i].solid_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      grad_f[species_num_index] = (G0);

      species_num_index ++;
    }
    }

  // Solid solid phases;
  for (i = 0; i < solidsolution_phase_num; i ++)
    {
      double sumn = 0;
      double nssi[phases.solidsolution[i].species_num];

      for (j = 0; j < phases.solidsolution[i].species_num; j++)
    {
      nssi[j] = x[j + species_num_index];
    }

      for (j = 0; j < phases.solidsolution[i].species_num; j++)
    {
      sumn += nssi[j];
    }

      for (j = 0; j < phases.solidsolution[i].species_num; j ++)
    {
      double A = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[0];
      double B = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[1];
      double C = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[2];
      double D = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[3];
      double E = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[4];
      double F = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[5];
      double G = phases.solidsolution[i].solidsolution_species[j].thermo_coeff[6];

      double G0 = A * (t - t*log(t)) - B/2.0 *t*t - C/6.0 *t*t*t -
        D/12.0 * t*t*t*t - E/2.0/t + F - G * t;

      double aj = rkg_ai (phases.solidsolution[i], j, nssi, T, P);

      grad_f[species_num_index] = (G0 + R * T * log(aj) / 1000);

      species_num_index ++;
    }
    }

  return TRUE;
}

static Bool eval_g(Index n, Number* x, Bool new_x,
            Index m, Number* g, UserDataPtr user_data)
{
  assert(n == total_species_num);
  assert(m == total_comp_num + 1);

  int i, j, k;

  // mass balance
  for (i = 0; i < m - 1; i ++)
    {
      g[i] = 0;

      int species_num_index = 0;

      // Gas phases
      for(j = 0; j < gas_phase_num; j ++)
    {
      for (k = 0; k < phases.gas[j].species_num; k ++)
        {
          g[i] += phases.gas[j].gas_species[k].elements[i] * x[species_num_index];

          species_num_index ++;
        }
    }

      // Aqueous phases
      for(j = 0; j < aqueous_phase_num; j ++)
    {
      for (k = 0; k < phases.aqueous[j].species_num; k ++)
        {
          g[i] += phases.aqueous[j].aqueous_species[k].elements[i] * x[species_num_index];

          species_num_index ++;
        }
    }

      // Solid phases
      for(j = 0; j < solid_phase_num; j ++)
    {
      for (k = 0; k < phases.solids[j].species_num; k ++)
        {
          g[i] += phases.solids[j].solid_species[k].elements[i] * x[species_num_index];

          species_num_index ++;
        }
    }

      // Solid solution phases
      for(j = 0; j < solidsolution_phase_num; j ++)
    {
      for (k = 0; k < phases.solidsolution[j].species_num; k ++)
        {
          g[i] += phases.solidsolution[j].solidsolution_species[k].elements[i] * x[species_num_index];

          species_num_index ++;
        }
    }

      g[i] -= total_components[i];
    }

  // charge balance
  {
    g[i] = 0;

    int species_num_index = 0;

    // Gas phases
    for(j = 0; j < gas_phase_num; j ++)
      {
    for (k = 0; k < phases.gas[j].species_num; k ++)
      {
        g[i] += phases.gas[j].gas_species[k].charge * x[species_num_index];

        species_num_index ++;
      }
      }

    // Aqueous phases
    for(j = 0; j < aqueous_phase_num; j ++)
      {
    for (k = 0; k < phases.aqueous[j].species_num; k ++)
      {
        g[i] += phases.aqueous[j].aqueous_species[k].charge * x[species_num_index];

        species_num_index ++;
      }
      }

    // Solid phases
    for(j = 0; j < solid_phase_num; j ++)
      {
    for (k = 0; k < phases.solids[j].species_num; k ++)
      {
        g[i] += phases.solids[j].solid_species[k].charge * x[species_num_index];

        species_num_index ++;
      }
      }

    // Solid solution phases
    for(j = 0; j < solidsolution_phase_num; j ++)
      {
    for (k = 0; k < phases.solidsolution[j].species_num; k ++)
      {
        g[i] += phases.solidsolution[j].solidsolution_species[k].charge * x[species_num_index];

        species_num_index ++;
      }
      }

    g[i] -= system_charge;
  }

  return TRUE;
}

static Bool eval_jac_g(Index n, Number *x, Bool new_x,
                Index m, Index nele_jac,
                Index *iRow, Index *jCol, Number *values,
                UserDataPtr user_data)
{
  int i, j, k;

  if (values == NULL) 
    {
      /* return the structure of the jacobian */

      /* this particular jacobian is dense */
      for (i = 0; i < m; i ++)
    {
      for (j = 0; j < n; j ++)
        {
          iRow[i*n + j] = i;
          jCol[i*n + j] = j;
        }
    }
    }
  else 
    {
      /* return the values of the jacobian of the constraints */

      int species_num_index = 0;
      // Gas species
      for (i = 0; i < gas_phase_num; i ++)
    {
      for (j = 0; j < phases.gas[i].species_num; j ++)
        {
          for (k = 0; k < m - 1; k ++)
        {
          //[component i, species j]
          values[k * total_species_num + species_num_index] = phases.gas[i].gas_species[j].elements[k];
        }

          values[(m - 1) * total_species_num + species_num_index] = phases.gas[i].gas_species[j].charge;

          species_num_index ++;
        }
    }

      // aqueous species
      for (i = 0; i < aqueous_phase_num; i ++)
    {
      for (j = 0; j < phases.aqueous[i].species_num; j ++)
        {
          for (k = 0; k < m - 1; k ++)
        {
          values[k * total_species_num + species_num_index] = phases.aqueous[i].aqueous_species[j].elements[k];
        }

          values[(m - 1) * total_species_num + species_num_index] = phases.aqueous[i].aqueous_species[j].charge;

          species_num_index ++;
        }
    }

      // Solid species
      for (i = 0; i < solid_phase_num; i ++)
    {
      for (j = 0; j < phases.solids[i].species_num; j ++)
        {
          for (k = 0; k < total_comp_num; k ++)
        {
          values[k * total_species_num + species_num_index] = phases.solids[i].solid_species[j].elements[k];
        }

          values[k * total_species_num + species_num_index] = phases.solids[i].solid_species[j].charge;

          species_num_index ++;
        }
    }

      // Solid solution species
      for (i = 0; i < solidsolution_phase_num; i ++)
    {
      for (j = 0; j < phases.solidsolution[i].species_num; j ++)
        {
          for (k = 0; k < total_comp_num; k ++)
        {
          values[k * total_species_num + species_num_index] = phases.solidsolution[i].solidsolution_species[j].elements[k];
        }

          values[k * total_species_num + species_num_index] = phases.solidsolution[i].solidsolution_species[j].charge;

          species_num_index ++;
        }
    }
    }

  return TRUE;
}

static Bool eval_h(Index n, Number *x, Bool new_x, Number obj_factor,
            Index m, Number *lambda, Bool new_lambda,
            Index nele_hess, Index *iRow, Index *jCol,
            Number *values, UserDataPtr user_data)
{
  return FALSE;
}

static Bool intermediate_cb(Index alg_mod, Index iter_count, Number obj_value,
                     Number inf_pr, Number inf_du, Number mu, Number d_norm,
                     Number regularization_size, Number alpha_du,
                     Number alpha_pr, Index ls_trials, UserDataPtr user_data)
{
  printf("Testing intermediate callback in iteration %d\n", iter_count);
  if (inf_pr < 1e-4) return FALSE;

  return TRUE;
}