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