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 |
|
|
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; |
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){ |
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; |