Line data Source code
1 : module test_medakzo
2 :
3 : use num_def
4 : use num_lib
5 : use mtx_lib
6 : use mtx_def
7 : use test_int_support, only: i_nfcn, i_njac
8 : use utils_lib, only: mesa_error
9 : use bari_medakzo, only: medakzo_feval, medakzo_jeval, medakzo_init, medakzo_solut
10 :
11 : implicit none
12 :
13 : integer :: mljac, mujac
14 :
15 : contains
16 :
17 0 : subroutine medakzo_feval_for_blk_dble(neqn, t, y, yprime, f, ierr, rpar, ipar)
18 : integer :: neqn, ierr, ipar(:)
19 : double precision :: t, y(:), yprime(:), f(:), rpar(:)
20 :
21 : integer :: N, i, j
22 0 : double precision :: zeta, dzeta, dzeta2, k, c, phi, alpha, beta, gama, dum
23 : parameter(k=100d0, c=4d0)
24 :
25 : include 'formats'
26 :
27 0 : N = neqn/2
28 0 : dzeta = 1d0/dble(N)
29 0 : dzeta2 = dzeta*dzeta
30 0 : dum = (dzeta - 1d0)*(dzeta - 1d0)/c
31 0 : alpha = 2d0*(dzeta - 1d0)*dum/c
32 0 : beta = dum*dum
33 :
34 0 : if (t <= 5d0) then
35 : phi = 2d0
36 : else
37 0 : phi = 0d0
38 : end if
39 :
40 0 : f(1) = (phi - 2d0*y(1) + y(3))*beta/dzeta2 + alpha*(y(3) - phi)/(2d0*dzeta) - k*y(1)*y(2)
41 0 : f(2) = -k*y(1)*y(2)
42 :
43 0 : do j = 2, N - 1
44 0 : i = 2*j - 1
45 0 : zeta = j*dzeta
46 0 : dum = (zeta - 1d0)*(zeta - 1d0)/c
47 0 : alpha = 2d0*(zeta - 1d0)*dum/c
48 0 : beta = dum*dum
49 0 : gama = (y(i - 2) - 2d0*y(i) + y(i + 2))*beta/dzeta2 + alpha*(y(i + 2) - y(i - 2))/(2d0*dzeta)
50 0 : f(i) = gama - k*y(i)*y(i + 1)
51 0 : i = 2*j
52 0 : f(i) = -k*y(i)*y(i - 1)
53 : end do
54 :
55 0 : f(2*N - 1) = -k*y(2*N - 1)*y(2*N)
56 0 : f(2*N) = -k*y(2*N - 1)*y(2*N)
57 :
58 0 : return
59 : end subroutine medakzo_feval_for_blk_dble
60 :
61 0 : subroutine medakzo_jeval_for_blk_dble(ldim, neqn, t, y, yprime, dfdy, ierr, rpar, ipar)
62 : integer :: ldim, neqn, ierr, ipar(:)
63 : double precision :: t, y(:), yprime(:), dfdy(:, :), rpar(:)
64 :
65 : integer :: N, i, j
66 0 : double precision :: zeta, dzeta, dzeta2, alpha, beta, k, c, dum, bz
67 : parameter(k=100d0, c=4d0)
68 :
69 0 : do j = 1, neqn
70 0 : do i = 1, 5
71 0 : dfdy(i, j) = 0d0
72 : end do
73 : end do
74 :
75 0 : N = neqn/2
76 0 : dzeta = 1d0/dble(N)
77 0 : dzeta2 = dzeta*dzeta
78 0 : dum = (dzeta - 1d0)*(dzeta - 1d0)/c
79 0 : alpha = 2d0*(dzeta - 1d0)*dum/c
80 0 : beta = dum*dum
81 :
82 0 : dfdy(3, 1) = -beta*2d0/dzeta2 - k*y(2)
83 0 : dfdy(1, 3) = beta/dzeta2 + alpha/(2d0*dzeta)
84 0 : dfdy(2, 2) = -k*y(1)
85 0 : dfdy(4, 1) = -k*y(2)
86 0 : dfdy(3, 2) = -k*y(1)
87 :
88 0 : do j = 2, N - 1
89 0 : i = 2*j - 1
90 0 : zeta = j*dzeta
91 0 : dum = (zeta - 1d0)*(zeta - 1d0)/c
92 0 : alpha = 2d0*(zeta - 1d0)*dum/c
93 0 : beta = dum*dum
94 0 : bz = beta/dzeta2
95 0 : dfdy(5, i - 2) = bz - alpha/(2d0*dzeta)
96 0 : dfdy(3, i) = -2d0*bz - k*y(i + 1)
97 0 : dfdy(1, i + 2) = bz + alpha/(2d0*dzeta)
98 0 : dfdy(2, i + 1) = -k*y(i)
99 0 : i = 2*j
100 0 : dfdy(4, i - 1) = -k*y(i)
101 0 : dfdy(3, i) = -k*y(i - 1)
102 : end do
103 :
104 0 : dfdy(3, 2*N - 1) = -k*y(2*N)
105 0 : dfdy(2, 2*N) = -k*y(2*N - 1)
106 0 : dfdy(4, 2*N - 1) = -k*y(2*N)
107 0 : dfdy(3, 2*N) = -k*y(2*N - 1)
108 :
109 0 : return
110 : end subroutine medakzo_jeval_for_blk_dble
111 :
112 0 : subroutine medakzo_fcn_blk_dble(n, caller_id, nvar, nz, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
113 : use const_def, only: dp
114 : integer, intent(in) :: n, caller_id, nvar, nz, lrpar, lipar
115 : real(dp), intent(in) :: x, h
116 : real(dp), intent(inout), pointer :: y(:) ! (n)
117 : real(dp), intent(inout), pointer :: f(:) ! (n) ! dy/dx
118 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
119 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
120 : integer, intent(out) :: ierr
121 0 : real(dp), target :: yprime_ary(n)
122 : real(dp), pointer :: yprime(:)
123 0 : ierr = 0
124 0 : ipar(i_nfcn) = ipar(i_nfcn) + 1
125 : yprime => yprime_ary
126 0 : call medakzo_feval_for_blk_dble(n, x, y, yprime, f, ierr, rpar, ipar)
127 0 : end subroutine medakzo_fcn_blk_dble
128 :
129 30720 : subroutine medakzo_derivs(n, x, h, vars, dvars_dx, lrpar, rpar, lipar, ipar, ierr)
130 : integer, intent(in) :: n, lrpar, lipar
131 : real(dp), intent(in) :: x, h
132 : real(dp), intent(inout) :: vars(:) ! (n)
133 : real(dp), intent(inout) :: dvars_dx(:) ! (n)
134 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
135 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
136 12318720 : real(dp) :: yprime(n)
137 : integer, intent(out) :: ierr
138 : include 'formats'
139 30720 : ierr = 0
140 30720 : ipar(i_nfcn) = ipar(i_nfcn) + 1
141 30720 : call medakzo_feval(n, x, vars, yprime, dvars_dx, ierr, rpar, ipar)
142 30720 : end subroutine medakzo_derivs
143 :
144 0 : subroutine medakzo_jac_blk_dble(n, caller_id, nvar, nz, x, h, y, f, lblk1, dblk1, ublk1, lrpar, rpar, lipar, ipar, ierr)
145 : use const_def, only: dp
146 : integer, intent(in) :: n, caller_id, nvar, nz, lrpar, lipar
147 : real(dp), intent(in) :: x, h
148 : real(dp), intent(inout), pointer :: y(:) ! (n)
149 : real(dp), intent(inout), pointer :: f(:) ! (n) ! dy/dx
150 : real(dp), dimension(:), pointer, intent(inout) :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nz)
151 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
152 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
153 : integer, intent(out) :: ierr
154 :
155 0 : real(dp), dimension(:, :, :), pointer :: lblk, dblk, ublk ! =(nvar,nvar,nz)
156 : integer, parameter :: ld_dfdy = 5 ! for medakzo
157 0 : real(dp), target :: dfdy1(ld_dfdy*n)
158 : real(dp), pointer :: dfdy(:, :)
159 : !real(dp), pointer :: banded(:,:,:)
160 : integer :: i, k
161 0 : ierr = 0
162 0 : dfdy(1:ld_dfdy, 1:n) => dfdy1(1:ld_dfdy*n)
163 :
164 : ierr = 0
165 0 : ipar(i_njac) = ipar(i_njac) + 1
166 0 : call medakzo_jeval_for_blk_dble(ld_dfdy, n, x, y, f, dfdy, ierr, rpar, ipar)
167 0 : if (ierr == 0) call medakzo_fcn_blk_dble(n, caller_id, nvar, nz, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
168 :
169 : !banded(1:ld_dfdy,1:nvar,1:nz) => dfdy1(1:ld_dfdy*n)
170 0 : lblk(1:nvar, 1:nvar, 1:nz) => lblk1(1:nvar*nvar*nz)
171 0 : dblk(1:nvar, 1:nvar, 1:nz) => dblk1(1:nvar*nvar*nz)
172 0 : ublk(1:nvar, 1:nvar, 1:nz) => ublk1(1:nvar*nvar*nz)
173 :
174 : ! convert from banded to block tridiagonal
175 : ! lblk(:,:,1) is not used; ublk(:,:,nz) is not used.
176 0 : k = 1
177 0 : dblk(1, 1, k) = dfdy(3, 1) ! partial of f(1,k) wrt var(1,k) dfdy(3,i)
178 0 : dblk(1, 2, k) = dfdy(2, 2) ! partial of f(1,k) wrt var(2,k) dfdy(2,i+1)
179 0 : dblk(2, 1, k) = dfdy(4, 1) ! partial of f(2,k) wrt var(1,k) dfdy(4,i)
180 0 : dblk(2, 2, k) = dfdy(3, 2) ! partial of f(2,k) wrt var(2,k) dfdy(3,i+1)
181 0 : ublk(1, 1, k) = dfdy(1, 3) ! partial of f(1,k) wrt var(1,k+1) dfdy(1,i+2)
182 :
183 : !dfdy(1,i+2) partial of f(1,k) wrt var(1,k+1)
184 : !dfdy(2,i+1) partial of f(1,k) wrt var(2,k)
185 : !dfdy(3,i) partial of f(1,k) wrt var(1,k)
186 : !dfdy(3,i+1) partial of f(2,k) wrt var(2,k)
187 : !dfdy(4,i) partial of f(2,k) wrt var(1,k)
188 : !dfdy(5,i-2) partial of f(1,k) wrt var(1,k-1)
189 :
190 0 : do k = 2, nz - 1
191 0 : i = 2*k - 1
192 : ! set lblk
193 0 : lblk(1, 1, k) = dfdy(5, i - 2) ! partial of f(1,k) wrt var(1,k-1)
194 0 : lblk(1, 2, k) = 0 ! partial of f(1,k) wrt var(2,k-1)
195 0 : lblk(2, 1, k) = 0 ! partial of f(2,k) wrt var(1,k-1)
196 0 : lblk(2, 2, k) = 0 ! partial of f(2,k) wrt var(2,k-1)
197 : ! set dblk
198 0 : dblk(1, 1, k) = dfdy(3, i) ! partial of f(1,k) wrt var(1,k) dfdy(3,i)
199 0 : dblk(1, 2, k) = dfdy(2, i + 1) ! partial of f(1,k) wrt var(2,k) dfdy(2,i+1)
200 0 : dblk(2, 1, k) = dfdy(4, i) ! partial of f(2,k) wrt var(1,k) dfdy(4,i)
201 0 : dblk(2, 2, k) = dfdy(3, i + 1) ! partial of f(2,k) wrt var(2,k) dfdy(3,i+1)
202 : ! set ublk
203 0 : ublk(1, 1, k) = dfdy(1, i + 2) ! partial of f(1,k) wrt var(1,k+1) dfdy(1,i+2)
204 0 : ublk(2, 1, k) = 0 ! partial of f(2,k) wrt var(1,k+1)
205 0 : ublk(1, 2, k) = 0 ! partial of f(1,k) wrt var(2,k+1)
206 0 : ublk(2, 2, k) = 0 ! partial of f(2,k) wrt var(2,k+1)
207 : end do
208 :
209 0 : k = nz
210 0 : i = 2*k - 1
211 0 : dblk(1, 1, k) = dfdy(3, i) ! partial of f(1,k) wrt var(1,k)
212 0 : dblk(1, 2, k) = dfdy(2, i + 1) ! partial of f(1,k) wrt var(2,k)
213 0 : dblk(2, 1, k) = dfdy(4, i) ! partial of f(2,k) wrt var(1,k)
214 0 : dblk(2, 2, k) = dfdy(3, i + 1) ! partial of f(2,k) wrt var(2,k)
215 :
216 0 : end subroutine medakzo_jac_blk_dble
217 :
218 2430 : subroutine medakzo_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
219 : integer, intent(in) :: n, ld_dfdy, lrpar, lipar
220 : real(dp), intent(in) :: x, h
221 : real(dp), intent(inout) :: y(:)
222 : real(dp), intent(inout) :: f(:), dfdy(:, :)
223 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
224 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
225 974430 : real(dp) :: yprime(n)
226 : integer, intent(out) :: ierr
227 : include 'formats'
228 2430 : ierr = 0
229 2430 : ipar(i_njac) = ipar(i_njac) + 1
230 2430 : call medakzo_jeval(ld_dfdy, n, x, y, yprime, dfdy, ierr, rpar, ipar)
231 2430 : if (ierr == 0) call medakzo_derivs(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
232 : !write(*,*)
233 : !write(*,2)'medakzo_jacob', ipar(i_njac), x, y(1), dfdy(3,1:2)
234 2430 : end subroutine medakzo_jacob
235 :
236 0 : subroutine medakzo_sjac(n, x, h, y, f, nzmax, ia, ja, values, lrpar, rpar, lipar, ipar, ierr)
237 : ! sparse jacobian. format either compressed row or compressed column.
238 : use mtx_lib, only: band_to_row_sparse_with_diag, band_to_col_sparse_with_diag, mtx_rcond_banded
239 : use test_int_support, only: ipar_sparse_format
240 : integer, intent(in) :: n, nzmax, lrpar, lipar
241 : real(dp), intent(in) :: x, h
242 : real(dp), intent(inout) :: y(:) ! (n)
243 : real(dp), intent(inout) :: f(:) ! (n) ! dy/dx
244 : integer, intent(inout) :: ia(:) ! (n+1)
245 : integer, intent(inout) :: ja(:) ! (nzmax)
246 : real(dp), intent(inout) :: values(:) ! (nzmax)
247 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
248 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
249 : integer, intent(out) :: ierr ! nonzero means terminate integration
250 :
251 0 : real(dp) :: dfdy(n, n)
252 : integer :: ld_dfdy, nz
253 0 : ld_dfdy = n
254 : ierr = 0
255 0 : call medakzo_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
256 0 : if (ierr /= 0) return
257 0 : if (ipar(ipar_sparse_format) == 0) then
258 0 : call band_to_row_sparse_with_diag(n, mljac, mujac, dfdy, ld_dfdy, nzmax, nz, ia, ja, values, ierr)
259 : else
260 0 : call band_to_col_sparse_with_diag(n, mljac, mujac, dfdy, ld_dfdy, nzmax, nz, ia, ja, values, ierr)
261 : end if
262 0 : end subroutine medakzo_sjac
263 :
264 8 : subroutine do_test_medakzo(which_solver, which_decsol, numerical_jacobian, show_all, quiet)
265 0 : use test_support, only: show_results, show_statistics, check_results
266 : use test_int_support, only: do_test_stiff_int
267 : integer, intent(in) :: which_solver, which_decsol
268 : logical, intent(in) :: numerical_jacobian, show_all, quiet
269 :
270 : integer, parameter :: nvar = 2, nz = 200
271 : integer, parameter :: n = nvar*nz ! the number of variables in the "medakzo" system of ODEs
272 9608 : real(dp), target :: y_ary(n), yprime(n), yexact(n)
273 : real(dp), pointer :: y(:)
274 : integer, parameter :: lrpar = 1, lipar = 3, iout = 1
275 : logical :: consis
276 : integer, parameter :: ndisc = 1, n_soln = 11
277 232 : real(dp) :: result(n_soln), soln(n_soln), h0, t(0:ndisc + 1), atol(1), rtol(1)
278 : integer :: j, k, matrix_type_spec, ierr, imas, mlmas, mumas, m1, m2, itol, nstep
279 16 : real(dp), target :: rpar_ary(lrpar)
280 : integer, target :: ipar_ary(lipar)
281 : real(dp), pointer :: rpar(:)
282 : integer, pointer :: ipar(:)
283 : integer :: caller_id, nvar_blk_dble, nz_blk_dble
284 8 : real(dp), dimension(:), pointer :: lblk, dblk, ublk ! =(nvar,nvar,nz)
285 8 : real(dp), dimension(:), pointer :: uf_lblk, uf_dblk, uf_ublk ! =(nvar,nvar,nz)
286 : logical, parameter :: dbg = .false.
287 :
288 : include 'formats'
289 :
290 8 : rpar => rpar_ary
291 8 : ipar => ipar_ary
292 8 : y => y_ary
293 :
294 8 : if (.not. quiet) write (*, *) 'medakzo'
295 :
296 8 : nullify (lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk)
297 8 : caller_id = 0
298 8 : nvar_blk_dble = 0
299 8 : nz_blk_dble = 0
300 :
301 8 : t(0) = 0d0
302 : if (dbg) then
303 : t(1) = 0.05d0
304 : t(2) = 0.20d0
305 : else
306 8 : t(1) = 5d0
307 8 : t(2) = 20d0
308 : end if
309 :
310 8 : itol = 0 ! scalar tolerances
311 16 : rtol = 1d-6
312 16 : atol = 1d-6
313 8 : h0 = 1d-9 ! initial step size
314 :
315 8 : matrix_type_spec = banded_matrix_type
316 8 : mljac = 2
317 8 : mujac = 2
318 :
319 8 : imas = 0
320 8 : mlmas = 0
321 8 : mumas = 0
322 :
323 8 : m1 = 0
324 8 : m2 = 0
325 :
326 8 : call medakzo_init(n, y, yprime, consis)
327 8 : nstep = 0
328 :
329 : if (nvar_blk_dble == 0) then
330 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
331 : medakzo_derivs, medakzo_jacob, medakzo_sjac, medakzo_solout, iout, &
332 : null_fcn_blk_dble, null_jac_blk_dble, &
333 : caller_id, nvar_blk_dble, nz_blk_dble, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
334 : n, ndisc, mljac, mujac, matrix_type_spec, null_mas, imas, mlmas, mumas, m1, m2, &
335 8 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
336 : else
337 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
338 : null_fcn, null_jac, null_sjac, medakzo_solout, iout, &
339 : medakzo_fcn_blk_dble, medakzo_jac_blk_dble, &
340 : caller_id, nvar_blk_dble, nz_blk_dble, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
341 : n, ndisc, mljac, mujac, matrix_type_spec, null_mas, imas, mlmas, mumas, m1, m2, &
342 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
343 : end if
344 8 : if (ierr /= 0) then
345 0 : write (*, *) 'test_medakzo ierr', ierr
346 0 : call mesa_error(__FILE__, __LINE__)
347 : end if
348 :
349 8 : call medakzo_solut(n, 0d0, yexact)
350 8 : j = 1
351 8 : do k = 1, n/2, max(1, (n/2 - 1)/11)
352 96 : if (j > n_soln) exit
353 88 : result(j) = y(1 + 2*(k - 1))
354 88 : soln(j) = yexact(1 + 2*(k - 1))
355 96 : j = j + 1
356 : end do
357 :
358 : if (.not. dbg) then
359 8 : call check_results(n, y, yexact, rtol(1)*50, ierr)
360 8 : if (ierr /= 0) then
361 0 : write (*, *) 'check results ierr', ierr
362 : !call mesa_error(__FILE__,__LINE__) ! do_test_medakzo
363 : end if
364 : end if
365 :
366 8 : if (quiet) return
367 :
368 8 : call show_results(n_soln, result, soln, show_all)
369 8 : call show_statistics(ipar(i_nfcn), ipar(i_njac), nstep, show_all)
370 :
371 16 : end subroutine do_test_medakzo
372 :
373 4876 : subroutine medakzo_solout(nr, xold, x, n, y, rwork, iwork, interp_y, lrpar, rpar, lipar, ipar, irtrn)
374 : ! nr is the step number.
375 : ! x is the current x value; xold is the previous x value.
376 : ! y is the current y value.
377 : ! irtrn negative means terminate integration.
378 : ! rwork and iwork hold info for
379 : integer, intent(in) :: nr, n, lrpar, lipar
380 : real(dp), intent(in) :: xold, x
381 : real(dp), intent(inout) :: y(:) ! (n)
382 : real(dp), intent(inout), target :: rwork(*)
383 : integer, intent(inout), target :: iwork(*)
384 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
385 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
386 : interface
387 : ! this subroutine can be called from your solout routine.
388 : ! it computes interpolated values for y components during the just completed step.
389 : real(dp) function interp_y(i, s, rwork, iwork, ierr)
390 : use const_def, only: dp
391 : implicit none
392 : integer, intent(in) :: i ! result is interpolated approximation of y(i) at x=s.
393 : real(dp), intent(in) :: s ! interpolation x value (between xold and x).
394 : real(dp), intent(inout), target :: rwork(*)
395 : integer, intent(inout), target :: iwork(*)
396 : integer, intent(out) :: ierr
397 : end function interp_y
398 : end interface
399 : integer, intent(out) :: irtrn
400 : integer :: ierr
401 : include 'formats'
402 : !if (mod(nr,10) == 0) write(*,2) 'step', nr, x, y(1:2)
403 : !if (nr >= 100) stop
404 4876 : ierr = 0
405 4876 : irtrn = 0
406 8 : end subroutine medakzo_solout
407 :
408 : end module test_medakzo
|
