Line data Source code
1 : !-----------------------------------------------------------------------
2 : !
3 : ! This file is part of the Test Set for IVP solvers
4 : ! http://www.dm.uniba.it/~testset/
5 : !
6 : ! Chemical Akzo Nobel problem
7 : ! ODE of dimension 6
8 : !
9 : ! DISCLAIMER: see
10 : ! http://www.dm.uniba.it/~testset/disclaimer.php
11 : !
12 : ! The most recent version of this source file can be found at
13 : ! http://www.dm.uniba.it/~testset/src/problems/chemakzo.f
14 : !
15 : ! This is revision
16 : ! $Id: chemakzo.F,v 1.2 2006/10/02 10:29:14 testset Exp $
17 : !
18 : !-----------------------------------------------------------------------
19 : module bari_chemakzo
20 :
21 : implicit none
22 :
23 : contains
24 :
25 0 : subroutine chemakzo_prob(fullnm, problm, type, neqn, ndisc, t, numjac, mljac, mujac, nummas, mlmas, mumas, ind)
26 : character(len=*) :: fullnm, problm, type
27 : integer :: neqn, ndisc, mljac, mujac, mlmas, mumas, ind(*)
28 : double precision :: t(0:*)
29 : logical :: numjac, nummas
30 :
31 : integer :: i
32 :
33 0 : fullnm = 'Chemical Akzo Nobel problem'
34 0 : problm = 'chemakzo'
35 0 : type = 'DAE'
36 0 : neqn = 6
37 0 : ndisc = 0
38 0 : t(0) = 0d0
39 0 : t(1) = 180d0
40 0 : numjac = .false.
41 0 : mljac = neqn
42 0 : mujac = neqn
43 0 : mlmas = 0
44 0 : mumas = 0
45 0 : do i = 1, neqn
46 0 : ind(i) = 0
47 : end do
48 :
49 0 : end subroutine chemakzo_prob
50 : !-----------------------------------------------------------------------
51 0 : subroutine chemakzo_init(neqn, y, yprime, consis)
52 : integer :: neqn
53 : double precision :: y(neqn), yprime(neqn)
54 : logical :: consis
55 :
56 : integer :: ierr, ipar(0)
57 : double precision :: rpar(0)
58 :
59 : double precision :: ks
60 : parameter(ks=115.83d0)
61 :
62 0 : y(1) = 0.444d0
63 0 : y(2) = 0.00123d0
64 0 : y(3) = 0d0
65 0 : y(4) = 0.007d0
66 0 : y(5) = 0d0
67 0 : y(6) = ks*y(1)*y(4)
68 :
69 0 : consis = .true.
70 :
71 0 : call chemakzo_feval(neqn, 0d0, y, y, yprime, ierr, rpar, ipar)
72 :
73 0 : end subroutine chemakzo_init
74 : !-----------------------------------------------------------------------
75 0 : subroutine chemakzo_feval(neqn, t, y, yprime, f, ierr, rpar, ipar)
76 : integer :: neqn, ierr, ipar(*)
77 : double precision :: t, y(neqn), yprime(neqn), f(neqn), rpar(*)
78 :
79 : double precision :: k1, k2, k3, k4, kbig, kla, po2, hen, ks
80 : parameter(k1=18.7d0, k2=0.58d0, k3=0.09d0, k4=0.42d0, kbig=34.4d0, kla=3.3d0, ks=115.83d0, po2=0.9d0, hen=737d0)
81 0 : double precision :: r1, r2, r3, r4, r5, fin, sqy2
82 : include 'formats'
83 :
84 0 : if (y(2) < 0d0) then
85 : !write(*,*) 'y(2)', y(2)
86 0 : ierr = -1
87 0 : return
88 : end if
89 :
90 0 : sqy2 = sqrt(y(2))
91 0 : r1 = k1*y(1)*y(1)*y(1)*y(1)*sqy2
92 0 : r2 = k2*y(3)*y(4)
93 0 : r3 = k2/kbig*y(1)*y(5)
94 0 : r4 = k3*y(1)*y(4)*y(4)
95 0 : r5 = k4*(y(6)**2)*sqy2
96 0 : fin = kla*(po2/hen - y(2))
97 :
98 0 : f(1) = ((-2d0*r1 + r2) - r3) - r4
99 0 : f(2) = ((-0.5d0*r1 - r4) - 0.5d0*r5) + fin
100 0 : f(3) = (r1 - r2) + r3
101 0 : f(4) = (-r2 + r3) - 2d0*r4
102 0 : f(5) = (r2 - r3) + r5
103 0 : f(6) = ks*y(1)*y(4) - y(6)
104 :
105 : end subroutine chemakzo_feval
106 : !-----------------------------------------------------------------------
107 0 : subroutine chemakzo_jeval(ldim, neqn, t, y, yprime, dfdy, ierr, rpar, ipar)
108 : integer :: ldim, neqn, ierr, ipar(*)
109 : double precision :: t, y(neqn), yprime(neqn), dfdy(ldim, neqn), rpar(*)
110 :
111 : double precision :: k1, k2, k3, k4, kbig, kla, ks
112 : parameter(k1=18.7d0, k2=0.58d0, k3=0.09d0, k4=0.42d0, kbig=34.4d0, kla=3.3d0, ks=115.83d0)
113 : integer :: i, j
114 0 : double precision :: r11, r12, r23, r24, r31, r35, r41, r44, r52, r56, fin2
115 0 : double precision :: y13
116 :
117 0 : if (y(2) < 0d0) then
118 0 : ierr = -1
119 0 : return
120 : end if
121 :
122 0 : do j = 1, neqn
123 0 : do i = 1, neqn
124 0 : dfdy(i, j) = 0d0
125 : end do
126 : end do
127 :
128 0 : y13 = y(1)*y(1)*y(1)
129 0 : r11 = 4d0*k1*y13*sqrt(y(2))
130 0 : r12 = 0.5d0*k1*y(1)*y13/sqrt(y(2))
131 0 : r23 = k2*y(4)
132 0 : r24 = k2*y(3)
133 0 : r31 = (k2/kbig)*y(5)
134 0 : r35 = (k2/kbig)*y(1)
135 0 : r41 = k3*y(4)**2
136 0 : r44 = 2d0*k3*y(1)*y(4)
137 0 : r52 = 0.5d0*k4*(y(6)**2)/sqrt(y(2))
138 0 : r56 = 2d0*k4*y(6)*sqrt(y(2))
139 0 : fin2 = -kla
140 :
141 0 : dfdy(1, 1) = (-2d0*r11 - r31) - r41
142 0 : dfdy(1, 2) = -2d0*r12
143 0 : dfdy(1, 3) = r23
144 0 : dfdy(1, 4) = r24 - r44
145 0 : dfdy(1, 5) = -r35
146 0 : dfdy(2, 1) = -0.5d0*r11 - r41
147 0 : dfdy(2, 2) = -0.5d0*r12 - 0.5d0*r52 + fin2
148 0 : dfdy(2, 4) = -r44
149 0 : dfdy(2, 6) = -0.5d0*r56
150 0 : dfdy(3, 1) = r11 + r31
151 0 : dfdy(3, 2) = r12
152 0 : dfdy(3, 3) = -r23
153 0 : dfdy(3, 4) = -r24
154 0 : dfdy(3, 5) = r35
155 0 : dfdy(4, 1) = r31 - 2d0*r41
156 0 : dfdy(4, 3) = -r23
157 0 : dfdy(4, 4) = -r24 - 2d0*r44
158 0 : dfdy(4, 5) = r35
159 0 : dfdy(5, 1) = -r31
160 0 : dfdy(5, 2) = r52
161 0 : dfdy(5, 3) = r23
162 0 : dfdy(5, 4) = r24
163 0 : dfdy(5, 5) = -r35
164 0 : dfdy(5, 6) = r56
165 0 : dfdy(6, 1) = ks*y(4)
166 0 : dfdy(6, 4) = ks*y(1)
167 0 : dfdy(6, 6) = -1d0
168 :
169 : end subroutine chemakzo_jeval
170 : !-----------------------------------------------------------------------
171 0 : subroutine chemakzo_meval(ldim, neqn, t, yprime, dfddy, ierr, rpar, ipar)
172 : integer :: ldim, neqn, ierr, ipar(*)
173 : double precision :: t, yprime(neqn), dfddy(ldim, neqn), rpar(*)
174 :
175 : integer :: i
176 0 : ierr = 0
177 :
178 0 : do i = 1, neqn - 1
179 0 : dfddy(1, i) = 1d0
180 : end do
181 :
182 0 : dfddy(1, neqn) = 0d0
183 :
184 0 : end subroutine chemakzo_meval
185 : !-----------------------------------------------------------------------
186 0 : subroutine chemakzo_solut(neqn, t, y)
187 : integer :: neqn
188 : double precision :: t, y(neqn)
189 :
190 : ! computed using true double precision on a Cray C90
191 :
192 : ! Solving Chemical Akzo Nobel problem using PSIDE
193 :
194 : ! User input:
195 :
196 : ! give relative error tolerance: 1d-19
197 : ! give absolute error tolerance: 1d-19
198 :
199 : ! Integration characteristics:
200 :
201 : ! number of integration steps 5535
202 : ! number of accepted steps 5534
203 : ! number of f evaluations 100411
204 : ! number of Jacobian evaluations 7
205 : ! number of LU decompositions 368
206 :
207 : ! CPU-time used: 30.77 sec
208 : !
209 :
210 0 : y(1) = 0.1150794920661702d0
211 0 : y(2) = 0.1203831471567715d-2
212 0 : y(3) = 0.1611562887407974d0
213 0 : y(4) = 0.3656156421249283d-3
214 0 : y(5) = 0.1708010885264404d-1
215 0 : y(6) = 0.4873531310307455d-2
216 :
217 0 : end subroutine chemakzo_solut
218 :
219 : end module bari_chemakzo
|
