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