/[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 927 by johnpye, Thu Nov 9 13:57:02 2006 UTC revision 928 by johnpye, Wed Nov 22 10:32:18 2006 UTC
# Line 174  int integrator_ida_solve( Line 174  int integrator_ida_solve(
174          return 0; /* failure */          return 0; /* failure */
175      }      }
176    
177        CONSOLE_DEBUG("RETRIEVING t0");
178    
179      /* retrieve initial values from the system */      /* retrieve initial values from the system */
180      t0 = samplelist_get(blsys->samples,start_index);      t0 = samplelist_get(blsys->samples,start_index);
181    
182        CONSOLE_DEBUG("RETRIEVING y0");
183    
184      y0 = N_VNew_Serial(size);      y0 = N_VNew_Serial(size);
185      integrator_get_y(blsys,NV_DATA_S(y0));      integrator_get_y(blsys,NV_DATA_S(y0));
186    
187        CONSOLE_DEBUG("RETRIEVING yp0");
188    
189      yp0 = N_VNew_Serial(size);      yp0 = N_VNew_Serial(size);
190      integrator_get_ydot(blsys,NV_DATA_S(yp0));      integrator_get_ydot(blsys,NV_DATA_S(yp0));
191    
# Line 426  int integrator_ida_jvex(realtype tt, N_V Line 432  int integrator_ida_jvex(realtype tt, N_V
432      char *relname, *varname;      char *relname, *varname;
433      int status;      int status;
434      double Jv_i;      double Jv_i;
435        int var_yindex;
436    
437      int *variables;      int *variables;
438      double *derivatives;      double *derivatives;
439      var_filter_t filter;      var_filter_t filter;
440      int count;      int count;
441    
442        fprintf(stderr,"\n--------------\n");
443      CONSOLE_DEBUG("EVALUTING JACOBIAN...");      CONSOLE_DEBUG("EVALUTING JACOBIAN...");
444    
445      blsys = (IntegratorSystem *)jac_data;      blsys = (IntegratorSystem *)jac_data;
# Line 467  int integrator_ida_jvex(realtype tt, N_V Line 475  int integrator_ida_jvex(realtype tt, N_V
475                  i< enginedata->nrels && relptr != NULL;                  i< enginedata->nrels && relptr != NULL;
476                  ++i, ++relptr                  ++i, ++relptr
477          ){          ){
478                fprintf(stderr,"\n");
479                relname = rel_make_name(blsys->system, *relptr);
480                CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
481                ASC_FREE(relname);
482    
483              /* get derivatives for this particular relation */              /* get derivatives for this particular relation */
484              status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);              status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
485                CONSOLE_DEBUG("Got derivatives against %d matching variables", count);
486    
487              for(j=0;j<count;++j){              for(j=0;j<count;++j){
488                  varname = var_make_name(blsys->system, enginedata->varlist[blsys->y_id[variables[i]]]);                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
489                  CONSOLE_DEBUG("diff var[%d] is %s",j,varname);                  CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);
490                  ASC_FREE(varname);                  ASC_FREE(varname);
491              }              }
492    
             CONSOLE_DEBUG("Got derivatives against %d matching variables", count);  
   
             relname = rel_make_name(blsys->system, *relptr);  
493              if(!status){              if(!status){
494                  CONSOLE_DEBUG("Derivatives for relation %d '%s' OK",i,relname);                  CONSOLE_DEBUG("Derivatives for relation %d OK",i);
495              }else{              }else{
496                  CONSOLE_DEBUG("ERROR calculating derivatives for relation %d '%s'",i,relname);                  CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
497                  break;                  break;
498              }              }
499              ASC_FREE(relname);  
500                /*
501                    Now we have the derivatives wrt each alg/diff variable in the
502                    present equation. variables[] points into the varlist. need
503                    a mapping from the varlist to the y and ydot lists.
504                */
505    
506              Jv_i = 0;              Jv_i = 0;
507              for(j=0; j < count; ++j){              for(j=0; j < count; ++j){
# Line 498  int integrator_ida_jvex(realtype tt, N_V Line 515  int integrator_ida_jvex(realtype tt, N_V
515                  }                  }
516                                    
517                  /* this bit is still not right!!! */                  /* this bit is still not right!!! */
518                  Jv_i += derivatives[j] * NV_Ith_S(v,blsys->y_id[variables[j]]);                  var_yindex = blsys->y_id[variables[j]];
519    
520                    if(var_yindex >= 0){
521                        CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));
522                        Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
523                    }else{
524                        var_yindex = -var_yindex-1;
525                        CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));
526                        Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
527                    }
528              }              }
529    
530              NV_Ith_S(Jv,i) = Jv_i;              NV_Ith_S(Jv,i) = Jv_i;

Legend:
Removed from v.927  
changed lines
  Added in v.928

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