Created
February 4, 2016 16:26
-
-
Save ehermes/04ee9aef38015f37da40 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| module solve_ida | |
| implicit none | |
| integer*8 :: neq = 4 | |
| integer*8 :: iout(25) | |
| real*8 :: rout(10) | |
| real*8 :: y0(4), yp0(4) | |
| real*8 :: diff(4), mas(4, 4) | |
| real*8 :: jac(4, 4) | |
| end module solve_ida | |
| program mkm | |
| implicit none | |
| integer*8, parameter :: nt = 100 | |
| integer*8, parameter :: nrates = 3 | |
| integer*8, parameter :: neq = 4 | |
| real*8 :: y(neq) | |
| real*8 :: u1(neq, nt), du1(neq, nt), r1(nrates, nt), t1(nt) | |
| t1 = 0.d0 | |
| u1 = 0.d0 | |
| du1 = 0.d0 | |
| r1 = 0.d0 | |
| y = (/2.d0, 0.d0, 0.9d0, 0.1d0/) | |
| print *, 'about to call initialize' | |
| call initialize(neq, y, 1d-9, (/1.d-11, 1.d-11, 1.d-11, 1.d-11/), (/0/), (/0.d0/), (/1.d0, 1.d0, 1.d0, 0.d0/)) | |
| print *, 'about to call solve' | |
| call solve(neq, nrates, nt, 36000.d0, t1, u1, du1, r1) | |
| print *, 'finished solving' | |
| print *, u1(:, nt) | |
| print *, du1(:, nt) | |
| print *, r1(:, nt) | |
| end program mkm | |
| subroutine initialize(neqin, y0in, rtol, atol, ipar, rpar, id_vec) | |
| use solve_ida, only: neq, iout, rout, y0, yp0, mas, diff | |
| implicit none | |
| integer*8, intent(in) :: neqin, ipar(*) | |
| real*8, intent(in) :: y0in(neqin), rtol, atol(*) | |
| real*8, intent(in) :: rpar(*) | |
| real*8, intent(in) :: id_vec(neqin) | |
| real*8 :: constr_vec(neqin) | |
| real*8 :: t0, yptmp(neqin) | |
| integer*8 :: nthreads, iatol, ier | |
| integer*8 :: i | |
| integer, external :: OMP_GET_NUM_THREADS | |
| iatol = 2 | |
| constr_vec = 1.d0 | |
| nthreads = OMP_GET_NUM_THREADS() | |
| y0 = y0in | |
| yp0 = 0 | |
| yptmp = 0 | |
| diff = id_vec | |
| mas = 0 | |
| t0 = 0 | |
| do i = 1, neq | |
| mas(i, i) = id_vec(i) | |
| enddo | |
| ! Calculate yp | |
| call fidaresfun(0.d0, y0, yptmp, yp0, ipar, rpar, ier) | |
| ! initialize Sundials | |
| call fnvinits(2, neq, ier) | |
| !call fnvinitomp(2, neq, nthreads, ier) | |
| ! allocate memory | |
| call fidamalloc(t0, y0, yp0, iatol, rtol, atol, iout, & | |
| rout, ipar, rpar, ier) | |
| ! set maximum number of steps | |
| call fidasetiin('MAX_NSTEPS', 5000000, ier) | |
| ! set algebraic variables | |
| call fidasetvin('ID_VEC', id_vec, ier) | |
| ! set constraints (all yi >= 0.) | |
| call fidasetvin('CONSTR_VEC', constr_vec, ier) | |
| ! ! initialize the solver | |
| ! call fidalapackdense(neq, ier) | |
| ! ! enable the jacobian | |
| ! call fidalapackdensesetjac(1, ier) | |
| ! initialize the solver | |
| call fidaspgmr(0, 0, 0, 0, 0, ier) | |
| ! call fidaspbcg(0, 0, 0, ier) | |
| ! call fidasptfqmr(0, 0, 0, ier) | |
| ! enable the jacobian | |
| ! call fidaspilssetjac(1, ier) | |
| ! enable the preconditioner | |
| call fidaspilssetprec(1, ier) | |
| end subroutine initialize | |
| subroutine solve(neqin, nrates, nt, tfinal, t1, u1, du1, r1) | |
| use solve_ida, only: y0, yp0 | |
| implicit none | |
| integer*8, intent(in) :: neqin, nt, nrates | |
| real*8, intent(in) :: tfinal | |
| real*8 :: yp(neqin) | |
| real*8 :: rpar(1) | |
| integer*8 :: ipar(1) | |
| real*8, intent(out) :: t1(nt) | |
| real*8, intent(out) :: u1(neqin, nt), du1(neqin, nt) | |
| real*8, intent(out) :: r1(nrates, nt) | |
| real*8 :: dt, tout | |
| integer*8 :: itask, ier | |
| integer*8 :: i | |
| yp = 0.d0 | |
| itask = 1 | |
| dt = tfinal / (nt - 1) | |
| tout = 0.0d0 | |
| u1 = 0 | |
| du1 = 0 | |
| u1(:, 1) = y0 | |
| du1(:, 1) = yp0 | |
| t1(1) = 0.d0 | |
| call ratecalc(4, 3, u1(:, 1), r1(:, 1)) | |
| do i = 2, nt | |
| tout = tout + dt | |
| call fidasolve(tout, t1(i), u1(:, i), du1(:, i), itask, ier) | |
| ! call fidaresfun(t1(i), u1(:, i), yp, du1(:, i), ipar, rpar, ier) | |
| call ratecalc(4, 3, u1(:, i), r1(:, i)) | |
| end do | |
| end subroutine solve | |
| subroutine finalize | |
| ! use solve_ida, only: y0, yp0, jac, mas, diff | |
| implicit none | |
| call fidafree | |
| end subroutine finalize | |
| subroutine fidaresfun(tres, yin, ypin, res, ipar, rpar, reserr) | |
| use solve_ida, only: neq, diff | |
| implicit none | |
| integer*8, intent(in) :: ipar(*) | |
| integer*8, intent(out) :: reserr | |
| real*8, intent(in) :: tres, rpar(*) | |
| real*8, intent(in) :: yin(neq), ypin(neq) | |
| real*8, intent(out) :: res(neq) | |
| real*8 :: x(1) | |
| res = 0 | |
| res(1) = -3.12034517121819d-10*yin(1)*yin(2)*exp(35.1649924759196d0 & | |
| *yin(1)) - 1.71569483215224d-16*yin(1)* & | |
| yin(4)*exp(35.1649924759196d0*yin(1)) + & | |
| 243626660.931473d0*yin(3) + 243658599.808889d0* & | |
| yin(4) | |
| res(2) = -3.12034517121819d-10*yin(1)*yin(2)*exp(35.1649924759196d0 & | |
| *yin(1)) + 243658599.808889d0*yin(4) | |
| res(3) = 1.71569483215224d-16*yin(1)*yin(4)*exp(35.1649924759196d0* & | |
| yin(1)) - 1.86861519121121d-7*yin(3)**2 - & | |
| 243626660.931473d0*yin(3) + 1729692473943.15d0* & | |
| yin(4)**2 | |
| res(4) = -yin(2) - yin(3) - yin(4) + 1 | |
| res = res - diff * ypin | |
| reserr = 0 | |
| end subroutine fidaresfun | |
| subroutine fidadjac(neqin, t, yin, ypin, r, jac, cj, ewt, h, & | |
| ipar, rpar, wk1, wk2, wk3, djacerr) | |
| use solve_ida, only: mas | |
| implicit none | |
| integer*8 :: neqin, ipar(*) | |
| integer*8 :: djacerr | |
| real*8 :: t, h, cj, rpar(*) | |
| real*8 :: yin(neqin), ypin(neqin), r(neqin), ewt(*), jac(neqin, neqin) | |
| real*8 :: wk1(*), wk2(*), wk3(*) | |
| jac = 0 | |
| jac(1, 1) = -1.0972691446816d-8*yin(1)*yin(2)*exp(35.1649924759196d0* & | |
| yin(1)) - 6.03323958636075d-15*yin(1)*yin(4) & | |
| *exp(35.1649924759196d0*yin(1)) - 3.12034517121819d-10* & | |
| yin(2)*exp(35.1649924759196d0*yin(1)) - & | |
| 1.71569483215224d-16*yin(4)*exp(35.1649924759196d0* & | |
| yin(1)) | |
| jac(2, 1) = -1.0972691446816d-8*yin(1)*yin(2)*exp(35.1649924759196d0* & | |
| yin(1)) - 3.12034517121819d-10*yin(2)*exp( & | |
| 35.1649924759196d0*yin(1)) | |
| jac(3, 1) = 6.03323958636075d-15*yin(1)*yin(4)*exp(35.1649924759196d0* & | |
| yin(1)) + 1.71569483215224d-16*yin(4)*exp( & | |
| 35.1649924759196d0*yin(1)) | |
| jac(1, 2) = -3.12034517121819d-10*yin(1)*exp(35.1649924759196d0*yin(1) & | |
| ) | |
| jac(2, 2) = -3.12034517121819d-10*yin(1)*exp(35.1649924759196d0*yin(1) & | |
| ) | |
| jac(4, 2) = -1 | |
| jac(1, 3) = 243626660.931473d0 | |
| jac(3, 3) = -3.73723038242241d-7*yin(3) - 243626660.931473d0 | |
| jac(4, 3) = -1 | |
| jac(1, 4) = -1.71569483215224d-16*yin(1)*exp(35.1649924759196d0*yin(1) & | |
| ) + 243658599.808889d0 | |
| jac(2, 4) = 243658599.808889d0 | |
| jac(3, 4) = 1.71569483215224d-16*yin(1)*exp(35.1649924759196d0*yin(1) & | |
| ) + 3459384947886.3d0*yin(4) | |
| jac(4, 4) = -1 | |
| jac = jac - cj * mas | |
| djacerr = 0 | |
| end subroutine fidadjac | |
| subroutine ratecalc(neqin, nrates, yin, rates) | |
| implicit none | |
| integer, intent(in) :: neqin, nrates | |
| real*8, intent(in) :: yin(neqin) | |
| real*8, intent(out) :: rates(nrates) | |
| rates = 0 | |
| rates(1) = -3.12034517121819d-10*yin(1)*yin(2)*exp(35.1649924759196d0 & | |
| *yin(1)) + 243658599.808889d0*yin(4) | |
| rates(2) = -9.34307595605603d-8*yin(3)**2 + 864846236971.575d0*yin(4) & | |
| **2 | |
| rates(3) = -1.71569483215224d-16*yin(1)*yin(4)*exp(35.1649924759196d0 & | |
| *yin(1)) + 243626660.931473d0*yin(3) | |
| end subroutine ratecalc | |
| subroutine fidajtimes(tres, yin, ypin, res, vin, fjv, cj, ewt, h, ipar, rpar, wk1, wk2, ier) | |
| use solve_ida, only: neq | |
| implicit none | |
| real*8, intent(in) :: tres, yin(neq), ypin(neq), res(neq), vin(neq), cj, h | |
| real*8 :: ewt(*), wk1(*), wk2(*), rpar(*) | |
| integer*8 :: ipar(*) | |
| integer :: i | |
| real*8, intent(out) :: fjv(neq) | |
| integer*8, intent(out) :: ier | |
| real*8 :: jac(neq, neq) | |
| call fidadjac(neq, tres, yin, ypin, res, jac, cj, ewt, h, & | |
| ipar, rpar, wk1, wk2, wk2, ier) | |
| do i = 1, neq | |
| fjv(i) = dot_product(vin, jac(i, :)) | |
| enddo | |
| end subroutine fidajtimes | |
| subroutine fidapsol(tres, yin, ypin, res, rvin, zv, cj, delta, ewt, ipar, rpar, wk1, ier) | |
| use solve_ida, only: neq, jac | |
| implicit none | |
| real*8, intent(in) :: tres, yin(neq), ypin(neq), res(neq), rvin(neq) | |
| real*8, intent(in) :: cj, delta, ewt(*), rpar(*) | |
| integer*8, intent(in) :: ipar(*) | |
| real*8 :: wk1(*) | |
| integer*8 :: ier | |
| real*8, intent(out) :: zv(neq) | |
| real*8 :: jaclu(neq, neq) | |
| integer :: ipiv(neq) | |
| integer :: i | |
| ! call fidadjac(neq, tres, yin, ypin, res, jaclu, cj, ewt, 1, & | |
| ! ipar, rpar, wk1, wk1, wk1, ier) | |
| zv = rvin | |
| jaclu = jac | |
| call dgesv(neq, 1, jaclu, neq, ipiv, zv, neq, ier) | |
| end subroutine fidapsol | |
| subroutine fidapset(tres, yin, ypin, res, cj, ewt, h, ipar, rpar, wk1, wk2, wk3, ier) | |
| use solve_ida, only: neq, jac | |
| implicit none | |
| real*8, intent(in) :: tres, yin(neq), ypin(neq), res(neq) | |
| real*8, intent(in) :: cj, ewt(*), h, rpar(*) | |
| real*8 :: wk1(*), wk2(*), wk3(*) | |
| integer*8, intent(in) :: ipar(*) | |
| integer*8, intent(out) :: ier | |
| call fidadjac(neq, tres, yin, ypin, res, jac, cj, ewt, h, & | |
| ipar, rpar, wk1, wk2, wk3, ier) | |
| end subroutine fidapset | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment