/[ascend]/trunk/base/generic/solver/ida.c
ViewVC logotype

Diff of /trunk/base/generic/solver/ida.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 945 by johnpye, Sat Nov 25 12:41:03 2006 UTC revision 946 by johnpye, Sat Nov 25 15:28:56 2006 UTC
# Line 225  int integrator_ida_params_default(Integr Line 225  int integrator_ida_params_default(Integr
225              ,(SlvParameterInitChar){{"linsolver"              ,(SlvParameterInitChar){{"linsolver"
226              ,"Linear solver",1              ,"Linear solver",1
227              ,"See IDA manual, section 5.5.3."              ,"See IDA manual, section 5.5.3."
228          }, "DENSE"}, (char *[]){"DENSE","SPGMR",NULL}          }, "SPGMR"}, (char *[]){"DENSE","SPGMR",NULL}
229      );      );
230    
231      asc_assert(p->num_parms == IDA_PARAMS_SIZE);      asc_assert(p->num_parms == IDA_PARAMS_SIZE);
# Line 329  int integrator_ida_solve( Line 329  int integrator_ida_solve(
329      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));      IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
330      /* there's no capability for setting *minimum* step size in IDA */      /* there's no capability for setting *minimum* step size in IDA */
331    
     CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");  
332    
333      /* attach linear solver module, using the default value of maxl */      /* attach linear solver module, using the default value of maxl */
334      linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);      linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
335        CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
336      if(strcmp(linsolver,"SPGMR")==0){      if(strcmp(linsolver,"SPGMR")==0){
337            CONSOLE_DEBUG("USING 'SCALED PRECONDITIONER GMRES' LINEAR SOLVER");
338          flag = IDASpgmr(ida_mem, 0);          flag = IDASpgmr(ida_mem, 0);
339          if(flag==IDASPILS_MEM_NULL){          if(flag==IDASPILS_MEM_NULL){
340              ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");              ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
# Line 358  int integrator_ida_solve( Line 359  int integrator_ida_solve(
359              CONSOLE_DEBUG("USING NUMERICAL DIFF");              CONSOLE_DEBUG("USING NUMERICAL DIFF");
360          }          }
361      }else if(strcmp(linsolver,"DENSE")==0){      }else if(strcmp(linsolver,"DENSE")==0){
362          CONSOLE_DEBUG("USING IDADEENSE SOLVER, size = %d",size);          CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
363          flag = IDADense(ida_mem, size);          flag = IDADense(ida_mem, size);
364          switch(flag){          switch(flag){
365              case IDADENSE_SUCCESS: break;              case IDADENSE_SUCCESS: break;
# Line 557  int integrator_ida_djex(long int Neq, re Line 558  int integrator_ida_djex(long int Neq, re
558  ){  ){
559      IntegratorSystem *blsys;      IntegratorSystem *blsys;
560      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
561      realtype *yval;      char *relname, *varname;
562        int status;
563        struct rel_relation **relptr;
564        int i;
565        var_filter_t filter = {VAR_SVAR, VAR_SVAR};
566        double *derivatives;
567        int *variables;
568        int count, j, var_yindex;
569    
570      blsys = (IntegratorSystem *)jac_data;      blsys = (IntegratorSystem *)jac_data;
571      enginedata = integrator_ida_enginedata(blsys);      enginedata = integrator_ida_enginedata(blsys);
572    
573        /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
574        variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
575        derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
576    
577      /* pass the values of everything back to the compiler */      /* pass the values of everything back to the compiler */
578      integrator_set_t(blsys, (double)tt);      integrator_set_t(blsys, (double)tt);
579      integrator_set_y(blsys, NV_DATA_S(yy));      integrator_set_y(blsys, NV_DATA_S(yy));
580      integrator_set_ydot(blsys, NV_DATA_S(yp));      integrator_set_ydot(blsys, NV_DATA_S(yp));
581    
582      yval = NV_DATA_S(yy);      /* print vars */
583        for(i=0; i < blsys->n_y; ++i){
584            varname = var_make_name(blsys->system, blsys->y[i]);
585            CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
586            ASC_FREE(varname);
587        }
588    
589        /* print derivatives */
590        for(i=0; i < blsys->n_y; ++i){
591            if(blsys->ydot[i]){
592                varname = var_make_name(blsys->system, blsys->ydot[i]);
593                CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yp,i));
594                ASC_FREE(varname);
595            }else{
596                varname = var_make_name(blsys->system, blsys->y[i]);
597                CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
598                ASC_FREE(varname);
599            }
600        }
601    
602        /* print step size */
603        CONSOLE_DEBUG("<c_j> = %f",c_j);
604    
605        for(i=0, relptr = enginedata->rellist;
606                i< enginedata->nrels && relptr != NULL;
607                ++i, ++relptr
608        ){
609            relname = rel_make_name(blsys->system, *relptr);
610            CONSOLE_DEBUG("%d: '%s'",i,relname);
611        }
612        
613        /* build up the dense jacobian matrix... */
614        status = 0;
615        for(i=0, relptr = enginedata->rellist;
616                i< enginedata->nrels && relptr != NULL;
617                ++i, ++relptr
618        ){
619            relname = rel_make_name(blsys->system, *relptr);
620            CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
621            ASC_FREE(relname);
622    
623            /* get derivatives for this particular relation */
624            status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
625            /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
626    
627            for(j=0;j<count;++j){
628                varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
629                /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */
630                ASC_FREE(varname);
631            }
632    
633            if(status){
634                relname = rel_make_name(blsys->system, *relptr);
635                CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
636                ASC_FREE(relname);
637                break;
638            }
639    
640            for(j=0; j < enginedata->nrels; ++j){
641                DENSE_ELEM(Jac,i,j) = 0;
642            }
643    
644            for(j=0; j < count; ++j){
645                /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
646                varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
647                if(varname){
648                    CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
649                    ASC_FREE(varname);
650                }else{
651                    CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
652                }
653                
654                var_yindex = blsys->y_id[variables[j] - 1];
655                CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex);
656    
657                if(var_yindex >= 0){
658                    DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
659                    CONSOLE_DEBUG("New value (%d,%d) is %f", i, var_yindex, DENSE_ELEM(Jac,i,var_yindex));
660                }else{
661                    DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
662                    CONSOLE_DEBUG("New value (%d,%d) is %f (deriv)", i, -var_yindex-1, DENSE_ELEM(Jac,i,-var_yindex-1));
663                }
664            }
665        }
666    
667        CONSOLE_DEBUG("PRINTING JAC");
668        for(i=0; i < blsys->n_y; ++i){
669            for(j=0; j < blsys->n_y; ++j){
670                fprintf(stderr,"%8.2e\t",DENSE_ELEM(Jac,i,j));
671            }
672            fprintf(stderr,"\n");
673        }
674    
675    #if 1
676        double *yval = NV_DATA_S(yy);
677      /* AWFUL HACK! In an attempt to prove that 'IDADense' works as expected,      /* AWFUL HACK! In an attempt to prove that 'IDADense' works as expected,
678      I've pulled the jacobians straight from idadenx.c (an IDA example) */      I've pulled the jacobians straight from idadenx.c (an IDA example) */
679  # define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1)  #define A(REL,I,J,DESIRED) if(fabs((DENSE_ELEM(Jac,I,J) - DESIRED)/DESIRED) > 1e-5){ \
680      IJth(Jac,1,1) = RCONST(-0.04) - c_j;          CONSOLE_DEBUG("Value for ('%s'=%d,%d) got %e, expected %e",REL,I,J,DENSE_ELEM(Jac,I,J),DESIRED); \
681      IJth(Jac,2,1) = RCONST(0.04);          status=1; \
682      IJth(Jac,3,1) = 1.0;      }else{ \
683      IJth(Jac,1,2) = RCONST(1.0e4)*yval[2];          CONSOLE_DEBUG("Value for (%d,%d) is OK",I,J); \
684      IJth(Jac,2,2) = RCONST(-1.0e4)*yval[2] - RCONST(6.0e7)*yval[1] - c_j;      }
685      IJth(Jac,3,2) = 1.0;      char *eqs[] = {"eq1","eq2","eq3"};
686      IJth(Jac,1,3) = RCONST(1.0e4)*yval[1];      for(i=0;i<enginedata->nrels;++i){
687      IJth(Jac,2,3) = RCONST(-1.0e4)*yval[1];          relname = rel_make_name(blsys->system, enginedata->rellist[i]);
688      IJth(Jac,3,3) = 1.0;          if(strcmp(relname,"eq1")==0){
689  #undef IJth              A(relname,i,0, -0.04 - c_j   );
690      return(0);              A(relname,i,1, 1.0e4*yval[2] );
691                A(relname,i,2, 1.0e4*yval[1] );
692            }else if(strcmp(relname,"eq2")==0){
693                A(relname,i,0, 0.04);
694                A(relname,i,1, RCONST(-1.0e4)*yval[2] - RCONST(6.0e7)*yval[1] - c_j);
695                A(relname,i,2, RCONST(-1.0e4)*yval[1]);
696            }else if(strcmp(relname,"eq3")==0){
697                A(relname,i,0, -1.0);
698                A(relname,i,1, -1.0);
699                A(relname,i,2, -1.0);
700            }
701        }  
702    #undef A
703    #endif      
704    
705        if(status){
706            ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
707            return 1;
708        }
709    
710        return 0;
711  }  }
712    
713  /**  /**
# Line 663  int integrator_ida_jvex(realtype tt, N_V Line 788  int integrator_ida_jvex(realtype tt, N_V
788              /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */              /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
789    
790              for(j=0;j<count;++j){              for(j=0;j<count;++j){
791                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);                  /* varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
792                  /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */                  CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);
793                  ASC_FREE(varname);                  ASC_FREE(varname); */
794              }              }
795    
796              if(!status){              if(!status){
# Line 683  int integrator_ida_jvex(realtype tt, N_V Line 808  int integrator_ida_jvex(realtype tt, N_V
808    
809              Jv_i = 0;              Jv_i = 0;
810              for(j=0; j < count; ++j){              for(j=0; j < count; ++j){
811                  /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */                  /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
812                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]); */
813                  if(varname){                  if(varname){
814                      /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]); */                      /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
815                      ASC_FREE(varname);                      ASC_FREE(varname); */
816                  }else{                  }else{
817                      CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);                      CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
818                  }                  }
819                                    
820                  var_yindex = blsys->y_id[variables[j]];                  var_yindex = blsys->y_id[variables[j] - 1];
821                  /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */                  /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
822    
823                  if(var_yindex >= 0){                  if(var_yindex >= 0){
# Line 707  int integrator_ida_jvex(realtype tt, N_V Line 832  int integrator_ida_jvex(realtype tt, N_V
832    
833              NV_Ith_S(Jv,i) = Jv_i;              NV_Ith_S(Jv,i) = Jv_i;
834              /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */              /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */
835            }
836              if(status){          if(status){
837                  /* presumably some error_reporter will already have been made*/              CONSOLE_DEBUG("JVEX ERROR RESULT");
838                  is_error = 1;              /* presumably some error_reporter will already have been made*/
839              }              is_error = 1;
840          }          }
841      }else{      }else{
842          CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);          CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);

Legend:
Removed from v.945  
changed lines
  Added in v.946

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22