Commit b7e74b78 authored by Olivier Lantsoght's avatar Olivier Lantsoght
Browse files

[MbsEquil] Full code (and subcode) return error and always free memory, doc updated.

parent d0c7cf65
......@@ -13,6 +13,7 @@
*/
#include <stdio.h>
#include "mbs_errors_names.h"
#include "mbs_equil.h"
#include "mbs_project_interface.h"
#include "mbs_equil_struct.h"
......@@ -208,7 +209,7 @@ void mbs_delete_equil(MbsEquil* eq, MbsData* s)
int mbs_run_equil(MbsEquil* eq, MbsData* s)
{
int err;
int err = 0;
// 1. Initialize the simulation
// - - - - - - - - - - - - - -
err = mbs_equil_init(eq, s);
......@@ -218,13 +219,16 @@ int mbs_run_equil(MbsEquil* eq, MbsData* s)
// 2. Run the simulation
// - - - - - - - - - - -
mbs_equil_loop(eq, s); // another name might be more appropriate ??Q
err = mbs_equil_loop(eq, s); // another name might be more appropriate
if (err < 0){
return err;
}
// 3. Finish the simulation
// - - - - - - - - - - - -
mbs_equil_finish(eq, s);
err = mbs_equil_finish(eq, s);
return 0;
return err;
}
int mbs_equil_save(MbsEquil* eq, MbsData *s, int iter)
......@@ -594,29 +598,34 @@ int mbs_equil_init(MbsEquil* eq, MbsData* s)
return 0;
}
void mbs_equil_loop(MbsEquil* eq, MbsData* s)
int mbs_equil_loop(MbsEquil* eq, MbsData* s)
{
int err = 0;
switch (eq->options->method)
{
case 1:
{
mbs_equil_fsolvepk(&(mbs_equil_fct), eq, eq->aux, s);
err = mbs_equil_fsolvepk(&(mbs_equil_fct), eq, eq->aux, s);
break;
}
default:
{
err = _MBS_ERR_MOD_EQUIL + _MBS_ERR_INIT;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Invalid option value: MbsEquil->options->method=%d.\n", eq->options->method);
mbs_error_msg(">>EQUIL>> [%d] Choose a method listed in the documentation.\n", -200);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] Choose a method listed in the documentation.\n", err);
return(err);
}
}
return err;
}
void mbs_equil_finish(MbsEquil* eq, MbsData* s)
int mbs_equil_finish(MbsEquil* eq, MbsData* s)
{
int err = 0;
if (!eq){
return;
return err;
}
if (eq->options->save2file)
......@@ -667,16 +676,22 @@ void mbs_equil_finish(MbsEquil* eq, MbsData* s)
// memory has been released, raise error
if (err1 < 0 && err2 <0)
{
mbs_error_msg("[%d] and [%d] in mbs_equil_finish !! \n", -200 + err1, -200 + err2);
exit(1);
err = _MBS_ERR_MOD_EQUIL + err1; // arbitrary choice
mbs_error_msg("[%d] and [%d] in mbs_equil_finish !! \n", err, _MBS_ERR_MOD_EQUIL + err2);
return(err);
}
else if (err1 < 0){
mbs_error_msg("[%d] in mbs_equil_finish !! \n", -200 + err1);
exit(1);
err = _MBS_ERR_MOD_EQUIL + err1;
mbs_error_msg("[%d] in mbs_equil_finish !! \n", err);
return(err);
}
else if (err2 < 0){
mbs_error_msg("[%d] in mbs_equil_finish !! \n", -200 + err2);
exit(1);
err = _MBS_ERR_MOD_EQUIL + err2;
mbs_error_msg("[%d] in mbs_equil_finish !! \n", err);
return(err);
}
}
......@@ -710,6 +725,7 @@ void mbs_equil_finish(MbsEquil* eq, MbsData* s)
{
mbs_msg("\n>>EQUIL>> Equilibrium analysis finished.\n");
}
return err;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
......@@ -807,7 +823,7 @@ int mbs_equil_fct(double x[], int nx, double fx[], MbsEquil *eq, MbsAux *aux, Mb
return err;
}
void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*, MbsData*), MbsEquil *eq, MbsAux *aux, MbsData *s)
int mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*, MbsData*), MbsEquil *eq, MbsAux *aux, MbsData *s)
{
int i, j;
int err = 0;
......@@ -815,18 +831,19 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
void* grad_fun_ptr = NULL; // in case of Analythical expression for the gradient
int nx, nf;
double *xd, *f, *fd; // xd is used for different calculations.
double **grad;
double *xd = NULL, *f = NULL, *fd = NULL; // xd is used for different calculations.
double **grad = NULL;
double normR;
int r_cou;
double **grads = NULL;
double **grads;
double *xds, *xs, *pgrads, *fs, *xdsR, *xsR; // s stands for sensitive , R stands for Relaxation
// s stands for sensitive , R stands for Relaxation
double *xds = NULL, *xs = NULL, *pgrads = NULL, *fs = NULL, *xdsR = NULL, *xsR = NULL;
double norm1, norm2, inorm2, dnorm2, norm3, sc_prdct;
int *indx;
int *indx = NULL;
double d;
// 0. initialization (and allocation)
......@@ -853,10 +870,15 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
}
if (err < 0) // Error management
{
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error at initial gradient computation \n");
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
......@@ -907,9 +929,15 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
if (eq->nxs == 0)
{
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
err = _MBS_ERR_MOD_EQUIL + _MBS_ERR_MOD_SPEC_12;
mbs_msg(">>EQUIL>> mbs_equil_fsolvepk: none of the proposed parameter are sensitive !\n");
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk \n", -212);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk \n", err);
return(err);
}
// 2. Second initialization (and allocation)
......@@ -927,11 +955,26 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
err = (*fun_ptr)(eq->x, nx, f, eq, aux, s);
if (err < 0) // Error management
{
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error at initial function evaluatoin !\n");
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
user_equil_loop(s, eq);
eq->norm_pk = norm_dvec_0(f, nf);
// Save results of the initial step to the buffer
......@@ -943,10 +986,24 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
}
if (err < 0) // Error management
{
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error during files saving !\n");
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
eq->iter = 0;
......@@ -982,10 +1039,24 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
}
if (err < 0) // Error management
{
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error during equilibrium process at iteration %d !\n", eq->iter);
mbs_error_msg("[%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg("[%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
}
// Gradient Resizing (sensitive variables only) : [nfs * nx]
......@@ -1004,10 +1075,26 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
err = ludcmp_0(grads, eq->nxs, indx, &d); // those functions make the matrix division (inversion)
if (err < 0) // Error management
{
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
free_ivec_0(indx);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error during equilibrium process at iteration %d !\n", eq->iter);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
lubksb_0(grads, eq->nxs, indx, xds);
free_ivec_0(indx);
......@@ -1097,10 +1184,24 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
err = (*fun_ptr)(eq->x, nx, f, eq, aux, s);
if (err < 0) // Error management
{
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error during relaxation step %d for equilibrium iteration %d !\n", r_cou, eq->iter);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return err;
}
eq->norm_pk = norm_dvec_0(f, nf);
......@@ -1137,31 +1238,72 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
}
if (err < 0) // Error management
{
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error during files saving !\n");
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
}
// the equil function error is not caught directly to allow the relaxation process to do something
if (err < 0) // Error management
{
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
err += _MBS_ERR_MOD_EQUIL;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Equilibrium function error during equilibrium process at iteration %d !\n", eq->iter);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
}
if (eq->iter >= eq->options->itermax) // error management
{
err = -11;
free_dmat_0(grads);
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_dvec_0(xds);
free_dvec_0(xs);
free_dvec_0(pgrads);
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
err = _MBS_ERR_MOD_EQUIL + _MBS_ERR_MOD_SPEC_11;
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Equilibrium process has not converged after %d iterations !\n", eq->options->itermax);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", -200 + err);
exit(1);
mbs_error_msg(">>EQUIL>> [%d] in mbs_equil_fsolvepk !! \n", err);
return(err);
}
else
{
......@@ -1187,6 +1329,8 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
free_dvec_0(fs);
free_dvec_0(xdsR);
free_dvec_0(xsR);
return(err);
}
int mbs_equil_grad_lpk(MbsEquil *eq, MbsData *s)
......@@ -1279,8 +1423,8 @@ int mbs_equil_grad_dev(MbsEquil *eq, MbsData *s)
int err = 0;
int isSensitive;
int nx, nf;
double *xd, *f, *fd; // xd is used for different calculations.
int *xs, *xns;
double *xd = NULL, *f = NULL, *fd = NULL; // xd is used for different calculations.
int *xs = NULL, *xns = NULL;
int nxs, nxns;
// 0. initialization (and allocation)
......@@ -1298,6 +1442,13 @@ int mbs_equil_grad_dev(MbsEquil *eq, MbsData *s)
err = mbs_equil_fct(eq->x, nx, f, eq, eq->aux, s);
if (err < 0) // Error management
{
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_ivec_0(xs);
free_ivec_0(xns);
mbs_msg("\t in EQUIL GRAD, computation of grad_Rr = delta Rr/delta x \n");
return err;
}
......@@ -1316,6 +1467,13 @@ int mbs_equil_grad_dev(MbsEquil *eq, MbsData *s)
err = mbs_equil_fct(xd, nx, fd, eq, eq->aux, s);
if (err < 0) // Error management
{
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_ivec_0(xs);
free_ivec_0(xns);
mbs_msg("\t in EQUIL GRAD, computation of grad_Rr = delta Rr/delta x \n");
return err;
}
......
......@@ -71,8 +71,13 @@ int mbs_equil_init(MbsEquil* mbs_equil,MbsData* mbs_data);
/**
* Process the equibrium : find the equilibrium variables, x that satisfy f(x)= [Fruc ; fxe] = 0
* Based on the solvepk algorithm (see Numerics)
*
* \param[in, out] mbs_equil the MbsEquil to be run.
* \param[in, out] mbs_data the MbsData structure of the model for which the equilibrium is computed.
*
* \return Error status, <0 in case of failure.
*/
void mbs_equil_loop(MbsEquil* mbs_equil,MbsData* mbs_data);
int mbs_equil_loop(MbsEquil* mbs_equil,MbsData* mbs_data);
/**
* Set equilibrum flag to done in MbsData structure
......@@ -80,7 +85,7 @@ void mbs_equil_loop(MbsEquil* mbs_equil,MbsData* mbs_data);
* offer the possibility for user to finish operation through the use of user_equil_finish
*
*/
void mbs_equil_finish(MbsEquil* mbs_equil,MbsData* mbs_data);
int mbs_equil_finish(MbsEquil* mbs_equil,MbsData* mbs_data);
/**
......@@ -164,32 +169,41 @@ void mbs_print_equil(MbsEquil* eq);
//-----------------------------------------------------------------
/**
* Solve equation f(x)=0 based on a given initial value for the variables, x0
* Algorithm developped by Paul Fisette and ?Q
* \brief Solve equation f(x)=0 based on given initial values.
*
* Algorithm developped by Paul Fisette and ...
*
* \p fun_ptr vector containing the expression of f(x) as function of x
* \p eq ?Q to change for a more general use of function of fsolvepk
* \p aux ?Q to change for a more general use of function of fsolvepk
* \p s ?Q to change for a more general use of function of fsolvepk
* \param[in] fun_ptr pointer to the error function f(x) to be solved.
* \param[in, out] eq the MbsEquil to be run.
* \param[in, out] aux the MbsAux structure related to the model for which the
* equilibrium is computed.
* \param[in, out] s the MbsData structure of the model for which the equilibrium
* is computed.
*
* \return Error status, <0 in case of failure.
*/
void mbs_equil_fsolvepk(int (*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*, MbsData*), MbsEquil *eq, MbsAux *aux, MbsData *s );
int mbs_equil_fsolvepk(int (*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*, MbsData*), MbsEquil *eq, MbsAux *aux, MbsData *s );
/**
* Compute the gradient of f(x) for a given configuration x based on the lpk parabolic fitting method
*
* \p eq ?Q to change for a more general use of function of fsolvepk
* \p s ?Q to change for a more general use of function of fsolvepk
*/
* \brief Compute the gradient of f(x) for a given configuration based on the lpk parabolic fitting method.
*
* \param[in, out] eq the MbsEquil to be run.
* \param[in, out] s the MbsData structure of the model for which the equilibrium
* is computed.
*
* \return Error status, <0 in case of failure.
*/
int mbs_equil_grad_lpk(MbsEquil *eq, MbsData *s);
/**
* Compute the gradient of f(x) for a given configuration x based on a deviation computation
*
* \p eq ?Q to change for a more general use of function of fsolvepk
* \p s ?Q to change for a more general use of function of fsolvepk
* \brief Compute the gradient of f(x) for a given configuration based on a deviation computation.
*
* \param[in, out] eq the MbsEquil to be run.
* \param[in, out] s the MbsData structure of the model for which the equilibrium
* is computed.
*
* \return Error status, <0 in case of failure.
*/
int mbs_equil_grad_dev(MbsEquil *eq, MbsData *s);
#endif
......@@ -9,6 +9,7 @@
// Last update : 14/06/2016
//---------------------------
#include "mbs_errors_names.h"
#include "mbs_linearipk.h"
#include "mbs_message.h"
#include "mbs_project_interface.h"
......@@ -243,9 +244,10 @@ int mbs_linearipk(double **GK, MbsLpk *lpk, MbsAux *aux, MbsData *s, int EqChoic
// waiting column k is still not acepted
else // Error management
{
err = _MBS_ERR_MOD_SPEC_10;
mbs_msg(">>LINPK>> Column %d/%d no convergence after %d iterations\n", lpk->diverge_ind[k] + 1, lpk->nx, lpk->itermax);
mbs_msg(">>LINPK>> Restart with other parameters !\n");
err = -10;
return err;
}
}
......@@ -290,8 +292,10 @@ int mbs_lineari_fct(double *x, double *Fx, MbsLpk *lpk, MbsAux *aux, MbsData *s,
{
if (isnan(s->uxd[i + 1])) // Error management
{
err = _MBS_ERR_MID_UXD + _MBS_ERR_LOW_NAN;
mbs_msg(">>LINPK>> User derivative uxd[%d] is Nan\n", i + 1);
return -47;
return err;
}
Fx[i] = s->uxd[i + 1];
}
......@@ -304,14 +308,18 @@ int mbs_lineari_fct(double *x, double *Fx, MbsLpk *lpk, MbsAux *aux, MbsData *s,
{
if (isnan(Fx[i])) // Error management
{
err = _MBS_ERR_MOD_SPEC_17;
mbs_msg(">>LINPK>> User equil Fxe[%d] is Nan\n", i + 1);
return -17; // equil specific code
return err;
}
}
}
else{
err = _MBS_ERR_MOD_SPEC_16;
mbs_msg(">>LINPK>> Invalid equation choice: %d\n", EqChoice);
return -16; // equil specific code
return err; // equil specific code
}
return err;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment