Line data Source code
1 : module test_diffusion
2 : use num_def
3 : use num_lib
4 : use math_lib
5 : use test_int_support, only: i_nfcn, i_njac
6 : use utils_lib, only: mesa_error
7 :
8 : implicit none
9 :
10 : integer :: mljac, mujac, nstep
11 :
12 : integer, parameter :: nz = 48
13 : integer, parameter :: diff_mujac = 1, diff_mljac = 1, diff_ldjac = diff_mujac + diff_mljac + 1
14 :
15 : real(dp) :: y0(nz), yexact(nz), dfdy2(4, nz), rcond, work(3*nz)
16 : integer :: ipiv(nz), iwork(nz), info
17 :
18 : real(dp) :: sig_dm(nz)
19 :
20 : contains
21 :
22 7847 : subroutine diffusion_op(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
23 : integer, intent(in) :: n, ld_dfdy, lrpar, lipar
24 : real(dp), intent(in) :: x, h
25 : real(dp), intent(inout) :: y(n)
26 : real(dp), intent(inout) :: f(n), dfdy(ld_dfdy, n)
27 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
28 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
29 : integer, intent(out) :: ierr
30 :
31 7847 : real(dp) :: sig1, sig2
32 : integer :: i
33 7847 : ierr = 0
34 :
35 1204894 : f = 0; dfdy = 0
36 :
37 : sig2 = 0
38 384503 : do i = 1, n
39 376656 : sig1 = sig2
40 376656 : if (i < n) then
41 368809 : sig2 = sig_dm(i + 1)
42 : else
43 : sig2 = 0
44 : end if
45 376656 : if (ld_dfdy > 0) then
46 145296 : if (i > 1) then
47 142269 : dfdy(3, i - 1) = sig1
48 : end if
49 145296 : dfdy(2, i) = -(sig1 + sig2)
50 145296 : if (i < n) then
51 142269 : dfdy(1, i + 1) = sig2
52 : end if
53 : end if
54 384503 : if (i == n) then
55 7847 : f(i) = sig1*(y(i - 1) - y(i))
56 368809 : else if (i == 1) then
57 7847 : f(i) = -sig2*(y(i) - y(i + 1))
58 : else
59 360962 : f(i) = sig1*(y(i - 1) - y(i)) - sig2*(y(i) - y(i + 1))
60 : end if
61 : end do
62 :
63 7847 : end subroutine diffusion_op
64 :
65 4820 : subroutine diffusion_derivs(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
66 : integer, intent(in) :: n, lrpar, lipar
67 : real(dp), intent(in) :: x, h
68 : real(dp), intent(inout) :: y(:) ! (n)
69 : real(dp), intent(inout) :: f(:) ! (n)
70 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
71 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
72 : integer, intent(out) :: ierr
73 :
74 9640 : real(dp) :: dfdy(0, n)
75 :
76 4820 : ierr = 0
77 4820 : ipar(i_nfcn) = ipar(i_nfcn) + 1
78 4820 : call diffusion_op(n, x, h, y, f, dfdy, 0, lrpar, rpar, lipar, ipar, ierr)
79 4820 : end subroutine diffusion_derivs
80 :
81 3027 : subroutine diffusion_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
82 : use mtx_lib, only: mtx_rcond_banded
83 : integer, intent(in) :: n, ld_dfdy, lrpar, lipar
84 : real(dp), intent(in) :: x, h
85 : real(dp), intent(inout) :: y(:)
86 : real(dp), intent(inout) :: f(:), dfdy(:, :)
87 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
88 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
89 : integer, intent(out) :: ierr
90 : include 'formats'
91 3027 : ierr = 0
92 3027 : ipar(i_njac) = ipar(i_njac) + 1
93 3027 : call diffusion_op(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
94 :
95 : return
96 :
97 : dfdy2(2, 1:n) = dfdy(1, 1:n)
98 : dfdy2(3, 1:n) = dfdy(2, 1:n) - 1
99 : dfdy2(4, 1:n) = dfdy(3, 1:n)
100 :
101 : call mtx_rcond_banded('N', n, n, 1, 1, dfdy2, 4, ipiv, rcond, work, iwork, info)
102 : write (*, 2) 'diffusion_jacob rcond', info, x, safe_log10(rcond)
103 :
104 : end subroutine diffusion_jacob
105 :
106 0 : subroutine diffusion_sjac(n, x, h, y, f, nzmax, ia, ja, values, lrpar, rpar, lipar, ipar, ierr)
107 : ! sparse jacobian. format either compressed row or compressed column.
108 : use mtx_lib, only: band_to_row_sparse_with_diag, band_to_col_sparse_with_diag, mtx_rcond_banded
109 : use test_int_support, only: ipar_sparse_format
110 : integer, intent(in) :: n, nzmax, lrpar, lipar
111 : real(dp), intent(in) :: x, h
112 : real(dp), intent(inout) :: y(:) ! (n)
113 : real(dp), intent(inout) :: f(:) ! (n) ! dy/dx
114 : integer, intent(inout) :: ia(:) ! (n+1)
115 : integer, intent(inout) :: ja(:) ! (nzmax)
116 : real(dp), intent(inout) :: values(:) ! (nzmax)
117 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
118 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
119 : integer, intent(out) :: ierr ! nonzero means terminate integration
120 :
121 0 : real(dp) :: dfdy(n, n)
122 : integer :: ld_dfdy, nz
123 0 : ld_dfdy = n
124 : ierr = 0
125 0 : call diffusion_jacob(n, x, h, y, f, dfdy, ld_dfdy, lrpar, rpar, lipar, ipar, ierr)
126 0 : if (ierr /= 0) return
127 0 : if (ipar(ipar_sparse_format) == 0) then
128 0 : call band_to_row_sparse_with_diag(n, mljac, mujac, dfdy, ld_dfdy, nzmax, nz, ia, ja, values, ierr)
129 : else
130 0 : call band_to_col_sparse_with_diag(n, mljac, mujac, dfdy, ld_dfdy, nzmax, nz, ia, ja, values, ierr)
131 : end if
132 0 : end subroutine diffusion_sjac
133 :
134 3034 : subroutine diffusion_solout(nr, xold, x, n, y, rwork, iwork, interp_y, lrpar, rpar, lipar, ipar, irtrn)
135 : ! nr is the step number.
136 : ! x is the current x value; xold is the previous x value.
137 : ! y is the current y value.
138 : ! irtrn negative means terminate integration.
139 : ! rwork and iwork hold info for
140 : integer, intent(in) :: nr, n, lrpar, lipar
141 : real(dp), intent(in) :: xold, x
142 : real(dp), intent(inout) :: y(:) ! (n)
143 : real(dp), intent(inout), target :: rwork(*)
144 : integer, intent(inout), target :: iwork(*)
145 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
146 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
147 : interface
148 : ! this subroutine can be called from your solout routine.
149 : ! it computes interpolated values for y components during the just completed step.
150 : real(dp) function interp_y(i, s, rwork, iwork, ierr)
151 : use const_def, only: dp
152 : implicit none
153 : integer, intent(in) :: i ! result is interpolated approximation of y(i) at x=s.
154 : real(dp), intent(in) :: s ! interpolation x value (between xold and x).
155 : real(dp), intent(inout), target :: rwork(*)
156 : integer, intent(inout), target :: iwork(*)
157 : integer, intent(out) :: ierr
158 : end function interp_y
159 : end interface
160 : integer, intent(out) :: irtrn
161 3034 : irtrn = 0
162 0 : end subroutine diffusion_solout
163 :
164 7 : subroutine do_test_diffusion(which_solver, which_decsol, numerical_jacobian, show_all, quiet)
165 : use test_support, only: show_results, show_statistics, check_results
166 : use test_int_support, only: do_test_stiff_int
167 : integer, intent(in) :: which_solver, which_decsol
168 : logical, intent(in) :: numerical_jacobian, show_all, quiet
169 :
170 : real(dp), parameter :: sig = 1d-2
171 :
172 : real(dp), parameter :: ystart = 1d0, tend = 1d4
173 : integer, parameter :: n = nz
174 : integer, parameter :: lrpar = 1, lipar = 3, iout = 1
175 : integer, parameter :: ndisc = 0, n_soln = 10
176 343 : real(dp), target :: y_ary(n)
177 : real(dp), pointer :: y(:)
178 196 : real(dp) :: result(n_soln), soln(n_soln), h0, atol(1), rtol(1), t(0:ndisc + 1)
179 : integer :: i, k, matrix_type_spec, ierr, imas, mlmas, mumas, m1, m2, itol
180 14 : real(dp), target :: rpar_ary(lrpar)
181 : integer, target :: ipar_ary(lipar)
182 : real(dp), pointer :: rpar(:)
183 : integer, pointer :: ipar(:)
184 : integer :: caller_id, nvar_blk_dble, nz_blk_dble
185 7 : real(dp), dimension(:), pointer :: lblk, dblk, ublk ! =(nvar,nvar,nz)
186 7 : real(dp), dimension(:), pointer :: uf_lblk, uf_dblk, uf_ublk ! =(nvar,nvar,nz)
187 :
188 7 : nullify (lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk)
189 7 : caller_id = 0
190 7 : nvar_blk_dble = 0
191 7 : nz_blk_dble = 0
192 :
193 7 : rpar => rpar_ary
194 7 : ipar => ipar_ary
195 7 : y => y_ary
196 :
197 7 : if (.not. quiet) write (*, *) 'diffusion'
198 :
199 7 : t(0) = 0d0
200 7 : t(1) = tend
201 :
202 7 : itol = 0 ! scalar tolerances
203 14 : rtol = 1d-6
204 14 : atol = 1d-6
205 7 : h0 = atol(1)*1d-1 ! initial step size
206 :
207 7 : matrix_type_spec = banded_matrix_type
208 7 : mljac = diff_mljac
209 7 : mujac = diff_mujac
210 :
211 7 : imas = 0
212 7 : mlmas = 0
213 7 : mumas = 0
214 :
215 7 : m1 = 0
216 7 : m2 = 0
217 :
218 7 : k = nz/2
219 175 : y(1:k) = 0
220 182 : y(k + 1:nz) = ystart
221 :
222 350 : y0 = y
223 :
224 343 : do k = 1, nz
225 343 : if (k == 1) then
226 7 : sig_dm(k) = 0
227 329 : else if (5*k <= n) then
228 56 : sig_dm(k) = 0
229 273 : else if (5*k >= 4*n) then
230 70 : sig_dm(k) = 0
231 : else
232 203 : sig_dm(k) = sig
233 : end if
234 : end do
235 :
236 7 : nstep = 0
237 :
238 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
239 : diffusion_derivs, diffusion_jacob, diffusion_sjac, diffusion_solout, iout, &
240 : null_fcn_blk_dble, null_jac_blk_dble, &
241 : caller_id, nvar_blk_dble, nz_blk_dble, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
242 : n, ndisc, mljac, mujac, matrix_type_spec, null_mas, imas, mlmas, mumas, m1, m2, &
243 7 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
244 7 : if (ierr /= 0) then
245 0 : write (*, *) 'test_diffusion ierr', ierr
246 0 : call mesa_error(__FILE__, __LINE__)
247 : end if
248 :
249 : !call write_diffusion_results
250 :
251 7 : call set_yexact
252 7 : i = 0
253 7 : do k = 10, nz - 10, (nz - 10)/12
254 70 : i = i + 1
255 : if (i > n_soln) exit
256 70 : soln(i) = yexact(k)
257 70 : result(i) = y(k)
258 : end do
259 :
260 7 : call show_results(n_soln, result, soln, show_all)
261 14 : call show_statistics(ipar(i_nfcn), ipar(i_njac), nstep, show_all)
262 :
263 : contains
264 :
265 7 : subroutine set_yexact
266 : ! for nz=48, sig = 1d-2, ystart = 1d0, tend = 1d4
267 7 : yexact(1) = 0.0000000000000000D+00
268 7 : yexact(2) = 0.0000000000000000D+00
269 7 : yexact(3) = 0.0000000000000000D+00
270 7 : yexact(4) = 0.0000000000000000D+00
271 7 : yexact(5) = 0.0000000000000000D+00
272 7 : yexact(6) = 0.0000000000000000D+00
273 7 : yexact(7) = 0.0000000000000000D+00
274 7 : yexact(8) = 0.0000000000000000D+00
275 7 : yexact(9) = 2.5602881146875722D-01
276 7 : yexact(10) = 2.5830832258765279D-01
277 7 : yexact(11) = 2.6284365888436573D-01
278 7 : yexact(12) = 2.6958764648712957D-01
279 7 : yexact(13) = 2.7847002063795223D-01
280 7 : yexact(14) = 2.8939802361837114D-01
281 7 : yexact(15) = 3.0225720584222654D-01
282 7 : yexact(16) = 3.1691243229910382D-01
283 7 : yexact(17) = 3.3320909577496399D-01
284 7 : yexact(18) = 3.5097453679942375D-01
285 7 : yexact(19) = 3.7001966806159553D-01
286 7 : yexact(20) = 3.9014079814263225D-01
287 7 : yexact(21) = 4.1112164592868455D-01
288 7 : yexact(22) = 4.3273553313243934D-01
289 7 : yexact(23) = 4.5474773813802244D-01
290 7 : yexact(24) = 4.7691799008736224D-01
291 7 : yexact(25) = 4.9900307794850141D-01
292 7 : yexact(26) = 5.2075954544461744D-01
293 7 : yexact(27) = 5.4194643935567555D-01
294 7 : yexact(28) = 5.6232807598367041D-01
295 7 : yexact(29) = 5.8167678861286254D-01
296 7 : yexact(30) = 5.9977561767405640D-01
297 7 : yexact(31) = 6.1642090507187908D-01
298 7 : yexact(32) = 6.3142475475257254D-01
299 7 : yexact(33) = 6.4461732303943897D-01
300 7 : yexact(34) = 6.5584890447885325D-01
301 7 : yexact(35) = 6.6499178183713092D-01
302 7 : yexact(36) = 6.7194181237134964D-01
303 7 : yexact(37) = 6.7661972646501933D-01
304 7 : yexact(38) = 6.7897211907369082D-01
305 7 : yexact(39) = 1.0000000000000000D+00
306 7 : yexact(40) = 1.0000000000000000D+00
307 7 : yexact(41) = 1.0000000000000000D+00
308 7 : yexact(42) = 1.0000000000000000D+00
309 7 : yexact(43) = 1.0000000000000000D+00
310 7 : yexact(44) = 1.0000000000000000D+00
311 7 : yexact(45) = 1.0000000000000000D+00
312 7 : yexact(46) = 1.0000000000000000D+00
313 7 : yexact(47) = 1.0000000000000000D+00
314 7 : yexact(48) = 1.0000000000000000D+00
315 7 : end subroutine set_yexact
316 :
317 : subroutine write_diffusion_results
318 : use utils_lib, only: mkdir
319 : use const_def
320 : character(len=100) :: filename, dir
321 : integer :: k, ierr, iounit
322 : dir = 'plot_data'
323 : call mkdir(dir)
324 : filename = trim(dir)//'/'//'diffusion.data'
325 : ierr = 0
326 : open (newunit=iounit, file=trim(filename), action='write', status='replace', iostat=ierr)
327 : if (ierr == 0) then
328 : write (*, *) 'write burn results to '//trim(filename)
329 : write (iounit, '(a)') 'burn log'
330 : write (iounit, '(99(a,1x))') 'y', 'y0', 'sigma'
331 : do k = 1, nz
332 : write (iounit, '(99e24.10)') y(k), y0(k), sig_dm(k)
333 : end do
334 : close (iounit)
335 : else
336 : write (*, *) 'failed to open internals file '//trim(filename)
337 : end if
338 : end subroutine write_diffusion_results
339 :
340 : end subroutine do_test_diffusion
341 :
342 : end module test_diffusion
|
