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