Line data Source code
1 : module test_vdpol
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_vdpol, only: vdpol_feval, vdpol_jeval, vdpol_init, vdpol_solut
9 :
10 : implicit none
11 :
12 : logical, parameter :: dbg = .false.
13 :
14 : integer :: cnt = 0
15 :
16 : contains
17 :
18 14 : subroutine vdpol_fcn_blk_dble(n, caller_id, nvar, nz, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
19 : use const_def, only: dp
20 : integer, intent(in) :: n, caller_id, nvar, nz, lrpar, lipar
21 : real(dp), intent(in) :: x, h
22 : real(dp), intent(inout), pointer :: y(:) ! (n)
23 : real(dp), intent(inout), pointer :: f(:) ! (n) ! dy/dx
24 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
25 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
26 : integer, intent(out) :: ierr
27 : include 'formats'
28 0 : call vdpol_derivs(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
29 : !write(*,2) 'vdpol_fcn_blk_dble', ipar(i_nfcn), x
30 0 : end subroutine vdpol_fcn_blk_dble
31 :
32 59826 : subroutine vdpol_derivs(n, x, h, vars, dvars_dx, lrpar, rpar, lipar, ipar, ierr)
33 : integer, intent(in) :: n, lrpar, lipar
34 : real(dp), intent(in) :: x, h
35 : real(dp), intent(inout) :: vars(:) ! (n)
36 : real(dp), intent(inout) :: dvars_dx(:) ! (n)
37 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
38 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
39 179478 : real(dp) :: yprime(n)
40 : integer, intent(out) :: ierr
41 59826 : ierr = 0
42 59826 : ipar(i_nfcn) = ipar(i_nfcn) + 1
43 59826 : call vdpol_feval(n, x, vars, yprime, dvars_dx, ierr, rpar, ipar)
44 :
45 : if (dbg) then
46 : cnt = cnt + 1
47 : write (*, *) 'func cnt', cnt
48 : write (*, *) 'func n', n
49 : write (*, *) 'func x', x
50 : write (*, *) 'func y(1)', vars(1)
51 : write (*, *) 'func y(2)', vars(2)
52 : write (*, *) 'func f(1)', dvars_dx(1)
53 : write (*, *) 'func f(2)', dvars_dx(2)
54 : write (*, *)
55 : if (cnt == 5) stop 'vdpol_derivs'
56 : end if
57 :
58 59826 : end subroutine vdpol_derivs
59 :
60 0 : subroutine vdpol_jac_blk_dble(n, caller_id, nvar, nz, x, h, y, f, lblk1, dblk1, ublk1, lrpar, rpar, lipar, ipar, ierr)
61 : use const_def, only: dp
62 : integer, intent(in) :: n, caller_id, nvar, nz, lrpar, lipar
63 : real(dp), intent(in) :: x, h
64 : real(dp), intent(inout), pointer :: y(:) ! (n)
65 : real(dp), intent(inout), pointer :: f(:) ! (n) ! dy/dx
66 : real(dp), dimension(:), pointer, intent(inout) :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nz)
67 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
68 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
69 : integer, intent(out) :: ierr
70 :
71 0 : real(dp), dimension(:, :, :), pointer :: lblk, dblk, ublk ! =(nvar,nvar,nz)
72 : integer, parameter :: ld_dfdy = 2 ! for vdpol
73 0 : real(dp), target :: dfdy1(ld_dfdy*n)
74 : real(dp), pointer :: dfdy(:, :)
75 : include 'formats'
76 0 : ierr = 0
77 0 : dfdy(1:ld_dfdy, 1:n) => dfdy1(1:ld_dfdy*n)
78 0 : call vdpol_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
79 : !write(*,2) 'vdpol_jac_blk_dble', ipar(i_nfcn), x
80 0 : if (ierr /= 0) return
81 0 : lblk(1:nvar, 1:nvar, 1:nz) => lblk1(1:nvar*nvar*nz)
82 0 : dblk(1:nvar, 1:nvar, 1:nz) => dblk1(1:nvar*nvar*nz)
83 0 : ublk(1:nvar, 1:nvar, 1:nz) => ublk1(1:nvar*nvar*nz)
84 : ! convert from square to block tridiagonal
85 0 : dblk(:, :, 1) = dfdy(:, :)
86 0 : ublk(:, :, 1) = 0
87 0 : lblk(:, :, 1) = 0
88 0 : end subroutine vdpol_jac_blk_dble
89 :
90 8221 : subroutine vdpol_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
91 : integer, intent(in) :: n, ld_dfdy, lrpar, lipar
92 : real(dp), intent(in) :: x, h
93 : real(dp), intent(inout) :: y(:)
94 : real(dp), intent(inout) :: f(:), dfdy(:, :)
95 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
96 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
97 24663 : real(dp) :: yprime(n)
98 : integer, intent(out) :: ierr
99 8221 : ierr = 0
100 8221 : ipar(i_njac) = ipar(i_njac) + 1
101 8221 : call vdpol_jeval(ld_dfdy, n, x, y, yprime, dfdy, ierr, rpar, ipar)
102 8221 : if (ierr == 0) call vdpol_derivs(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
103 :
104 : if (.false.) then
105 : write (*, *) 'jac_fcn y(1)', y(1)
106 : write (*, *) 'jac_fcn y(2)', y(2)
107 : write (*, *) 'jac_fcn dfdy(1,1)', dfdy(1, 1)
108 : write (*, *) 'jac_fcn dfdy(1,2)', dfdy(1, 2)
109 : write (*, *) 'jac_fcn dfdy(2,1)', dfdy(2, 1)
110 : write (*, *) 'jac_fcn dfdy(2,2)', dfdy(2, 2)
111 : write (*, *)
112 : end if
113 :
114 8221 : end subroutine vdpol_jacob
115 :
116 0 : subroutine vdpol_sjac(n, x, h, y, f, nzmax, ia, ja, values, lrpar, rpar, lipar, ipar, ierr)
117 : ! sparse jacobian. format either compressed row or compressed column.
118 : use mtx_lib, only: dense_to_row_sparse_with_diag, dense_to_col_sparse_with_diag
119 : use test_int_support, only: ipar_sparse_format
120 : integer, intent(in) :: n, nzmax, lrpar, lipar
121 : real(dp), intent(in) :: x, h
122 : real(dp), intent(inout) :: y(:) ! (n)
123 : real(dp), intent(inout) :: f(:) ! (n) ! dy/dx
124 : integer, intent(inout) :: ia(:) ! (n+1)
125 : integer, intent(inout) :: ja(:) ! (nzmax)
126 : real(dp), intent(inout) :: values(:) ! (nzmax)
127 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
128 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
129 : integer, intent(out) :: ierr ! nonzero means terminate integration
130 0 : real(dp) :: dfdy(n, n)
131 : integer :: ld_dfdy, nz
132 0 : ld_dfdy = n
133 : ierr = 0
134 0 : call vdpol_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
135 0 : if (ierr /= 0) return
136 0 : if (ipar(ipar_sparse_format) == 0) then
137 0 : call dense_to_row_sparse_with_diag(n, n, dfdy, nzmax, nz, ia, ja, values, ierr)
138 : else
139 0 : call dense_to_col_sparse_with_diag(n, n, dfdy, nzmax, nz, ia, ja, values, ierr)
140 : end if
141 0 : end subroutine vdpol_sjac
142 :
143 0 : subroutine vdpol_solout(nr, xold, x, n, y, rwork, iwork, interp_y, lrpar, rpar, lipar, ipar, irtrn)
144 : ! nr is the step number.
145 : ! x is the current x value; xold is the previous x value.
146 : ! y is the current y value.
147 : ! irtrn negative means terminate integration.
148 : ! rwork and iwork hold info for
149 : integer, intent(in) :: nr, n, lrpar, lipar
150 : real(dp), intent(in) :: xold, x
151 : real(dp), intent(inout) :: y(:) ! (n)
152 : real(dp), intent(inout), target :: rwork(*)
153 : integer, intent(inout), target :: iwork(*)
154 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
155 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
156 : interface
157 : ! this subroutine can be called from your solout routine.
158 : ! it computes interpolated values for y components during the just completed step.
159 : real(dp) function interp_y(i, s, rwork, iwork, ierr)
160 : use const_def, only: dp
161 : implicit none
162 : integer, intent(in) :: i ! result is interpolated approximation of y(i) at x=s.
163 : real(dp), intent(in) :: s ! interpolation x value (between xold and x).
164 : real(dp), intent(inout), target :: rwork(*)
165 : integer, intent(inout), target :: iwork(*)
166 : integer, intent(out) :: ierr
167 : end function interp_y
168 : end interface
169 : integer, intent(out) :: irtrn
170 0 : real(dp) :: xout, y1, y2
171 : integer :: ierr
172 :
173 : if (dbg .and. nr > 450) stop
174 :
175 0 : ierr = 0
176 0 : irtrn = 0
177 0 : xout = rpar(1)
178 0 : if (nr == 1) then
179 0 : write (6, 99) x, y(1), y(2), nr - 1
180 0 : xout = 0.2d0
181 : else
182 0 : do
183 0 : if (x >= xout) then
184 0 : y1 = interp_y(1, xout, rwork, iwork, ierr)
185 0 : if (ierr /= 0) exit
186 0 : y2 = interp_y(2, xout, rwork, iwork, ierr)
187 0 : write (6, 99) xout, y1, y2, nr - 1
188 0 : if (ierr /= 0) exit
189 0 : xout = xout + 0.2d0
190 : cycle
191 : end if
192 0 : exit
193 : end do
194 : end if
195 0 : if (ierr /= 0) then
196 0 : write (*, *) 'problem with interp_y in vdpol_solout'
197 0 : irtrn = -1
198 : end if
199 0 : rpar(1) = xout
200 : 99 format(1x, 'x =', f5.2, ' y =', 2e18.10, ' nstep =', i8)
201 0 : end subroutine vdpol_solout
202 :
203 14 : subroutine do_test_vdpol(which_solver, which_decsol, numerical_jacobian, show_all, quiet)
204 : use test_support, only: show_results, show_statistics, check_results
205 : use test_int_support, only: do_test_stiff_int
206 : integer, intent(in) :: which_solver, which_decsol
207 : logical, intent(in) :: numerical_jacobian, show_all, quiet
208 :
209 : integer, parameter :: nvar = 2, nz = 1
210 : integer, parameter :: n = nvar*nz ! the number of variables in the "vdpol" system of ODEs
211 98 : real(dp), target :: y_ary(n), yprime(n), yexact(n)
212 : real(dp), pointer :: y(:)
213 : integer, parameter :: lrpar = 1, lipar = 3
214 : logical :: consis
215 : integer, parameter :: ndisc = 0
216 84 : real(dp) :: h0, t(0:ndisc + 1), atol(1), rtol(1)
217 : integer :: mujac, mljac, matrix_type_spec, ierr, imas, mlmas, mumas, m1, m2, itol, iout, nstep
218 28 : real(dp), target :: rpar_ary(lrpar)
219 : integer, target :: ipar_ary(lipar)
220 : real(dp), pointer :: rpar(:)
221 : integer, pointer :: ipar(:)
222 : integer :: caller_id, nvar_blk_dble, nz_blk_dble
223 14 : real(dp), dimension(:), pointer :: lblk, dblk, ublk ! =(nvar,nvar,nz)
224 14 : real(dp), dimension(:), pointer :: uf_lblk, uf_dblk, uf_ublk ! =(nvar,nvar,nz)
225 :
226 14 : rpar => rpar_ary
227 14 : ipar => ipar_ary
228 14 : y => y_ary
229 :
230 14 : if (.not. quiet) write (*, *) 'vdpol'
231 :
232 14 : nullify (lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk)
233 14 : caller_id = 0
234 14 : nvar_blk_dble = 0
235 14 : nz_blk_dble = 0
236 :
237 14 : t(0) = 0
238 14 : t(1) = 2d0
239 :
240 14 : itol = 0 ! scalar tolerances
241 : rtol(1) = 1d-8
242 : atol(1) = 1d-8
243 : h0 = 1d-10 ! initial step size
244 :
245 : rtol(1) = 1d-6
246 : atol(1) = 1d-6
247 : h0 = 1d-8 ! initial step size
248 :
249 14 : rtol(1) = 1d-4
250 14 : atol(1) = 1d-4
251 14 : h0 = 1d-4 ! initial step size
252 :
253 14 : mljac = n ! square matrix
254 14 : mujac = n
255 14 : matrix_type_spec = square_matrix_type
256 :
257 14 : imas = 0
258 14 : mlmas = 0
259 14 : mumas = 0
260 :
261 14 : m1 = 0
262 14 : m2 = 0
263 :
264 14 : if (show_all) then
265 0 : iout = 1
266 : else
267 14 : iout = 0
268 : end if
269 :
270 56 : ipar = 0
271 :
272 14 : call vdpol_init(n, y, yprime, consis)
273 14 : nstep = 0
274 :
275 : if (nvar_blk_dble == 0) then
276 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
277 : vdpol_derivs, vdpol_jacob, vdpol_sjac, vdpol_solout, iout, &
278 : null_fcn_blk_dble, null_jac_blk_dble, &
279 : caller_id, nvar_blk_dble, nz_blk_dble, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
280 : n, ndisc, mljac, mujac, matrix_type_spec, null_mas, imas, mlmas, mumas, m1, m2, &
281 14 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
282 : else
283 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
284 : null_fcn, null_jac, null_sjac, vdpol_solout, iout, &
285 : vdpol_fcn_blk_dble, vdpol_jac_blk_dble, &
286 : caller_id, nvar_blk_dble, nz_blk_dble, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
287 : n, ndisc, mljac, mujac, matrix_type_spec, null_mas, imas, mlmas, mumas, m1, m2, &
288 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
289 : end if
290 14 : if (ierr /= 0) then
291 0 : write (*, *) 'test_vdpol ierr', ierr
292 0 : call mesa_error(__FILE__, __LINE__)
293 : end if
294 :
295 14 : call vdpol_solut(n, 0d0, yexact)
296 : !call check_results(n,y,yexact,rtol(1)*2,ierr)
297 14 : if (ierr /= 0) then
298 0 : write (*, *) 'check results ierr', ierr
299 0 : call mesa_error(__FILE__, __LINE__) ! do_test_vdpol
300 : end if
301 :
302 14 : if (quiet) return
303 :
304 14 : call show_results(n, y, yexact, show_all)
305 14 : call show_statistics(ipar(i_nfcn), ipar(i_njac), nstep, show_all)
306 :
307 28 : end subroutine do_test_vdpol
308 :
309 : end module test_vdpol
|
