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