Hi everybody,
I am working on the C++ implementation of an IVP problem described by
a DAE system composed of 153 equations (NEQ), 39 among them being
differential, and the relationships are non-linear. A Nightmare to
start with, and of course I am running into troubles with IDASolve().
I want to state that:
1) I've run the Robertson kinetics DAE serial example problem, and I
successfully modified some part of it
2) I got the same system of equation coded in Matlab, but in a tricky
way using ODE15s and a ODE formulation, with recalculation of all the
algebraic part (this code is not mine)
GENERAL INFORMATION:
I am on a MacOsX machine and I am using Xcode.
EXAMPLE OF THE CODE
I attach an extract of the code, with all the relevant part. For the
sake of simplicity I've removed very long list of numbers and the
actual residual function. If anyone wants to help I would be happy
sending the actual code.
The initial conditions are provided using the matlab code to generate
the set of values, so they shouldn't be the problem.
I suspect the problem could be with the formulation of the residual
function, but I would like an opinion, before going to find an error
in the formulation of so many equations.
static int resGPB(realtype time, N_Vector y, N_Vector ydot, N_Vector
resval,
void *user_data);
int main()
{
void *mem;
N_Vector y, ydot;
N_Vector idv, avtol;
int iout, retval;
realtype rtol, atol, t0, tout, tret;
y = N_VNew_Serial(NEQ);
if(check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
ydot = N_VNew_Serial(NEQ);
if(check_flag((void *)ydot, "N_VNew_Serial", 0)) return(1);
// avtol = N_VNew_Serial(NEQ);
// if(check_flag((void *)avtol, "N_VNew_Serial", 0)) return(1);
idv = N_VNew_Serial(NEQ);
if(check_flag((void *)idv, "id", 0)) return(1);
realtype id_vect[] = {...};
realtype ydot_vect[] = {...};
realtype IC[] = {...};
for(int i=1; i<=NEQ; i++){
Ith(idv,i)=id_vect[i-1];
Ith(y,i)=IC[i-1];
Ith(ydot,i) = ydot_vect[i-1];
// Ith(avtol,i) = ATOL;
}
GPBdata = (UserData) malloc(sizeof *GPBdata);
/* Set remaining inputs to IDAMalloc. */
t0 = 0;
rtol = 1.0e-5;
atol = 1.0e-12;
// all IDACreate and IDAMalloc to initialize IDA.
mem = IDACreate();
if(check_flag((void *)mem, "IDACreate", 0)) return(1);
// retval = IDASetUserData(mem, GPBdata);
// if(check_flag(&retval, "IDASetUserData", 1)) return(1);
retval = IDASetId(mem, idv);
if(check_flag(&retval, "IDASetId", 1)) return(1);
retval = IDAInit(mem, resGPB, t0, y, ydot);
if(check_flag(&retval, "IDAInit", 1)) return(1);
retval = IDASStolerances(mem, rtol, atol);
if(check_flag(&retval, "IDASStolerances", 1)) return(1);
retval = IDADense(mem, NEQ);
if(check_flag(&retval, "IDADense", 1)) return(1);
// retval = IDASetSuppressAlg(mem, TRUE);
// if(check_flag(&retval, "IDASetSuppressAlg", 1)) return(1);
// Call IDACalcIC (with default options) to correct the initial values.
tout = RCONST(1);
for (iout = 1; iout <= 10; iout++) {
retval = IDASolve(mem, tout, &tret, y, ydot, IDA_NORMAL);
if(check_flag(&retval, "IDASolve", 1)) return(retval);
PrintOutput(mem, tret, y);
if (iout < 3) tout *= 1; else tout += 1;
}
N_VDestroy_Serial(y);
N_VDestroy_Serial(ydot);
N_VDestroy_Serial(idv);
IDAFree(&mem);
free(GPBdata);
return(0);
}
/*
* Define the system residual function.
*/
static int resGPB(realtype time, N_Vector y, N_Vector ydot, N_Vector rr,
void *user_data){
...
}
ERROR
I obtain the error
[IDAS ERROR] IDASolve
At t = 0 and h = 0, the corrector convergence failed repeatedly or
with |h| = hmin.
SUNDIALS_ERROR: IDASolve() failed with flag = -4
Any help will be very much appreciated.
Thanks in advance
Francesco S. Pasqualini
