| 1 |
#include "idaanalyse.h" |
| 2 |
|
| 3 |
#include <utilities/ascPanic.h> |
| 4 |
|
| 5 |
#ifdef ASC_IDA_NEW_ANALYSE |
| 6 |
# include <solver/diffvars.h> |
| 7 |
# include <solver/slv_stdcalls.h> |
| 8 |
# include <solver/block.h> |
| 9 |
# include <solver/diffvars_impl.h> |
| 10 |
#include <solver/slvDOF.h> |
| 11 |
#endif |
| 12 |
|
| 13 |
#define ANALYSE_DEBUG |
| 14 |
|
| 15 |
#define VARMSG(MSG) \ |
| 16 |
varname = var_make_name(sys->system,v); \ |
| 17 |
CONSOLE_DEBUG(MSG,varname); \ |
| 18 |
ASC_FREE(varname) |
| 19 |
|
| 20 |
static int integrator_ida_check_partitioning(IntegratorSystem *sys); |
| 21 |
static int integrator_ida_check_diffindex(IntegratorSystem *sys); |
| 22 |
/* static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys); */ |
| 23 |
|
| 24 |
/** |
| 25 |
A var can be non-incident. If it *is* non incident and we're going to |
| 26 |
keep it, it will have to have derivative that *is* incident, and that |
| 27 |
meets the following filter. |
| 28 |
|
| 29 |
If it doesn't have a valid derivative (eg the derivative is fixed, or |
| 30 |
the variable doesn't HAVE a derivative), we will mark the non-deriv |
| 31 |
var non-ACTIVE, so anyway it will end up meeting this filter after we've |
| 32 |
run integrator_ida_check_vars. |
| 33 |
*/ |
| 34 |
const var_filter_t integrator_ida_nonderiv = { |
| 35 |
VAR_SVAR | VAR_ACTIVE | VAR_FIXED | VAR_DERIV, |
| 36 |
VAR_SVAR | VAR_ACTIVE | 0 | 0 |
| 37 |
}; |
| 38 |
|
| 39 |
const var_filter_t integrator_ida_deriv = { |
| 40 |
VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV, |
| 41 |
VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0 | VAR_DERIV |
| 42 |
}; |
| 43 |
|
| 44 |
const rel_filter_t integrator_ida_rel = { |
| 45 |
REL_INCLUDED | REL_EQUALITY | REL_ACTIVE, |
| 46 |
REL_INCLUDED | REL_EQUALITY | REL_ACTIVE |
| 47 |
}; |
| 48 |
|
| 49 |
/** |
| 50 |
Check derivative chains: if a var or derivative is inelligible, |
| 51 |
follow down the chain fixing & setting zero any derivatives. |
| 52 |
*/ |
| 53 |
static int integrator_ida_check_vars(IntegratorSystem *sys){ |
| 54 |
const SolverDiffVarCollection *diffvars; |
| 55 |
char *varname; |
| 56 |
int n_y = 0; |
| 57 |
int i, j; |
| 58 |
struct var_variable *v; |
| 59 |
SolverDiffVarSequence seq; |
| 60 |
int vok; |
| 61 |
|
| 62 |
CONSOLE_DEBUG("BEFORE CHECKING VARS"); |
| 63 |
integrator_ida_analyse_debug(sys,stderr); |
| 64 |
|
| 65 |
/* we shouldn't have allocated these yet: just be sure */ |
| 66 |
asc_assert(sys->y==NULL); |
| 67 |
asc_assert(sys->ydot==NULL); |
| 68 |
asc_assert(sys->n_y==0); |
| 69 |
|
| 70 |
/* get the the dervative chains from the system */ |
| 71 |
diffvars = system_get_diffvars(sys->system); |
| 72 |
if(diffvars==NULL){ |
| 73 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty"); |
| 74 |
return 1; |
| 75 |
} |
| 76 |
|
| 77 |
system_var_list_debug(sys->system); |
| 78 |
|
| 79 |
/* add the variables from the derivative chains */ |
| 80 |
for(i=0; i<diffvars->nseqs; ++i){ |
| 81 |
seq = diffvars->seqs[i]; |
| 82 |
asc_assert(seq.n >= 1); |
| 83 |
v = seq.vars[0]; |
| 84 |
|
| 85 |
vok = var_apply_filter(v,&integrator_ida_nonderiv); |
| 86 |
|
| 87 |
if(vok && !var_incident(v)){ |
| 88 |
VARMSG("good var '%s' is not incident"); |
| 89 |
/* var meets our filter, but perhaps it's not incident? */ |
| 90 |
if(seq.n == 1 || var_apply_filter(seq.vars[1],&integrator_ida_nonderiv)){ |
| 91 |
VARMSG("DEACTIVATING NON-INCIDENT VAR '%s' (NO DERIVATIVE)"); |
| 92 |
var_set_active(v,0); |
| 93 |
vok = 0; |
| 94 |
}else{ |
| 95 |
VARMSG("'%s' has a derivative that's OK"); |
| 96 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Non-incident var with an incident derivative. ASCEND can't handle this case at the moment, but we hope to fix it."); |
| 97 |
return 1; |
| 98 |
} |
| 99 |
} |
| 100 |
|
| 101 |
if(!vok){ |
| 102 |
VARMSG("'%s' fails non-deriv filter"); |
| 103 |
for(j=1;j<seq.n;++j){ |
| 104 |
v = seq.vars[j]; |
| 105 |
var_set_active(v,FALSE); |
| 106 |
var_set_value(v,0); |
| 107 |
VARMSG("Derivative '%s' SET ZERO AND INACTIVE"); |
| 108 |
} |
| 109 |
continue; |
| 110 |
} |
| 111 |
|
| 112 |
VARMSG("We will use var '%s'"); |
| 113 |
if(seq.n > 2){ |
| 114 |
varname = var_make_name(sys->system,v); |
| 115 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'",varname); |
| 116 |
ASC_FREE(varname); |
| 117 |
return 2; |
| 118 |
}else if(seq.n==2){ |
| 119 |
/* diff var */ |
| 120 |
if(var_apply_filter(seq.vars[1],&integrator_ida_deriv)){ |
| 121 |
/* add the diff & deriv vars to the lists */ |
| 122 |
n_y++; |
| 123 |
VARMSG("Added diff var '%s'"); |
| 124 |
v = seq.vars[1]; VARMSG("and its derivative '%s'"); |
| 125 |
continue; |
| 126 |
} |
| 127 |
VARMSG("Diff var '%s' being converted to alg var..."); |
| 128 |
v = seq.vars[1]; VARMSG("...because deriv var '%s' fails filter"); |
| 129 |
/* fall through */ |
| 130 |
} |
| 131 |
|
| 132 |
VARMSG("Adding '%s' to algebraic"); |
| 133 |
n_y++; |
| 134 |
} |
| 135 |
|
| 136 |
/* we assert that all vars in y meet the integrator_ida_nonderiv filter */ |
| 137 |
/* we assert that all vars in ydot meet the integrator_ida_deriv filter */ |
| 138 |
|
| 139 |
CONSOLE_DEBUG("Found %d good non-derivative vars", n_y); |
| 140 |
sys->n_y = n_y; |
| 141 |
|
| 142 |
return 0; |
| 143 |
} |
| 144 |
|
| 145 |
|
| 146 |
/** |
| 147 |
Sort the lists. First we will put the non-derivative vars, then we will |
| 148 |
put the derivative vars. Then we will put all the others. |
| 149 |
*/ |
| 150 |
static int integrator_ida_sort_rels_and_vars(IntegratorSystem *sys){ |
| 151 |
int ny1, nydot, nr; |
| 152 |
|
| 153 |
|
| 154 |
CONSOLE_DEBUG("BEFORE SORTING RELS AND VARS"); |
| 155 |
integrator_ida_analyse_debug(sys,stderr); |
| 156 |
|
| 157 |
/* we should not have allocated y or ydot yet */ |
| 158 |
asc_assert(sys->y==NULL && sys->ydot==NULL); |
| 159 |
|
| 160 |
/* but we should have found some variables (and know how many) */ |
| 161 |
asc_assert(sys->n_y); |
| 162 |
|
| 163 |
if(system_cut_vars(sys->system, 0, &integrator_ida_nonderiv, &ny1)){ |
| 164 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting non-derivs"); |
| 165 |
return 1; |
| 166 |
} |
| 167 |
|
| 168 |
CONSOLE_DEBUG("cut_vars: ny1 = %d, sys->n_y = %d",ny1,sys->n_y); |
| 169 |
asc_assert(ny1 == sys->n_y); |
| 170 |
|
| 171 |
if(system_cut_vars(sys->system, ny1, &integrator_ida_deriv, &nydot)){ |
| 172 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs"); |
| 173 |
return 1; |
| 174 |
} |
| 175 |
|
| 176 |
if(system_cut_rels(sys->system, 0, &integrator_ida_rel, &nr)){ |
| 177 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs"); |
| 178 |
return 1; |
| 179 |
} |
| 180 |
|
| 181 |
if(ny1 != nr){ |
| 182 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem is not square (ny = %d, nr = %d)",ny1,nr); |
| 183 |
return 2; |
| 184 |
} |
| 185 |
|
| 186 |
return 0; |
| 187 |
} |
| 188 |
|
| 189 |
/** |
| 190 |
Build up the lists 'ydot' and 'y_id. |
| 191 |
'ydot' allows us to find the derivative of a variable given its sindex. |
| 192 |
'y_id' allows us to locate the variable of a derivative given its sindex. |
| 193 |
|
| 194 |
Hence |
| 195 |
y_id[var_sindex(derivvar)-sys->n_y] = var_sindex(diffvar) |
| 196 |
ydot[var_sindex(diffvar)] = derivvar (or NULL if diffvar is algebraic) |
| 197 |
|
| 198 |
Note that there is 'shifting' required when looking up y_id. |
| 199 |
|
| 200 |
Note that we want to get rid of 'y' but for the moment we are also assigning |
| 201 |
that. We want to be rid of it because if the vars are ordered correctly |
| 202 |
it shouldn't be needed. |
| 203 |
*/ |
| 204 |
static int integrator_ida_create_lists(IntegratorSystem *sys){ |
| 205 |
const SolverDiffVarCollection *diffvars; |
| 206 |
int i, j; |
| 207 |
struct var_variable *v; |
| 208 |
char *varname; |
| 209 |
|
| 210 |
SolverDiffVarSequence seq; |
| 211 |
|
| 212 |
asc_assert(sys->y==NULL); |
| 213 |
asc_assert(sys->ydot==NULL); |
| 214 |
asc_assert(sys->y_id== NULL); |
| 215 |
|
| 216 |
/* get the the dervative chains from the system */ |
| 217 |
diffvars = system_get_diffvars(sys->system); |
| 218 |
if(diffvars==NULL){ |
| 219 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty"); |
| 220 |
return 1; |
| 221 |
} |
| 222 |
|
| 223 |
/* allocate space (more than enough) */ |
| 224 |
sys->y = ASC_NEW_ARRAY(struct var_variable *,sys->n_y); |
| 225 |
sys->y_id = ASC_NEW_ARRAY_CLEAR(int,sys->n_y); |
| 226 |
sys->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,sys->n_y); |
| 227 |
j = 0; |
| 228 |
|
| 229 |
CONSOLE_DEBUG("Passing through chains..."); |
| 230 |
|
| 231 |
/* add the variables from the derivative chains */ |
| 232 |
for(i=0; i<diffvars->nseqs; ++i){ |
| 233 |
|
| 234 |
CONSOLE_DEBUG("i = %d",i); |
| 235 |
|
| 236 |
seq = diffvars->seqs[i]; |
| 237 |
asc_assert(seq.n >= 1); |
| 238 |
v = seq.vars[0]; |
| 239 |
|
| 240 |
varname = var_make_name(sys->system, v); |
| 241 |
CONSOLE_DEBUG("alg '%s'",varname); |
| 242 |
ASC_FREE(varname); |
| 243 |
|
| 244 |
|
| 245 |
if(!var_apply_filter(v,&integrator_ida_nonderiv)){ |
| 246 |
continue; |
| 247 |
} |
| 248 |
|
| 249 |
varname = var_make_name(sys->system, v); |
| 250 |
CONSOLE_DEBUG("alg '%s' is GOOD",varname); |
| 251 |
ASC_FREE(varname); |
| 252 |
|
| 253 |
if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){ |
| 254 |
asc_assert(var_sindex(seq.vars[1]) >= sys->n_y); |
| 255 |
asc_assert(var_sindex(seq.vars[1])-sys->n_y < sys->n_y); |
| 256 |
|
| 257 |
varname = var_make_name(sys->system, seq.vars[1]); |
| 258 |
CONSOLE_DEBUG("diff '%s' IS GOOD",varname); |
| 259 |
ASC_FREE(varname); |
| 260 |
|
| 261 |
sys->y_id[var_sindex(seq.vars[1]) - sys->n_y] = j; |
| 262 |
sys->ydot[j] = seq.vars[1]; |
| 263 |
}else{ |
| 264 |
asc_assert(sys->ydot[j]==NULL); |
| 265 |
} |
| 266 |
|
| 267 |
sys->y[j] = v; |
| 268 |
j++; |
| 269 |
} |
| 270 |
|
| 271 |
return 0; |
| 272 |
} |
| 273 |
|
| 274 |
/** |
| 275 |
Perform additional problem analysis to prepare problem for integration with |
| 276 |
IDA. |
| 277 |
|
| 278 |
We assume that the analyse_generate_diffvars routine has been called, so |
| 279 |
that we just need to call slv_get_diffvars to access the derivative |
| 280 |
chains. |
| 281 |
|
| 282 |
We can also assume that the independent variable has been found. |
| 283 |
|
| 284 |
See mailing list, ~Jan 2007. |
| 285 |
|
| 286 |
Note, the stuff for identifying the static and output sub-problems should |
| 287 |
be part of integrator.c, not this file. We will assume this is handled |
| 288 |
|
| 289 |
@return 0 on success |
| 290 |
@see integrator_analyse |
| 291 |
*/ |
| 292 |
int integrator_ida_analyse(IntegratorSystem *sys){ |
| 293 |
int res; |
| 294 |
const SolverDiffVarCollection *diffvars; |
| 295 |
char *varname; |
| 296 |
int i; |
| 297 |
|
| 298 |
asc_assert(sys->engine==INTEG_IDA); |
| 299 |
|
| 300 |
#ifdef ANALYSE_DEBUG |
| 301 |
CONSOLE_DEBUG("Starting IDA analysis"); |
| 302 |
#endif |
| 303 |
|
| 304 |
res = integrator_ida_check_vars(sys); |
| 305 |
if(res){ |
| 306 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars"); |
| 307 |
return 1; |
| 308 |
} |
| 309 |
|
| 310 |
#ifdef ANALYSE_DEBUG |
| 311 |
CONSOLE_DEBUG("Sorting rels and vars"); |
| 312 |
#endif |
| 313 |
|
| 314 |
res = integrator_ida_sort_rels_and_vars(sys); |
| 315 |
if(res){ |
| 316 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars"); |
| 317 |
return 1; |
| 318 |
} |
| 319 |
|
| 320 |
#ifdef ANALYSE_DEBUG |
| 321 |
CONSOLE_DEBUG("Creating lists"); |
| 322 |
#endif |
| 323 |
|
| 324 |
CONSOLE_DEBUG("BEFORE MAKING LISTS"); |
| 325 |
integrator_ida_debug(sys,stderr); |
| 326 |
|
| 327 |
res = integrator_ida_create_lists(sys); |
| 328 |
if(res){ |
| 329 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists"); |
| 330 |
return 1; |
| 331 |
} |
| 332 |
|
| 333 |
#ifdef ANALYSE_DEBUG |
| 334 |
CONSOLE_DEBUG("Checking lists"); |
| 335 |
#endif |
| 336 |
|
| 337 |
asc_assert(sys->y); |
| 338 |
asc_assert(sys->ydot); |
| 339 |
asc_assert(sys->y_id); |
| 340 |
|
| 341 |
integrator_ida_debug(sys,stderr); |
| 342 |
|
| 343 |
if(integrator_ida_check_diffindex(sys)){ |
| 344 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex"); |
| 345 |
return 360; |
| 346 |
} |
| 347 |
|
| 348 |
/* the following causes TestIDA.testincidence3 to fail: |
| 349 |
if(sys->n_y != slv_get_num_solvers_rels(sys->system)){ |
| 350 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square"); |
| 351 |
return 2; |
| 352 |
}*/ |
| 353 |
|
| 354 |
#if 0 |
| 355 |
/* check structural singularity for the two IDACalcIC scenarios */ |
| 356 |
|
| 357 |
/* ...(1) FIX the derivatives */ |
| 358 |
CONSOLE_DEBUG("Checking system with derivatives fixed..."); |
| 359 |
for(i=0;i<n_y;++i){ |
| 360 |
if(sys->ydot[i])var_set_fixed(sys->ydot[i],1); |
| 361 |
} |
| 362 |
|
| 363 |
slv_block_partition(sys->system); |
| 364 |
res = integrator_ida_block_check(sys); |
| 365 |
if(res)return 100 + res; |
| 366 |
|
| 367 |
/* ...(2) FREE the derivatives, FIX the diffvars */ |
| 368 |
CONSOLE_DEBUG("Checking system with differential variables fixed..."); |
| 369 |
for(i=0;i<n_y;++i){ |
| 370 |
if(sys->ydot[i]){ |
| 371 |
var_set_fixed(sys->ydot[i],0); |
| 372 |
var_set_fixed(sys->y[i],1); |
| 373 |
} |
| 374 |
} |
| 375 |
slv_block_partition(sys->system); |
| 376 |
res = integrator_ida_block_check(sys); |
| 377 |
if(res)return 200 + res; |
| 378 |
|
| 379 |
/* ...(3) restore as it was, FREE the diffvars */ |
| 380 |
for(i=0;i<n_y;++i){ |
| 381 |
if(sys->ydot[i]){ |
| 382 |
var_set_fixed(sys->y[i],0); |
| 383 |
} |
| 384 |
} |
| 385 |
#endif |
| 386 |
|
| 387 |
/* get the the dervative chains from the system */ |
| 388 |
diffvars = system_get_diffvars(sys->system); |
| 389 |
|
| 390 |
/* check the indep var is same as was located elsewhere */ |
| 391 |
if(diffvars->nindep>1){ |
| 392 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)"); |
| 393 |
return 3; |
| 394 |
}else if(diffvars->nindep<1){ |
| 395 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)"); |
| 396 |
return 4; |
| 397 |
}else if(diffvars->indep[0]!=sys->x){ |
| 398 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match"); |
| 399 |
return 5; |
| 400 |
} |
| 401 |
|
| 402 |
#ifdef ANALYSE_DEBUG |
| 403 |
CONSOLE_DEBUG("Collecting observed variables"); |
| 404 |
#endif |
| 405 |
|
| 406 |
/* get the observations */ |
| 407 |
/** @TODO this should use a 'problem provider' API instead of static data |
| 408 |
from the system object */ |
| 409 |
sys->n_obs = diffvars->nobs; |
| 410 |
sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs); |
| 411 |
for(i=0;i<sys->n_obs;++i){ |
| 412 |
sys->obs[i] = diffvars->obs[i]; |
| 413 |
varname = var_make_name(sys->system,sys->obs[i]); |
| 414 |
CONSOLE_DEBUG("'%s' is observation",varname); |
| 415 |
ASC_FREE(varname); |
| 416 |
} |
| 417 |
|
| 418 |
return 0; |
| 419 |
} |
| 420 |
|
| 421 |
/*------------------------------------------------------------------------------ |
| 422 |
ANALYSIS ROUTINE (new implementation) |
| 423 |
*/ |
| 424 |
|
| 425 |
static int integrator_ida_check_partitioning(IntegratorSystem *sys){ |
| 426 |
long i, nv; |
| 427 |
int err=0; |
| 428 |
char *varname; |
| 429 |
struct var_variable **vlist, *v; |
| 430 |
nv = slv_get_num_solvers_vars(sys->system); |
| 431 |
vlist = slv_get_solvers_var_list(sys->system); |
| 432 |
var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED |
| 433 |
,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 }; |
| 434 |
for(i=0;i<nv;++i){ |
| 435 |
v=vlist[i]; |
| 436 |
if(!var_apply_filter(v,&vf))continue; |
| 437 |
varname = var_make_name(sys->system,v); |
| 438 |
if(!var_deriv(v)){ |
| 439 |
fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname); |
| 440 |
if(i>=sys->n_y){ |
| 441 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname); |
| 442 |
err++; |
| 443 |
} |
| 444 |
}else{ |
| 445 |
fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname); |
| 446 |
if(i<sys->n_y){ |
| 447 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)" |
| 448 |
,varname, sys->n_y, i |
| 449 |
); |
| 450 |
err++; |
| 451 |
} |
| 452 |
} |
| 453 |
ASC_FREE(varname); |
| 454 |
} |
| 455 |
if(err){ |
| 456 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning"); |
| 457 |
return 1; |
| 458 |
} |
| 459 |
return 0; |
| 460 |
} |
| 461 |
|
| 462 |
int integrator_ida_block_check(IntegratorSystem *sys){ |
| 463 |
int res; |
| 464 |
int dof; |
| 465 |
#ifdef ANALYSE_DEBUG |
| 466 |
int *vlist, *vp, i, nv, nv_ok; |
| 467 |
char *varname; |
| 468 |
struct var_variable **solversvars; |
| 469 |
|
| 470 |
var_filter_t vfilt = { |
| 471 |
VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED, |
| 472 |
VAR_ACTIVE | VAR_INCIDENT | 0 |
| 473 |
}; |
| 474 |
|
| 475 |
nv = slv_get_num_solvers_vars(sys->system); |
| 476 |
solversvars = slv_get_solvers_var_list(sys->system); |
| 477 |
CONSOLE_DEBUG("-------------- nv = %d -------------",nv); |
| 478 |
for(nv_ok=0, i=0; i < nv; ++i){ |
| 479 |
if(var_apply_filter(solversvars[i],&vfilt)){ |
| 480 |
varname = var_make_name(sys->system,solversvars[i]); |
| 481 |
fprintf(stderr,"%s\n",varname); |
| 482 |
ASC_FREE(varname); |
| 483 |
nv_ok++; |
| 484 |
} |
| 485 |
} |
| 486 |
CONSOLE_DEBUG("----------- got %d ok -------------",nv_ok); |
| 487 |
#endif |
| 488 |
|
| 489 |
if(!slvDOF_status(sys->system, &res, &dof)){ |
| 490 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status"); |
| 491 |
return -1; |
| 492 |
} |
| 493 |
switch(res){ |
| 494 |
case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break; |
| 495 |
case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */ |
| 496 |
case 3: CONSOLE_DEBUG("System is structurally singular"); break; |
| 497 |
case 4: CONSOLE_DEBUG("System is overspecified"); break; |
| 498 |
default: |
| 499 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status"); |
| 500 |
return -2; |
| 501 |
} |
| 502 |
|
| 503 |
#ifdef ANALYSE_DEBUG |
| 504 |
/* if it was underspecified, what vars could be fixed? */ |
| 505 |
if(res==1){ |
| 506 |
CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof); |
| 507 |
solversvars = slv_get_solvers_var_list(sys->system); |
| 508 |
if(!slvDOF_eligible(sys->system, &vlist)){ |
| 509 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list"); |
| 510 |
return -3; |
| 511 |
} |
| 512 |
for(vp=vlist;*vp!=-1;++vp){ |
| 513 |
varname = var_make_name(sys->system, solversvars[*vp]); |
| 514 |
CONSOLE_DEBUG("Fixable var: %s",varname); |
| 515 |
ASC_FREE(varname); |
| 516 |
} |
| 517 |
CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist)); |
| 518 |
return 1; |
| 519 |
} |
| 520 |
#endif |
| 521 |
|
| 522 |
return res; |
| 523 |
} |
| 524 |
|
| 525 |
/* return 0 on succes */ |
| 526 |
static int check_dups(struct var_variable **list,int n,int allownull){ |
| 527 |
int i,j; |
| 528 |
struct var_variable *v; |
| 529 |
for(i=0; i< n; ++i){ |
| 530 |
v=list[i]; |
| 531 |
if(v==NULL){ |
| 532 |
if(allownull)continue; |
| 533 |
else return 2; |
| 534 |
} |
| 535 |
for(j=0; j<i-1;++j){ |
| 536 |
if(list[j]==NULL)continue; |
| 537 |
if(v==list[j])return 1; |
| 538 |
} |
| 539 |
} |
| 540 |
return 0; |
| 541 |
} |
| 542 |
|
| 543 |
/* |
| 544 |
We are going to check that |
| 545 |
|
| 546 |
- n_y in range (0,n_vars) |
| 547 |
- no duplicates anywhere in the varlist |
| 548 |
- no duplicates in non-NULL elements of ydot |
| 549 |
|
| 550 |
- first n_y vars in solver's var list match those in vector y and match the non-deriv filter. |
| 551 |
- var_sindex for first n_y vars match their position the the solver's var list |
| 552 |
|
| 553 |
- ydot contains n_ydot non-NULL elements, each match the deriv filter. |
| 554 |
- vars in solver's var list positions [n_y,n_y+n_ydot) all match deriv filter |
| 555 |
|
| 556 |
- y_id contains n_ydot elements, with int values in range [0,n_y) |
| 557 |
- var_sindex(ydot[y_id[i]]) - n_y == i for i in [0,n_ydot) |
| 558 |
|
| 559 |
- all the vars [n_y+n_ydot,n_vars) fail the non-deriv filter and fail the deriv filter. |
| 560 |
*/ |
| 561 |
static int integrator_ida_check_diffindex(IntegratorSystem *sys){ |
| 562 |
int i, n_vars, n_ydot=0; |
| 563 |
struct var_variable **list, *v; |
| 564 |
char *varname; |
| 565 |
int diffindex; |
| 566 |
const char *msg; |
| 567 |
|
| 568 |
CONSOLE_DEBUG("Checking diffindex vector"); |
| 569 |
|
| 570 |
if(sys->y_id == NULL || sys->y == NULL || sys->ydot == NULL){ |
| 571 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"list(s) NULL"); |
| 572 |
return 1; |
| 573 |
} |
| 574 |
|
| 575 |
list = slv_get_solvers_var_list(sys->system); |
| 576 |
n_vars = slv_get_num_solvers_vars(sys->system); |
| 577 |
|
| 578 |
if(check_dups(list,n_vars,FALSE)){ |
| 579 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates or NULLs in solver's var list"); |
| 580 |
return 1; |
| 581 |
} |
| 582 |
|
| 583 |
if(check_dups(sys->ydot,n_vars,TRUE)){ |
| 584 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates in ydot vector"); |
| 585 |
return 1; |
| 586 |
} |
| 587 |
|
| 588 |
/* check n_y in range */ |
| 589 |
if(sys->n_y <= 0 || sys->n_y >= n_vars){ |
| 590 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"n_y = %d invalid (n_vars = %d)",sys->n_y, n_vars); |
| 591 |
return 1; |
| 592 |
} |
| 593 |
|
| 594 |
/* check first n_y vars */ |
| 595 |
for(i=0; i < sys->n_y; ++i){ |
| 596 |
v = list[i]; |
| 597 |
if(!var_apply_filter(v, &integrator_ida_nonderiv)){ |
| 598 |
msg = "'%s' (in first n_y vars) fails non-deriv filter"; goto finish; |
| 599 |
}else if(v != sys->y[i]){ |
| 600 |
msg = "'%s' not matched in y vector"; goto finish; |
| 601 |
}else if(var_sindex(v) != i){ |
| 602 |
msg = "'%s' has wrong var_sindex"; goto finish; |
| 603 |
} |
| 604 |
/* meets filter, matches in y vector, has correct var_sindex. */ |
| 605 |
} |
| 606 |
|
| 607 |
/* count non-NULLs in ydot */ |
| 608 |
for(i=0; i < sys->n_y; ++i){ |
| 609 |
v = sys->ydot[i]; |
| 610 |
if(v!=NULL){ |
| 611 |
if(var_sindex(v) < sys->n_y){ |
| 612 |
msg = "'%s' has var_sindex < n_y"; goto finish; |
| 613 |
}else if(!var_apply_filter(v,&integrator_ida_deriv)){ |
| 614 |
msg = "'%s' (in next n_ydot vars) fails deriv filter"; goto finish; |
| 615 |
} |
| 616 |
/* lies beyond first n_y vars, match deriv filter */ |
| 617 |
n_ydot++; |
| 618 |
} |
| 619 |
} |
| 620 |
|
| 621 |
/* check that vars [n_y, n_y+n_ydot) in solver's var list match the deriv filter */ |
| 622 |
for(i=sys->n_y; i< sys->n_y + n_ydot; ++i){ |
| 623 |
v = list[i]; |
| 624 |
if(!var_apply_filter(v,&integrator_ida_deriv)){ |
| 625 |
msg = "'%s', in range [n_y,n_y+n_ydot), fails deriv filter"; goto finish; |
| 626 |
} |
| 627 |
} |
| 628 |
|
| 629 |
/* check values in y_id are ints int range [0,n_y), and that they point to correct vars */ |
| 630 |
for(i=0; i<n_ydot; ++i){ |
| 631 |
if(sys->y_id[i] < 0 || sys->y_id[i] >= sys->n_y){ |
| 632 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"value of y_id[%d] is out of range",i); |
| 633 |
return 1; |
| 634 |
} |
| 635 |
v = sys->ydot[sys->y_id[i]]; |
| 636 |
if(var_sindex(v) - sys->n_y != i){ |
| 637 |
msg = "'%s' not index correctly from y_id"; goto finish; |
| 638 |
} |
| 639 |
} |
| 640 |
|
| 641 |
/* check remaining vars fail both filters */ |
| 642 |
for(i=sys->n_y + n_ydot; i<n_vars; ++i){ |
| 643 |
v = list[i]; |
| 644 |
if(var_apply_filter(v,&integrator_ida_nonderiv)){ |
| 645 |
msg = "Var '%s' at end meets non-deriv filter, but shouldn't"; goto finish; |
| 646 |
} |
| 647 |
if(var_apply_filter(v,&integrator_ida_deriv)){ |
| 648 |
CONSOLE_DEBUG("position = %d",i); |
| 649 |
msg = "Var '%s' at end meets deriv filter, but shouldn't"; goto finish; |
| 650 |
} |
| 651 |
} |
| 652 |
|
| 653 |
return 0; |
| 654 |
finish: |
| 655 |
varname = var_make_name(sys->system,v); |
| 656 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,msg,varname); |
| 657 |
ASC_FREE(varname); |
| 658 |
return 1; |
| 659 |
} |
| 660 |
|
| 661 |
/** |
| 662 |
@TODO provide a macro implementation of this |
| 663 |
*/ |
| 664 |
int integrator_ida_diffindex(const IntegratorSystem *sys, const struct var_variable *deriv){ |
| 665 |
int diffindex; |
| 666 |
#ifdef DIFFINDEX_DEBUG |
| 667 |
asc_assert( var_deriv (deriv)); |
| 668 |
asc_assert( var_active (deriv)); |
| 669 |
asc_assert( var_incident (deriv)); |
| 670 |
asc_assert( var_svar (deriv)); |
| 671 |
|
| 672 |
asc_assert(!var_fixed (deriv)); |
| 673 |
|
| 674 |
asc_assert(var_sindex(deriv) >= sys->n_y); |
| 675 |
asc_assert(diffindex == var_sindex(sys->y[diffindex])); |
| 676 |
#endif |
| 677 |
asc_assert(var_sindex(deriv) >= sys->n_y); |
| 678 |
diffindex = sys->y_id[var_sindex(deriv) - sys->n_y]; |
| 679 |
|
| 680 |
return diffindex; |
| 681 |
} |
| 682 |
|
| 683 |
|
| 684 |
/** |
| 685 |
This function will output the data structures provided to use BY THE |
| 686 |
SYSTEM -- not the ones we're working with here IN THE SOLVER. |
| 687 |
*/ |
| 688 |
int integrator_ida_analyse_debug(const IntegratorSystem *sys,FILE *fp){ |
| 689 |
return system_diffvars_debug(sys->system,fp); |
| 690 |
} |
| 691 |
|
Parent Directory
| 
Revision Log