Line data Source code
1 : module test_chemakzo
2 :
3 : use num_def
4 : use num_lib
5 : use test_int_support, only: i_nfcn, i_njac
6 : use utils_lib, only: mesa_error
7 : use bari_chemakzo, only: chemakzo_feval, chemakzo_jeval, chemakzo_init, chemakzo_meval, chemakzo_solut
8 :
9 : implicit none
10 :
11 : contains
12 :
13 0 : subroutine chemakzo_derivs(n, x, h, vars, dvars_dx, lrpar, rpar, lipar, ipar, ierr)
14 : integer, intent(in) :: n, lrpar, lipar
15 : real(dp), intent(in) :: x, h
16 : real(dp), intent(inout) :: vars(:) ! (n)
17 : real(dp), intent(inout) :: dvars_dx(:) ! (n)
18 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
19 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
20 0 : real(dp) :: yprime(n)
21 : integer, intent(out) :: ierr
22 0 : ierr = 0
23 0 : ipar(i_nfcn) = ipar(i_nfcn) + 1
24 0 : call chemakzo_feval(n, x, vars, yprime, dvars_dx, ierr, rpar, ipar)
25 0 : end subroutine chemakzo_derivs
26 :
27 0 : subroutine chemakzo_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
28 : integer, intent(in) :: n, ld_dfdy, lrpar, lipar
29 : real(dp), intent(in) :: x, h
30 : real(dp), intent(inout) :: y(:) ! (n)
31 : real(dp), intent(inout) :: f(:), dfdy(:, :) ! (ld_dfdy,n)
32 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
33 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
34 0 : real(dp) :: yprime(n)
35 : integer, intent(out) :: ierr
36 0 : ierr = 0
37 0 : ipar(i_njac) = ipar(i_njac) + 1
38 0 : call chemakzo_jeval(ld_dfdy, n, x, y, yprime, dfdy, ierr, rpar, ipar)
39 0 : if (ierr == 0) call chemakzo_derivs(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
40 0 : end subroutine chemakzo_jacob
41 :
42 0 : subroutine chemakzo_sjac(n, x, h, y, f, nzmax, ia, ja, values, lrpar, rpar, lipar, ipar, ierr)
43 : ! sparse jacobian. format either compressed row or compressed column.
44 : use mtx_lib, only: dense_to_row_sparse_with_diag, dense_to_col_sparse_with_diag
45 : use test_int_support, only: ipar_sparse_format
46 : integer, intent(in) :: n, nzmax, lrpar, lipar
47 : real(dp), intent(in) :: x, h
48 : real(dp), intent(inout) :: y(:) ! (n)
49 : real(dp), intent(inout) :: f(:) ! (n) ! dy/dx
50 : integer, intent(inout) :: ia(:) ! (n+1)
51 : integer, intent(inout) :: ja(:) ! (nzmax)
52 : real(dp), intent(inout) :: values(:) ! (nzmax)
53 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
54 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
55 : integer, intent(out) :: ierr ! nonzero means terminate integration
56 0 : real(dp) :: dfdy(n, n)
57 : integer :: ld_dfdy, nz
58 0 : ld_dfdy = n
59 : ierr = 0
60 0 : call chemakzo_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
61 0 : if (ierr /= 0) return
62 0 : if (ipar(ipar_sparse_format) == 0) then
63 0 : call dense_to_row_sparse_with_diag(n, n, dfdy, nzmax, nz, ia, ja, values, ierr)
64 : else
65 0 : call dense_to_col_sparse_with_diag(n, n, dfdy, nzmax, nz, ia, ja, values, ierr)
66 : end if
67 0 : end subroutine chemakzo_sjac
68 :
69 0 : subroutine chemakzo_solout(nr, xold, x, n, y, rwork, iwork, interp_y, lrpar, rpar, lipar, ipar, irtrn)
70 : ! nr is the step number.
71 : ! x is the current x value; xold is the previous x value.
72 : ! y is the current y value.
73 : ! irtrn negative means terminate integration.
74 : ! rwork and iwork hold info for
75 : integer, intent(in) :: nr, n, lrpar, lipar
76 : real(dp), intent(in) :: xold, x
77 : real(dp), intent(inout) :: y(:) ! (n)
78 : real(dp), intent(inout), target :: rwork(*)
79 : integer, intent(inout), target :: iwork(*)
80 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
81 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
82 : interface
83 : ! this subroutine can be called from your solout routine.
84 : ! it computes interpolated values for y components during the just completed step.
85 : real(dp) function interp_y(i, s, rwork, iwork, ierr)
86 : use const_def, only: dp
87 : implicit none
88 : integer, intent(in) :: i ! result is interpolated approximation of y(i) at x=s.
89 : real(dp), intent(in) :: s ! interpolation x value (between xold and x).
90 : real(dp), intent(inout), target :: rwork(*)
91 : integer, intent(inout), target :: iwork(*)
92 : integer, intent(out) :: ierr
93 : end function interp_y
94 : end interface
95 : integer, intent(out) :: irtrn
96 0 : real(dp) :: xout, val
97 : integer :: i, ierr
98 : include 'formats'
99 0 : irtrn = 0
100 : !write(*,2) 'x', nr, x
101 0 : xout = 100
102 0 : if (x >= xout .and. xout > xold) then
103 0 : write (*, *)
104 0 : write (*, 1) 'dense output for x=', xout
105 : ierr = 0
106 0 : do i = 1, n
107 0 : val = interp_y(i, xout, rwork, iwork, ierr)
108 0 : if (ierr /= 0) then
109 0 : write (*, 2) 'interp_y failed', i, val
110 : else
111 0 : write (*, 2) 'val', i, val
112 : end if
113 : end do
114 0 : write (*, *)
115 : end if
116 :
117 0 : end subroutine chemakzo_solout
118 :
119 0 : subroutine chemakzo_mas_band(n, am, lmas, lrpar, rpar, lipar, ipar)
120 : integer, intent(in) :: n, lmas, lrpar, lipar
121 : real(dp), intent(inout) :: am(:, :) ! (lmas,n)
122 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
123 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
124 0 : real(dp) :: yprime(n), t
125 : integer :: ierr
126 0 : ierr = 0
127 0 : call chemakzo_meval(lmas, n, t, yprime, am, ierr, rpar, ipar)
128 0 : end subroutine chemakzo_mas_band
129 :
130 0 : subroutine chemakzo_mas_full(n, am, lmas, lrpar, rpar, lipar, ipar)
131 : integer, intent(in) :: n, lmas, lrpar, lipar
132 : real(dp), intent(inout) :: am(:, :) ! (lmas,n)
133 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
134 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
135 0 : real(dp) :: yprime(n), t
136 : integer :: ierr, i
137 0 : ierr = 0
138 0 : call chemakzo_meval(lmas, n, t, yprime, am, ierr, rpar, ipar)
139 0 : do i = 2, n
140 0 : am(i, i) = am(1, i)
141 0 : am(1, i) = 0
142 : end do
143 0 : end subroutine chemakzo_mas_full
144 :
145 0 : subroutine do_test_chemakzo(which_solver, which_decsol, m_band, numerical_jacobian, show_all, quiet)
146 : use test_support, only: show_results, show_statistics, check_results
147 : use test_int_support, only: do_test_stiff_int
148 : integer, intent(in) :: which_solver, which_decsol
149 : logical, intent(in) :: m_band, numerical_jacobian, show_all, quiet
150 :
151 : integer, parameter :: n = 6 ! the number of variables in the "chemakzo" system of ODEs
152 0 : real(dp), target :: y_ary(n), yprime(n), yexact(n)
153 : real(dp), pointer :: y(:)
154 : integer, parameter :: lrpar = 1, lipar = 3, iout = 2 ! test dense output
155 : logical :: consis
156 : integer, parameter :: ndisc = 0
157 0 : real(dp) :: h0, t(0:ndisc + 1), atol(1), rtol(1)
158 : integer :: mujac, mljac, matrix_type_spec, ierr, imas, mlmas, mumas, m1, m2, itol, nstep
159 0 : real(dp), target :: rpar_ary(lrpar)
160 : integer, target :: ipar_ary(lipar)
161 : real(dp), pointer :: rpar(:)
162 : integer, pointer :: ipar(:)
163 : integer :: caller_id, nvar, nz
164 0 : real(dp), dimension(:), pointer :: lblk, dblk, ublk ! =(nvar,nvar,nz)
165 0 : real(dp), dimension(:), pointer :: uf_lblk, uf_dblk, uf_ublk ! =(nvar,nvar,nz)
166 :
167 0 : nullify (lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk)
168 0 : caller_id = 0
169 0 : nvar = 0
170 0 : nz = 0
171 :
172 0 : rpar => rpar_ary
173 0 : ipar => ipar_ary
174 0 : y => y_ary
175 :
176 0 : if (.not. quiet) write (*, *) 'chemakzo'
177 :
178 0 : t(0) = 0
179 0 : t(1) = 180d0
180 :
181 0 : itol = 0 ! scalar tolerances
182 0 : rtol(1) = 1d-8
183 0 : atol(1) = 1d-8
184 0 : h0 = 1d-10 ! initial step size
185 :
186 0 : mljac = n ! square matrix
187 0 : mujac = n
188 0 : matrix_type_spec = square_matrix_type
189 :
190 0 : imas = 1
191 0 : m1 = 0
192 0 : m2 = 0
193 :
194 0 : call chemakzo_init(n, y, yprime, consis)
195 0 : nstep = 0
196 0 : if (m_band) then
197 0 : write (*, *) 'M banded'
198 : ! mass matrix is diagonal
199 0 : mlmas = 0
200 0 : mumas = 0
201 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
202 : chemakzo_derivs, chemakzo_jacob, chemakzo_sjac, chemakzo_solout, iout, &
203 : null_fcn_blk_dble, null_jac_blk_dble, &
204 : caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
205 : n, ndisc, mljac, mujac, matrix_type_spec, chemakzo_mas_band, imas, mlmas, mumas, m1, m2, &
206 0 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
207 : else
208 0 : write (*, *) 'M full'
209 0 : mlmas = n
210 0 : mumas = n
211 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
212 : chemakzo_derivs, chemakzo_jacob, chemakzo_sjac, chemakzo_solout, iout, &
213 : null_fcn_blk_dble, null_jac_blk_dble, &
214 : caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
215 : n, ndisc, mljac, mujac, matrix_type_spec, chemakzo_mas_full, imas, mlmas, mumas, m1, m2, &
216 0 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
217 : end if
218 0 : if (ierr /= 0) then
219 0 : write (*, *) 'chemakzo ierr', ierr
220 0 : call mesa_error(__FILE__, __LINE__)
221 : end if
222 :
223 0 : call chemakzo_solut(n, 0d0, yexact)
224 0 : call check_results(n, y, yexact, rtol(1)*10, ierr)
225 0 : if (ierr /= 0) then
226 0 : write (*, *) 'check results ierr', ierr
227 0 : call mesa_error(__FILE__, __LINE__) ! do_test_vdpol
228 : end if
229 :
230 0 : if (quiet) return
231 :
232 0 : call show_results(n, y, yexact, show_all)
233 0 : call show_statistics(ipar(i_nfcn), ipar(i_njac), nstep, show_all)
234 :
235 0 : end subroutine do_test_chemakzo
236 :
237 : end module test_chemakzo
|
