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

[Francesco Pasqualini-3]

rather large set of equations
Date: Thu, 15 Oct 2009 00:08:39 +0200
Sender: owner-sundials-users@...
Reply-To: sundials-users@...

Hi Alan,
 > I don't see any actual errors in what you've supplied
Great!

 > , but
 > there are two questionable things:
 >
uchh :D

 > 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.)
 >
Well, the point is that my system is likely to be not index one, I am
developing my understanding of the basics in these same time period,
so I can't be precise. My actual guess is that my systems is likely to
be in "The Hessenberg form of size three" and suspecting so, I am
following the advice of the manual "For problems that do not fall into
either of these categories, the user is responsible for passing
consistent values, or risk failure in the numerical integration.

For the consistency of the initial conditions, at this point, I've
proceeded like that: I initialized my matlab routine, with a very
small time step so that it performs only one iteration, and then I use
the vector generated after this iteration to feed the C++ routine.
They should be consistent.

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

That is a mistake I made, cause I've some macro defined and used in
that fragment, but in an attempt to simplify, I changed with 1.... in
fact the real number was 0.3, but actually the whole loop can be made
clearer: I must admit I've copied it from another examples :D

I attach another implementation where the loop is 1, 2, 3, ... , 10
and the error is still the original one.

What do you suggest as best candidate to be the problem?

Thank you,
Francesco


 >>
 >>

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, iout*tout, &tret, y, ydot,
IDA_NORMAL);
                  if(check_flag(&retval, "IDASolve", 1)) return(retval);

                  PrintOutput(mem, tret, y);

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

Francesco S. Pasqualini




-