SloppyCell.daskr

99 """Integrate a system of ordinary differential algebraic equations with or 100 without events. 101 102 Description: 103 104 Solve a system of ordinary differential algebraic equations using the 105 FORTRAN library DDASKR. 106 107 Uses the backward differentiation formulas of orders one through five 108 to solve a system of the form G(T,Y,Y') = 0 where Y can be a vector. 109 110 Values for Y and Y' at the initial time must be given as input. These 111 values should be consistent, that is, if t0, y0, yp0 are the given 112 initial values, they should satisfy G(t0,y0,yp0) = 0. However, if 113 consistent values are not known, in many cases you can have DDASKR 114 solve for them if the appropriate option in the info array is set. 115 116 Normally, DDASKR solves the system from some intial t0 to a specified 117 tout. This wrapper takes a list of times t as input, and then solves 118 the system for the range of times. In the interval mode of operation, 119 only the solution values at time points in t will be returned. 120 Intermediate results can also be obtained easily by specifying info(3). 121 122 While integrating the given DAE system, DDASKR also searches for roots 123 of the given constraint functions Ri(T,Y,Y') given by RT. If DDASKR 124 detects a sign change in any Ri(T,Y,Y'), it will return the intermediate 125 value of T and Y for which Ri(T,Y,Y') = 0. 126 127 Caution: If some Ri has a root at or very near the initial time, DDASKR 128 may fail to find it, or may find extraneous roots there, because it does 129 not yet have a sufficient history of the solution. 130 131 Inputs: 132 133 res -- res(t, y, ...) computes the residual function for the 134 differential/algebraic system to be solved. 135 t -- a sequence of time points for which to solve for y. The intial 136 value of the dependent variable(s) y should correspond to the 137 first time in this sequence. 138 y0 -- this array must contain the initial values of the NEQ solution 139 components at the initial time point (note NEQ is defined as the 140 length of this array). 141 yp0 -- this array must contain the initial values of the NEQ first 142 derivatives of the solution components at the initial time point. 143 jac -- this is the name of a routine that you may supply that relates to 144 the Jacobian matrix of the nonlinear system that the code must 145 solve at each time step. 146 rtol -- array of relative error tolerances (must be non-negative). 147 atol -- array of absolute error tolerances (must be non-negative). 148 rt -- this is the name of the subroutine for defining the vector 149 Ri(T, Y, Y') of constraint functions. 150 nrt -- the number of constraint functions Ri(T, Y, Y'). 151 jac -- the name of a routine that you may supply (optionally) that 152 relactes to the Jacobian matrix of the nonlinear system that the 153 code must solve at each time step. If no Jacobian is passed then 154 the system will automatically use a numerical Jacobian. 155 args -- parameter for passing extra information to the wrapper. 156 tstop -- sometimes it is not possible to integrate past some point tstop 157 because the equation changes there or it is not defined past 158 tstop. if tstop is set the integrator will not integrate past 159 tstop. NOTE: If you select a tstop, the integrator may skip 160 some points in your list of time points, t. This is because it 161 may pass tstop before hitting all the points in the list of times 162 for an integration that is progressing quickly. 163 intermediate_output -- True to have all output returned. 164 ineq_constr -- True to have inequality constraint checking on. 165 If True, the code will check for non-negativity in Y during the 166 integration (not during initial condition calculations). 167 calculate_ic -- True to have the DDASKR automatically calculate 168 consistent initial conditions. Given Y_d, the code will 169 calculate Y_a and Y'_d (note: you must also pass var_types 170 to indicate which variables are algebraic and which are 171 differential). We do not use the functionality of daskr 172 that allows calculation of Y given Y'. 173 var_types -- this list tells whether each variable in the system is 174 differential or algebraic. Vor a variable Yi, var_types[i] = +1 if 175 Yi is differential, and var_types[i] = -1 if it is algebraic. 176 redir_output -- False to have print statements from DASKR printed to the 177 standard output. If redir_output is True, the output will be 178 redirected. 179 max_steps -- set this if you want the integrator to be able to take more 180 than 500 steps between time points when you're not in 181 "intermediate output" mode. If intermediate_output = True then 182 every time step is considered an output point so there are never 183 multiple steps between time points. If intermediate_output = 184 False (default) then the integrator will stop after taking 500 185 steps unless max_steps is set > 500. The default value for 186 max_steps is 500. 187 rpar -- If not None, this sequence is passed to the function being 188 evaluated and can be used to pass additional arguments. 189 hmax -- If not None, this sets an upper limit to the step sizes that 190 ddaskr will take during the integration. Setting this can help 191 integration for difficult systems, but it can also slow things 192 down. 193 194 Outputs: (yout, tout, ypout, t_root, y_root, i_root) 195 196 yout -- the numerically calculated list of solution points corresponding 197 to the times in tout. 198 tout -- the list of time points corresponding to yout. 199 ypout -- the numerically calculated list of solution points corresponding 200 to the times in tout. 201 t_root -- the time of a terminating event. 202 y_root -- the system state at time t_root. 203 i_root -- tells which event(s) fired and how they fired. 204 If i_root(i) = 0, then event i did not fire. If nonzero, 205 i_root(i) shows the direction of the sign change in Ri. 206 i_root(i) = 1 means Ri changed from negative to positive. 207 i_root(i) = -1 means Ri changed from positive to negative. 208 209 """ 210 211 if scipy.isscalar(rpar): 212 raise ValueError('rpar must be a sequence.') 213 214 if rpar is None or len(rpar) == 0: 215 # rpar needs to have at least length 1 or daskr will screw up arguments. 216 rpar = scipy.empty(1) 217 218 y = y0 219 yp = yp0 220 ires = 0 221 222 # We calculate the number of equations in the residual function from the 223 # length of the dependent variable array. 224 # The number of event functions are passed to daeint as nrt. 225 neq = len(y0) 226 # ipar is used to hold neq and the length of rpar. 227 ipar = [neq, len(rpar)] 228 229 # The size of the work arrays normally depends on the options we select. 230 # Since info(11)=0 always and info(8)=0 always in this wrapper, we 231 # can always set BASEr, BASEi, and MAXORD as follows: 232 MAXORD = 5 233 BASEr = 60+max(MAXORD+4,7)*neq+3*nrt 234 BASEi = 40 235 # Since info(5) = 0, add neq*2 to LRW 236 LRW = BASEr + (neq**2) 237 LIW = BASEi + 2*neq 238 239 # Make the work arrays 240 # Making these int32 and float64 avoids crashes on 64-bit. 241 rwork = scipy.zeros(LRW, scipy.float64) 242 iwork = scipy.zeros(LIW, scipy.int32) 243 244 # The info array is used to give the daskr code more details about how we 245 # want the problem solved. 246 # Set all the options for info, as well as related options. 247 info = scipy.zeros(20, scipy.int32) 248 249 # Now we go through the options one by one and set them according to the 250 # inputs. Note that some options are not used and kept at '0'. 251 252 # Is this the first call for the problem? (yes) 253 info[0] = 0 254 255 # Are both error tolerances rtol and atol scalars? (no, we use arrays of 256 # tolerances) 257 info[1] = 1 258 259 # Now we take care of the tolerance settings 260 if scipy.isscalar(rtol) or len(rtol) != neq: 261 raise ValueError('rtol must be an array or a vector with same length ' 262 'as number of equations.') 263 if scipy.isscalar(atol) or len(atol) != neq: 264 raise ValueError('atol must be an array or a vector with same length ' 265 'as number of equations.') 266 267 # Do you want the solution at intermediate steps (and not just at the 268 # specified points)? 269 if intermediate_output == False: 270 info[2] = 0 271 elif intermediate_output == True: 272 info[2] = 1 273 else: 274 raise ValueError('int_output must be True or False.') 275 276 # Can the integration be carried out without any restrictions on the 277 # independent variable t? 278 if tstop == None: 279 info[3] = 0 280 else: 281 # check to make sure that tstop is greater than or equal to the final 282 # time point requested 283 if tstop < t[-1]: 284 raise ValueError('tstop must be greater than or equal to the ' 285 'final time point requested.') 286 info[3] = 1 287 rwork[0] = tstop 288 # There appears to be a bug in ddaskr. If we're in intermediate output mode 289 # and we have events, the daskr may skip a requested time and instead 290 # return an event firing time. A reasonable fix appears to be using tstop 291 # to prevent this. 292 EVENT_INT_BUG_WORKAROUND = False 293 if intermediate_output and nrt: 294 EVENT_INT_BUG_WORKAROUND = True 295 296 # Do you want the code to evaluate the partial derivatives automatically by 297 # numerical differences? 298 if jac == None: 299 info[4] = 0 300 jac = dummy_func 301 else: 302 info[4] = 1 303 304 # info[5], info[6], info[7], info[8] are not used 305 if hmax is not None: 306 info[7] = 1 307 rwork[2] = hmax 308 309 # Do you want the code to solve the problem without invoking any special 310 # inequality constraints? 311 # LID is a variable used for setting up which variables are algebraic and 312 # which are differential. The value of LID depends on info(9). 313 # For now we only allow options 0 or 2 so that we don't have to worry 314 # about telling the integrator how each variable is constrained. 315 # Note: This wrapper doesn't don't support info[9] = 0 or 2 316 LID = 0 317 if ineq_constr == False: 318 info[9] = 0 319 LID = 40 320 #elif ineq_constr == 1: 321 # info[9] = 1 322 # LID = 40 + neq 323 elif ineq_constr == True: 324 info[9] = 2 325 LID = 40 326 #elif ineq_constr == 3: 327 # info[9] = 3 328 # LID = 40 + neq 329 else: 330 raise ValueError, 'ineq_constr must be True or False.' 331 332 # Is it necessary to calculate the initial condition? 333 # (i.e., are the initial conditions consistent?) 334 if calculate_ic == False: 335 info[10] = 0 336 # if they are inconsistent, given Y_d, calculate Y_a and Y'_d must also 337 # specify which components are algebraic and which are differential 338 elif calculate_ic == True: 339 info[10] = 1 340 # check the var_types parameter 341 if var_types == None: 342 raise ValueError('if calculate_ic = 1, the var_types array ' 343 'must be passed.') 344 elif len(var_types) != len(y0): 345 raise ValueError('var_types must be of length len(y0)') 346 else: 347 # +1 is for differential variables 348 # -1 is for algebraic variables 349 for ii in range(len(var_types)): 350 if var_types[ii] == +1: 351 iwork[LID+ii] = +1 352 elif var_types[ii] == -1: 353 iwork[LID+ii] = -1 354 else: 355 raise ValueError('the var_types array may only contain ' 356 'entries of +1 or -1.') 357 358 359 # info[11] and info[12] are not used 360 361 # Do you want to proceed to the integration after the initial condition 362 # calculation is done? 363 # (for now we will assume that we will always stop to capture the 364 # calculated initial condition) 365 # Note that once the initial condition is calculated, info[13] must 366 # be reset to 0 to allow the integration to proceed 367 info[13] = 1 368 369 # info[14] is not used 370 371 # Do you wish to control errors locally on all the variables? 372 # (for now we will assume the answer is yes) 373 info[15] = 0 374 375 # Are the intial condition heuristic controls to be given their default 376 # values? 377 # (for now we will assume yes) 378 info[16] = 0 379 380 #info[18] is not used. 381 382 # if no root function was passed in, we should assign the dummy function 383 if rt == None: 384 rt = dummy_func 385 386 # assign the dummy function to psol since we don't use it 387 psol = dummy_func 388 389 390 # Set up for first call 391 idid = 1 # First call 392 t0, tindex = t[0], 1 393 tout, yout = t0, y0 394 t_root, y_root, i_root = [], [], [] 395 396 # variables for storing the output. 397 tout_l = [] 398 yout_l = [] 399 ypout_l= [] 400 401 # add the initial point to the output array if the initial condition is 402 # consistent 403 if calculate_ic == False: 404 tout_l.append(t0) 405 yout_l.append(y0) 406 ypout_l.append(yp0) 407 408 # set the current time to the initial time 409 tcurrent = t0 410 411 # step_count is used to keep track of how many times idid = -1 occurs in a 412 # row 413 step_count = 0 414 415 if redir_output == True: 416 redir.start() 417 418 if hasattr(res, '_cpointer'): 419 res = res._cpointer 420 if hasattr(jac, '_cpointer'): 421 jac = jac._cpointer 422 if hasattr(rt, '_cpointer'): 423 rt = rt._cpointer 424 425 try: 426 while tcurrent < t[-1]: 427 # set the desired output time 428 twanted = t[tindex] 429 if EVENT_INT_BUG_WORKAROUND: 430 info[3] = 1 431 if tstop and tstop < twanted: 432 rwork[0] = tstop 433 else: 434 rwork[0] = twanted 435 436 # continue the integration 437 treached, y, yp, idid, jroot = \ 438 _daskr.ddaskr(res, tcurrent, y, yp, twanted, 439 info, rtol, atol, 440 rwork, iwork, 441 rpar, ipar, 442 jac, psol, rt, nrt) 443 444 # check for a negative value of idid so we know if there was a 445 # problem 446 if idid <= 0: 447 # set appropriate options 448 info[0]=1 449 # idid=-1 is handled below 450 if idid < -1: 451 messages = redir.stop() 452 # The task was interrupted 453 if messages is not None: 454 logger.warn(messages) 455 logger.warn(_msgs[idid]) 456 # if idid=-1, 500 steps have been taken sine the last time 457 # idid=-1. 458 # Whether we should continue depends on the value of max_steps. 459 elif idid == -1: 460 # if we haven't hit the value of max_steps, then we should 461 # continue the integration. Otherwise we should raise an 462 # exception. 463 if max_steps > step_count: 464 step_count += 500 465 continue 466 elif max_steps <= step_count: 467 messages = redir.stop() 468 # report the error message for idid = -1 469 if messages is not None: 470 logger.warn(messages) 471 logger.warn(_msgs[idid]) 472 473 474 # send what output was obtained 475 # If the integrator tried and failed to calculate the initial 476 # condition, then yout_l will be of length 0, and have the 477 # wrong shape. 478 # We replace it and ypout_l with arrays containing 0 rows, but 479 # the proper number of columns, for consistency with 480 # what we would get if it had succeeded. 481 if len(yout_l) == 0: 482 yout_l = scipy.zeros((0, len(y0)), scipy.float_) 483 ypout_l = scipy.zeros((0, len(y0)), scipy.float_) 484 outputs = (scipy.array(yout_l), scipy.array(tout_l), 485 scipy.array(ypout_l), 486 t_root, y_root, i_root) 487 raise daeintException(_msgs[idid], outputs) 488 489 # if there were no errors, continue to post output or take 490 # appropriate action if we triggered an event or tstop 491 else: 492 493 # any time a succesful idid occurs we should reset step_count 494 step_count = 0 495 496 # updating of output should be done whether or not we are in 497 # intermediate output mode 498 499 # update the output 500 tout_l.append(treached) 501 yout_l.append(y) 502 ypout_l.append(yp) 503 504 # update the current time 505 tcurrent = treached 506 507 # if we hit one of the specified times, update time index 508 if treached == t[tindex]: 509 # update the time index 510 tindex += 1 511 512 # if the time we reached is tstop, then stop the integration 513 if idid == 2: 514 break 515 516 # if an initial condition was calculated, update info[13] 517 # this allows the integration to proceed. Also set 518 # calculate_ic to False to avoid a redundant initial condition 519 # calculation 520 if idid == 4: 521 info[13] = 0 522 calculate_ic = False 523 524 # if the solution was successful because a root of R(T,Y,Y') was 525 # found at treached, then restart the integration. 526 elif idid == 5: 527 t_root = tcurrent 528 y_root = y 529 i_root = jroot 530 # include the following break if we want all events to be 531 # terminating 532 break 533 534 535 # format the output 536 tout_l = scipy.array(tout_l) 537 yout_l = scipy.array(yout_l) 538 ypout_l = scipy.array(ypout_l) 539 540 # Process outputs 541 outputs = (yout_l, tout_l, ypout_l, t_root, y_root, i_root) 542 543 return outputs 544 545 finally: 546 messages = redir.stop()

Source Code for Module SloppyCell.daskr