| 1 |
/* ASCEND modelling environment |
| 2 |
Copyright (C) 2006-2011 Carnegie Mellon University |
| 3 |
|
| 4 |
This program is free software; you can redistribute it and/or modify |
| 5 |
it under the terms of the GNU General Public License as published by |
| 6 |
the Free Software Foundation; either version 2, or (at your option) |
| 7 |
any later version. |
| 8 |
|
| 9 |
This program is distributed in the hope that it will be useful, |
| 10 |
but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 11 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 12 |
GNU General Public License for more details. |
| 13 |
|
| 14 |
You should have received a copy of the GNU General Public License |
| 15 |
along with this program; if not, write to the Free Software |
| 16 |
Foundation, Inc., 59 Temple Place - Suite 330, |
| 17 |
Boston, MA 02111-1307, USA. |
| 18 |
*//** |
| 19 |
@file |
| 20 |
Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes |
| 21 |
as part of the GPL-licensed SUNDIALS solver package from LLNL. |
| 22 |
|
| 23 |
IDA provides the non-linear parts, as well as a number of pluggable linear |
| 24 |
solvers: dense, banded and krylov types. |
| 25 |
|
| 26 |
We also implement here an EXPERIMENTAL direct sparse linear solver for IDA |
| 27 |
using the ASCEND linsolqr routines. |
| 28 |
|
| 29 |
@see http://www.llnl.gov/casc/sundials/ |
| 30 |
*//* |
| 31 |
by John Pye, May 2006 |
| 32 |
*/ |
| 33 |
|
| 34 |
#define _GNU_SOURCE |
| 35 |
|
| 36 |
#include <signal.h> |
| 37 |
#include <setjmp.h> |
| 38 |
#include <fenv.h> |
| 39 |
#include <math.h> |
| 40 |
|
| 41 |
/* SUNDIALS includes */ |
| 42 |
#ifdef ASC_WITH_IDA |
| 43 |
|
| 44 |
#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2 |
| 45 |
# include <sundials/sundials_config.h> |
| 46 |
# include <sundials/sundials_nvector.h> |
| 47 |
# include <ida/ida_spgmr.h> |
| 48 |
# include <ida.h> |
| 49 |
# include <nvector_serial.h> |
| 50 |
#else |
| 51 |
# include <sundials/sundials_config.h> |
| 52 |
# include <nvector/nvector_serial.h> |
| 53 |
# include <ida/ida.h> |
| 54 |
#endif |
| 55 |
|
| 56 |
# include <sundials/sundials_dense.h> |
| 57 |
# include <ida/ida_spgmr.h> |
| 58 |
# include <ida/ida_spbcgs.h> |
| 59 |
# include <ida/ida_sptfqmr.h> |
| 60 |
# include <ida/ida_dense.h> |
| 61 |
|
| 62 |
# ifndef IDA_SUCCESS |
| 63 |
# error "Failed to include SUNDIALS IDA header file" |
| 64 |
# endif |
| 65 |
#else |
| 66 |
# error "If you're building this file, you should have ASC_WITH_IDA" |
| 67 |
#endif |
| 68 |
|
| 69 |
#ifdef ASC_WITH_MMIO |
| 70 |
# include <mmio.h> |
| 71 |
#endif |
| 72 |
|
| 73 |
#include <ascend/general/platform.h> |
| 74 |
#include <ascend/utilities/error.h> |
| 75 |
#include <ascend/utilities/ascSignal.h> |
| 76 |
#include <ascend/general/panic.h> |
| 77 |
#include <ascend/compiler/instance_enum.h> |
| 78 |
|
| 79 |
#include <ascend/system/slv_client.h> |
| 80 |
#include <ascend/system/relman.h> |
| 81 |
#include <ascend/system/block.h> |
| 82 |
#include <ascend/system/slv_stdcalls.h> |
| 83 |
#include <ascend/system/jacobian.h> |
| 84 |
#include <ascend/system/bndman.h> |
| 85 |
|
| 86 |
#include <ascend/utilities/config.h> |
| 87 |
#include <ascend/integrator/integrator.h> |
| 88 |
|
| 89 |
#include "idalinear.h" |
| 90 |
#include "idaanalyse.h" |
| 91 |
#include "ida_impl.h" |
| 92 |
#include "idaprec.h" |
| 93 |
/* |
| 94 |
for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2 |
| 95 |
*/ |
| 96 |
#ifndef SUNDIALS_VERSION_MINOR |
| 97 |
# ifdef __GNUC__ |
| 98 |
# warning "GUESSING SUNDIALS VERSION 2.2" |
| 99 |
# endif |
| 100 |
# define SUNDIALS_VERSION_MINOR 2 |
| 101 |
#endif |
| 102 |
#ifndef SUNDIALS_VERSION_MAJOR |
| 103 |
# define SUNDIALS_VERSION_MAJOR 2 |
| 104 |
#endif |
| 105 |
|
| 106 |
/* SUNDIALS 2.4.0 introduces new DlsMat in place of DenseMat */ |
| 107 |
#if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==4 |
| 108 |
# define IDA_MTX_T DlsMat |
| 109 |
# define IDADENSE_SUCCESS IDADLS_SUCCESS |
| 110 |
# define IDADENSE_MEM_NULL IDADLS_MEM_NULL |
| 111 |
# define IDADENSE_ILL_INPUT IDADLS_ILL_INPUT |
| 112 |
# define IDADENSE_MEM_FAIL IDADLS_MEM_FAIL |
| 113 |
#else |
| 114 |
# define IDA_MTX_T DenseMat |
| 115 |
#endif |
| 116 |
|
| 117 |
/* #define FEX_DEBUG */ |
| 118 |
#define JEX_DEBUG |
| 119 |
/* #define DJEX_DEBUG */ |
| 120 |
#define SOLVE_DEBUG |
| 121 |
#define STATS_DEBUG |
| 122 |
#define PREC_DEBUG |
| 123 |
/* #define ROOT_DEBUG */ |
| 124 |
|
| 125 |
/* #define DIFFINDEX_DEBUG */ |
| 126 |
/* #define ANALYSE_DEBUG */ |
| 127 |
/* #define DESTROY_DEBUG */ |
| 128 |
/* #define MATRIX_DEBUG */ |
| 129 |
|
| 130 |
|
| 131 |
/*------ |
| 132 |
Full jacobian preconditioner -- experimental |
| 133 |
*/ |
| 134 |
|
| 135 |
static int integrator_ida_psetup_jacobian(realtype tt, |
| 136 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 137 |
realtype c_j, void *prec_data, |
| 138 |
N_Vector tmp1, N_Vector tmp2, |
| 139 |
N_Vector tmp3 |
| 140 |
); |
| 141 |
|
| 142 |
static int integrator_ida_psolve_jacobian(realtype tt, |
| 143 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 144 |
N_Vector rvec, N_Vector zvec, |
| 145 |
realtype c_j, realtype delta, void *prec_data, |
| 146 |
N_Vector tmp |
| 147 |
); |
| 148 |
|
| 149 |
static void integrator_ida_pcreate_jacobian(IntegratorSystem *integ); |
| 150 |
|
| 151 |
const IntegratorIdaPrec prec_jacobian = { |
| 152 |
integrator_ida_pcreate_jacobian |
| 153 |
, integrator_ida_psetup_jacobian |
| 154 |
, integrator_ida_psolve_jacobian |
| 155 |
}; |
| 156 |
|
| 157 |
/*------ |
| 158 |
Jacobi preconditioner -- experimental |
| 159 |
*/ |
| 160 |
|
| 161 |
static int integrator_ida_psetup_jacobi(realtype tt, |
| 162 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 163 |
realtype c_j, void *prec_data, |
| 164 |
N_Vector tmp1, N_Vector tmp2, |
| 165 |
N_Vector tmp3 |
| 166 |
); |
| 167 |
|
| 168 |
static int integrator_ida_psolve_jacobi(realtype tt, |
| 169 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 170 |
N_Vector rvec, N_Vector zvec, |
| 171 |
realtype c_j, realtype delta, void *prec_data, |
| 172 |
N_Vector tmp |
| 173 |
); |
| 174 |
|
| 175 |
static void integrator_ida_pcreate_jacobi(IntegratorSystem *integ); |
| 176 |
|
| 177 |
const IntegratorIdaPrec prec_jacobi = { |
| 178 |
integrator_ida_pcreate_jacobi |
| 179 |
, integrator_ida_psetup_jacobi |
| 180 |
, integrator_ida_psolve_jacobi |
| 181 |
}; |
| 182 |
|
| 183 |
/*---------------------------------------------- |
| 184 |
FULL JACOBIAN PRECONDITIONER -- EXPERIMENTAL. |
| 185 |
*/ |
| 186 |
|
| 187 |
static void integrator_ida_pcreate_jacobian(IntegratorSystem *integ){ |
| 188 |
IntegratorIdaData *enginedata =integ->enginedata; |
| 189 |
IntegratorIdaPrecDataJacobian *precdata; |
| 190 |
precdata = ASC_NEW(IntegratorIdaPrecDataJacobian); |
| 191 |
mtx_matrix_t P; |
| 192 |
asc_assert(integ->n_y); |
| 193 |
precdata->L = linsolqr_create_default(); |
| 194 |
|
| 195 |
/* allocate matrix to be used by linsolqr */ |
| 196 |
P = mtx_create(); |
| 197 |
mtx_set_order(P, integ->n_y); |
| 198 |
linsolqr_set_matrix(precdata->L, P); |
| 199 |
|
| 200 |
enginedata->pfree = &integrator_ida_pfree_jacobian; |
| 201 |
enginedata->precdata = precdata; |
| 202 |
CONSOLE_DEBUG("Allocated memory for Full Jacobian preconditioner"); |
| 203 |
} |
| 204 |
|
| 205 |
void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata){ |
| 206 |
mtx_matrix_t P; |
| 207 |
IntegratorIdaPrecDataJacobian *precdata; |
| 208 |
|
| 209 |
if(enginedata->precdata){ |
| 210 |
precdata = (IntegratorIdaPrecDataJacobian *)enginedata->precdata; |
| 211 |
P = linsolqr_get_matrix(precdata->L); |
| 212 |
mtx_destroy(P); |
| 213 |
linsolqr_destroy(precdata->L); |
| 214 |
ASC_FREE(precdata); |
| 215 |
enginedata->precdata = NULL; |
| 216 |
|
| 217 |
CONSOLE_DEBUG("Freed memory for Full Jacobian preconditioner"); |
| 218 |
} |
| 219 |
enginedata->pfree = NULL; |
| 220 |
} |
| 221 |
|
| 222 |
/** |
| 223 |
EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers |
| 224 |
|
| 225 |
'setup' function. |
| 226 |
*/ |
| 227 |
static int integrator_ida_psetup_jacobian(realtype tt, |
| 228 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 229 |
realtype c_j, void *p_data, |
| 230 |
N_Vector tmp1, N_Vector tmp2, |
| 231 |
N_Vector tmp3 |
| 232 |
){ |
| 233 |
int i, j, res; |
| 234 |
IntegratorSystem *integ; |
| 235 |
IntegratorIdaData *enginedata; |
| 236 |
IntegratorIdaPrecDataJacobian *precdata; |
| 237 |
linsolqr_system_t L; |
| 238 |
mtx_matrix_t P; |
| 239 |
struct rel_relation **relptr; |
| 240 |
|
| 241 |
integ = (IntegratorSystem *)p_data; |
| 242 |
enginedata = integ->enginedata; |
| 243 |
precdata = (IntegratorIdaPrecDataJacobian *)(enginedata->precdata); |
| 244 |
double *derivatives; |
| 245 |
struct var_variable **variables; |
| 246 |
int count, status; |
| 247 |
char *relname; |
| 248 |
mtx_coord_t C; |
| 249 |
|
| 250 |
L = precdata->L; |
| 251 |
P = linsolqr_get_matrix(L); |
| 252 |
mtx_clear(P); |
| 253 |
|
| 254 |
CONSOLE_DEBUG("Setting up Jacobian preconditioner"); |
| 255 |
|
| 256 |
variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2); |
| 257 |
derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2); |
| 258 |
|
| 259 |
/** |
| 260 |
@TODO FIXME here we are using the very inefficient and contorted approach |
| 261 |
of calculating the whole jacobian, then extracting just the diagonal elements. |
| 262 |
*/ |
| 263 |
|
| 264 |
for(i=0, relptr = enginedata->rellist; |
| 265 |
i< enginedata->nrels && relptr != NULL; |
| 266 |
++i, ++relptr |
| 267 |
){ |
| 268 |
/* get derivatives for this particular relation */ |
| 269 |
status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval); |
| 270 |
if(status){ |
| 271 |
relname = rel_make_name(integ->system, *relptr); |
| 272 |
CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname); |
| 273 |
ASC_FREE(relname); |
| 274 |
break; |
| 275 |
} |
| 276 |
/* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */ |
| 277 |
/* find the diagonal elements */ |
| 278 |
for(j=0; j<count; ++j){ |
| 279 |
if(var_deriv(variables[j])){ |
| 280 |
mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), c_j * derivatives[j]); |
| 281 |
}else{ |
| 282 |
mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), derivatives[j]); |
| 283 |
} |
| 284 |
} |
| 285 |
} |
| 286 |
|
| 287 |
mtx_assemble(P); |
| 288 |
|
| 289 |
if(status){ |
| 290 |
CONSOLE_DEBUG("Error found when evaluating derivatives"); |
| 291 |
res = 1; goto finish; /* recoverable */ |
| 292 |
} |
| 293 |
|
| 294 |
integrator_ida_write_incidence(integ); |
| 295 |
|
| 296 |
res = 0; |
| 297 |
finish: |
| 298 |
ASC_FREE(variables); |
| 299 |
ASC_FREE(derivatives); |
| 300 |
return res; |
| 301 |
}; |
| 302 |
|
| 303 |
/** |
| 304 |
EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers |
| 305 |
|
| 306 |
'solve' function. |
| 307 |
*/ |
| 308 |
static int integrator_ida_psolve_jacobian(realtype tt, |
| 309 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 310 |
N_Vector rvec, N_Vector zvec, |
| 311 |
realtype c_j, realtype delta, void *p_data, |
| 312 |
N_Vector tmp |
| 313 |
){ |
| 314 |
IntegratorSystem *integ; |
| 315 |
IntegratorIdaData *data; |
| 316 |
IntegratorIdaPrecDataJacobian *precdata; |
| 317 |
integ = (IntegratorSystem *)p_data; |
| 318 |
data = integ->enginedata; |
| 319 |
precdata = (IntegratorIdaPrecDataJacobian *)(data->precdata); |
| 320 |
linsolqr_system_t L = precdata->L; |
| 321 |
|
| 322 |
linsolqr_add_rhs(L,NV_DATA_S(rvec),FALSE); |
| 323 |
|
| 324 |
mtx_region_t R; |
| 325 |
R.row.low = R.col.low = 0; |
| 326 |
R.row.high = R.col.high = mtx_order(linsolqr_get_matrix(L)) - 1; |
| 327 |
linsolqr_set_region(L,R); |
| 328 |
|
| 329 |
linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L))); |
| 330 |
linsolqr_reorder(L, &R, linsolqr_rmethod(L)); |
| 331 |
|
| 332 |
/// @TODO more here |
| 333 |
|
| 334 |
linsolqr_remove_rhs(L,NV_DATA_S(rvec)); |
| 335 |
|
| 336 |
CONSOLE_DEBUG("Solving Jacobian preconditioner (c_j = %f)",c_j); |
| 337 |
return 0; |
| 338 |
}; |
| 339 |
|
| 340 |
|
| 341 |
/*---------------------------------------------- |
| 342 |
JACOBI PRECONDITIONER -- EXPERIMENTAL. |
| 343 |
*/ |
| 344 |
|
| 345 |
static void integrator_ida_pcreate_jacobi(IntegratorSystem *integ){ |
| 346 |
IntegratorIdaData *enginedata =integ->enginedata; |
| 347 |
IntegratorIdaPrecDataJacobi *precdata; |
| 348 |
precdata = ASC_NEW(IntegratorIdaPrecDataJacobi); |
| 349 |
|
| 350 |
asc_assert(integ->n_y); |
| 351 |
precdata->PIii = N_VNew_Serial(integ->n_y); |
| 352 |
|
| 353 |
enginedata->pfree = &integrator_ida_pfree_jacobi; |
| 354 |
enginedata->precdata = precdata; |
| 355 |
CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner"); |
| 356 |
} |
| 357 |
|
| 358 |
void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){ |
| 359 |
if(enginedata->precdata){ |
| 360 |
IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata; |
| 361 |
N_VDestroy_Serial(precdata->PIii); |
| 362 |
|
| 363 |
ASC_FREE(precdata); |
| 364 |
enginedata->precdata = NULL; |
| 365 |
CONSOLE_DEBUG("Freed memory for Jacobi preconditioner"); |
| 366 |
} |
| 367 |
enginedata->pfree = NULL; |
| 368 |
} |
| 369 |
|
| 370 |
/** |
| 371 |
EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers |
| 372 |
|
| 373 |
'setup' function. |
| 374 |
*/ |
| 375 |
static int integrator_ida_psetup_jacobi(realtype tt, |
| 376 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 377 |
realtype c_j, void *p_data, |
| 378 |
N_Vector tmp1, N_Vector tmp2, |
| 379 |
N_Vector tmp3 |
| 380 |
){ |
| 381 |
int i, j, res; |
| 382 |
IntegratorSystem *integ; |
| 383 |
IntegratorIdaData *enginedata; |
| 384 |
IntegratorIdaPrecDataJacobi *precdata; |
| 385 |
struct rel_relation **relptr; |
| 386 |
|
| 387 |
integ = (IntegratorSystem *)p_data; |
| 388 |
enginedata = integ->enginedata; |
| 389 |
precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata); |
| 390 |
double *derivatives; |
| 391 |
struct var_variable **variables; |
| 392 |
int count, status; |
| 393 |
char *relname; |
| 394 |
|
| 395 |
CONSOLE_DEBUG("Setting up Jacobi preconditioner"); |
| 396 |
|
| 397 |
variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2); |
| 398 |
derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2); |
| 399 |
|
| 400 |
/** |
| 401 |
@TODO FIXME here we are using the very inefficient and contorted approach |
| 402 |
of calculating the whole jacobian, then extracting just the diagonal elements. |
| 403 |
*/ |
| 404 |
|
| 405 |
for(i=0, relptr = enginedata->rellist; |
| 406 |
i< enginedata->nrels && relptr != NULL; |
| 407 |
++i, ++relptr |
| 408 |
){ |
| 409 |
|
| 410 |
/* get derivatives for this particular relation */ |
| 411 |
status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval); |
| 412 |
if(status){ |
| 413 |
relname = rel_make_name(integ->system, *relptr); |
| 414 |
CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname); |
| 415 |
ASC_FREE(relname); |
| 416 |
break; |
| 417 |
} |
| 418 |
/* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */ |
| 419 |
/* find the diagonal elements */ |
| 420 |
for(j=0; j<count; ++j){ |
| 421 |
if(var_sindex(variables[j])==i){ |
| 422 |
if(var_deriv(variables[j])){ |
| 423 |
NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]); |
| 424 |
}else{ |
| 425 |
NV_Ith_S(precdata->PIii, i) = 1./derivatives[j]; |
| 426 |
} |
| 427 |
|
| 428 |
} |
| 429 |
} |
| 430 |
#ifdef PREC_DEBUG |
| 431 |
CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i)); |
| 432 |
#endif |
| 433 |
} |
| 434 |
|
| 435 |
if(status){ |
| 436 |
CONSOLE_DEBUG("Error found when evaluating derivatives"); |
| 437 |
res = 1; goto finish; /* recoverable */ |
| 438 |
} |
| 439 |
|
| 440 |
integrator_ida_write_incidence(integ); |
| 441 |
|
| 442 |
res = 0; |
| 443 |
finish: |
| 444 |
ASC_FREE(variables); |
| 445 |
ASC_FREE(derivatives); |
| 446 |
return res; |
| 447 |
}; |
| 448 |
|
| 449 |
/** |
| 450 |
EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers |
| 451 |
|
| 452 |
'solve' function. |
| 453 |
*/ |
| 454 |
static int integrator_ida_psolve_jacobi(realtype tt, |
| 455 |
N_Vector yy, N_Vector yp, N_Vector rr, |
| 456 |
N_Vector rvec, N_Vector zvec, |
| 457 |
realtype c_j, realtype delta, void *p_data, |
| 458 |
N_Vector tmp |
| 459 |
){ |
| 460 |
IntegratorSystem *integ; |
| 461 |
IntegratorIdaData *data; |
| 462 |
IntegratorIdaPrecDataJacobi *precdata; |
| 463 |
integ = (IntegratorSystem *)p_data; |
| 464 |
data = integ->enginedata; |
| 465 |
precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata); |
| 466 |
|
| 467 |
CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j); |
| 468 |
N_VProd(precdata->PIii, rvec, zvec); |
| 469 |
return 0; |
| 470 |
}; |
| 471 |
|
| 472 |
|
Parent Directory
| 
Revision Log