Re: IDASolve() failed with flag = -4, for a

[Alan Hindmarsh]

rather large set of equations
Sender: owner-sundials-users@...
Reply-To: sundials-users@...

Hello Francesco,

I don't see any actual errors in what you've supplied, but
there are two questionable things:

1. You have a comment about calling IDACalcIC, but no such call.
Inconsistent IC is the most likely cause of the initial failure you have.
It would be very useful to know if IDACalcIC succeeds with whatever you
give for initial guesses.  (It still might fail if the guesses are not
good enough.)

2. Your iout loop appears to generate tout = 1, 1, 1, 2, 3, ... .  Do
you have some reason for the extra calls with tout = 1?

-Alan H



On Wed, 14 Oct 2009, Francesco Pasqualini wrote:

 > 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
 >
 >
 >
 >
 >

----------------------------------------------------------------------------------------------------------------