Line data Source code
1 :
2 : module test_newton
3 : use const_def, only: dp
4 : use num_def
5 : use num_lib
6 : use mtx_def
7 : use mtx_lib
8 : use math_lib
9 : use utils_lib, only: mesa_error
10 :
11 : implicit none
12 :
13 : real(dp), parameter :: one = 1
14 :
15 : integer, parameter :: nz = 1001, nvar = 2 !use odd number of zones for problem symmetry
16 : integer, parameter :: nsec = 0 ! number of secondaries per zone
17 : integer, parameter :: ldy = nz
18 :
19 : integer, parameter :: i_conc = 1, i_flux = 2, equ_conc = 1, equ_flux = 2
20 :
21 : integer :: matrix_type
22 : real(dp), pointer, dimension(:) :: equ1, x1, xold1, dx1, xscale1, y1
23 : real(dp), pointer, dimension(:, :) :: equ, x, xold, dx, xscale, y
24 : real(dp), pointer, dimension(:, :, :) :: ublk, dblk, lblk
25 :
26 : logical, parameter :: dbg = .false.
27 :
28 : contains
29 :
30 2 : subroutine do_test_newton( &
31 : do_numerical_jacobian, which_decsol_in)
32 : logical, intent(in) :: do_numerical_jacobian
33 : integer, intent(in) :: which_decsol_in
34 :
35 2 : real(dp) :: dt, kappa, alphat, tmax, time, actual
36 2 : real(dp), pointer, dimension(:) :: concentration, fluxes
37 2 : real(dp) :: xmin, xmax, delx
38 : integer :: ierr, which_decsol, numsteps, midpt, maxsteps, neq
39 : character(len=64) :: decsol_option_name
40 :
41 : real(dp), parameter :: expected = 2.9347118120566711D-02 ! using lapack
42 :
43 : include 'formats'
44 :
45 2 : which_decsol = which_decsol_in
46 2 : call decsol_option_str(which_decsol, decsol_option_name, ierr)
47 2 : if (ierr /= 0) return
48 :
49 1 : write (*, *) 'do_test_newton using '//trim(decsol_option_name)
50 :
51 : ! diffusion problem setup
52 1 : kappa = 1.0
53 1 : alphat = 1.d1 ! multiply explicit stability time step by this factor
54 1 : time = 0.0
55 1 : xmin = 0.0
56 1 : xmax = 1.0
57 1 : delx = (xmax - xmin)/float(nz) !use uniform spatial mesh
58 1 : tmax = pow2(10.0*delx)/kappa !maximum evolution time in units of stability time step
59 :
60 1 : allocate (concentration(nz), fluxes(nz), stat=ierr)
61 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
62 :
63 1 : neq = nvar*nz
64 : allocate ( &
65 : equ1(neq), x1(neq), xold1(neq), dx1(neq), &
66 1 : xscale1(neq), y1(ldy*nsec), stat=ierr)
67 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
68 :
69 1 : x(1:nvar, 1:nz) => x1(1:neq)
70 1 : xold(1:nvar, 1:nz) => xold1(1:neq)
71 1 : dx(1:nvar, 1:nz) => dx1(1:neq)
72 1 : equ(1:nvar, 1:nz) => equ1(1:neq)
73 1 : xscale(1:nvar, 1:nz) => xscale1(1:neq)
74 1 : y(1:ldy, 1:nsec) => y1(1:ldy*nsec)
75 :
76 1002 : concentration = 0.0
77 1002 : fluxes = 0.0
78 1 : midpt = ceiling(float(nz)/2)
79 1 : concentration(midpt) = 1.0 !delta function spike
80 1 : numsteps = 0
81 1 : maxsteps = 500
82 1 : dt = alphat*(delx*delx)/kappa ! explicit stability time step multiplied by alphat
83 11 : do while (time < tmax .and. numsteps < maxsteps)
84 10 : numsteps = numsteps + 1
85 : call do_1step_diffuse( &
86 : do_numerical_jacobian, which_decsol, &
87 10 : dt, kappa, nz, nvar, concentration, fluxes, delx, ierr)
88 10 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
89 10 : time = time + dt
90 : !write(*,2) 'diffusion step', numsteps, concentration(midpt), time/tmax
91 : end do
92 :
93 1 : actual = concentration(midpt)
94 1 : write (*, 1) 'expected', expected
95 1 : write (*, 1) 'actual', actual
96 : !write(*,1) '(actual - expected)/expected', (actual - expected)/expected
97 1 : write (*, *)
98 :
99 1 : deallocate (concentration, fluxes)
100 1 : deallocate (equ1, x1, xold1, dx1, xscale1, y1)
101 :
102 2 : end subroutine do_test_newton
103 :
104 10 : subroutine do_1step_diffuse( &
105 : do_numerical_jacobian, which_decsol, &
106 : dt, kappa, nz, nvar, concentration, fluxes, delx, ierr)
107 : logical, intent(in) :: do_numerical_jacobian
108 : integer, intent(in) :: which_decsol
109 : integer, intent(in) :: nz, nvar
110 : real(dp), intent(inout) :: dt, kappa
111 : real(dp), pointer, dimension(:), intent(inout) :: concentration, fluxes
112 : real(dp), intent(in) :: delx
113 : integer, intent(out) :: ierr
114 :
115 : integer :: liwork, lwork
116 10 : integer, dimension(:), pointer :: iwork
117 10 : real(dp), dimension(:), pointer :: work
118 :
119 : integer, parameter :: lrpar = 3, lipar = 1
120 : integer, target :: ipar_target(lipar)
121 40 : real(dp), target :: rpar_target(lrpar)
122 : integer, pointer :: ipar(:)
123 : real(dp), pointer :: rpar(:)
124 :
125 : integer :: lrd, lid
126 10 : integer, pointer :: ipar_decsol(:) ! (lid)
127 10 : real(dp), pointer :: rpar_decsol(:) ! (lrd)
128 :
129 : integer :: mljac, mujac
130 10 : real(dp) :: tol_correction_norm, tol_max_correction, tol_residual_norm
131 : logical :: nonconv
132 :
133 : include 'formats'
134 :
135 10 : ierr = 0
136 :
137 10 : ipar => ipar_target
138 10 : rpar => rpar_target
139 :
140 10 : rpar(1) = dt
141 10 : rpar(2) = kappa
142 10 : rpar(3) = delx
143 :
144 10 : if (do_numerical_jacobian) then
145 0 : ipar(1) = 1
146 : else
147 10 : ipar(1) = 0
148 : end if
149 :
150 10 : call set_fluxes(concentration, fluxes, kappa, delx)
151 20040 : xold(i_conc, :) = concentration ! starting model
152 20040 : xold(i_flux, :) = fluxes
153 30040 : dx = 0d0
154 60080 : x = xold
155 :
156 10 : tol_correction_norm = 1d-9 ! upper limit on magnitude of average scaled correction
157 10 : tol_max_correction = 1d99
158 10 : tol_residual_norm = 1d99
159 :
160 10 : mljac = 2*nvar - 1 ! number of subdiagonals
161 10 : mujac = mljac ! number of superdiagonals
162 10 : if (which_decsol == lapack) then
163 10 : call lapack_work_sizes(nz*nvar, lrd, lid)
164 10 : if (do_numerical_jacobian) then
165 0 : matrix_type = square_matrix_type
166 0 : mljac = nz*nvar - 1
167 0 : mujac = mljac
168 0 : if (nz*nvar > 51) then
169 0 : write (*, *) 'numerical jac is very slow for large nz*nvar'
170 0 : call mesa_error(__FILE__, __LINE__)
171 : end if
172 : else
173 10 : matrix_type = banded_matrix_type
174 : end if
175 : else
176 0 : write (*, *) 'bad value for which_decsol'
177 0 : call mesa_error(__FILE__, __LINE__)
178 : end if
179 :
180 10 : allocate (rpar_decsol(lrd), ipar_decsol(lid), stat=ierr)
181 10 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
182 :
183 : call newton_work_sizes( &
184 10 : mljac, mujac, nvar, nz, nsec, matrix_type, lwork, liwork, ierr)
185 10 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
186 :
187 10 : allocate (work(lwork), iwork(liwork), stat=ierr)
188 10 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
189 :
190 320600 : work = 0
191 20220 : iwork = 0
192 :
193 10 : iwork(i_try_really_hard) = 1 ! try really hard for first model
194 10 : iwork(i_model_number) = 1
195 :
196 : !iwork(i_debug) = 1
197 :
198 10 : if (which_decsol == lapack) then
199 : call do_newt( &
200 : lapack_decsol, null_decsolblk, &
201 10 : nonconv, ierr)
202 : else
203 0 : stop 'bad which_decsol'
204 : end if
205 :
206 10 : if (nonconv) then
207 0 : write (*, *) 'failed to converge'
208 0 : call mesa_error(__FILE__, __LINE__)
209 : end if
210 :
211 20040 : concentration = x(i_conc, :)
212 20040 : fluxes = x(i_flux, :)
213 :
214 10 : deallocate (iwork, work, rpar_decsol, ipar_decsol)
215 :
216 : contains
217 :
218 10 : subroutine do_newt(decsol, decsolblk, nonconv, ierr)
219 : interface
220 : include "mtx_decsol.dek"
221 : include "mtx_decsolblk_dble.dek"
222 : end interface
223 : logical, intent(out) :: nonconv
224 : integer, intent(out) :: ierr
225 :
226 10 : real(dp), pointer :: AF1(:)
227 10 : AF1 => null()
228 : call newton( &
229 : nz, nvar, x1, xold1, matrix_type, mljac, mujac, &
230 : decsol, decsolblk, &
231 : lrd, rpar_decsol, lid, ipar_decsol, which_decsol, tol_correction_norm, &
232 : default_set_primaries, default_set_secondaries, diffusion_set_xscale, &
233 : default_Bdomain, default_xdomain, eval_equations, &
234 : default_size_equ, default_sizeB, default_inspectB, &
235 : enter_setmatrix, exit_setmatrix, failed_in_setmatrix, &
236 : default_force_another_iter, xscale1, equ1, ldy, nsec, y1, &
237 : work, lwork, iwork, liwork, &
238 10 : AF1, lrpar, rpar, lipar, ipar, nonconv, ierr)
239 10 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
240 10 : if (associated(AF1)) deallocate (AF1)
241 10 : end subroutine do_newt
242 :
243 : end subroutine do_1step_diffuse
244 :
245 10 : subroutine set_fluxes(concentration, fluxes, kappa, delx)
246 : real(dp), intent(in) :: delx, kappa
247 : real(dp), pointer, dimension(:), intent(inout) :: concentration, fluxes
248 : integer :: i
249 10010 : do i = 2, nz
250 10010 : fluxes(i) = -kappa*(concentration(i) - concentration(i - 1))/delx
251 : end do
252 10 : fluxes(1) = 0
253 10 : end subroutine set_fluxes
254 :
255 61 : subroutine eval_equations(iter, nvar, nz, x, xscale, equ, lrpar, rpar, lipar, ipar, ierr)
256 : integer, intent(in) :: iter, nvar, nz
257 : real(dp), pointer, dimension(:, :) :: x, xscale, equ ! (nvar, nz)
258 : integer, intent(in) :: lrpar, lipar
259 : real(dp), intent(inout) :: rpar(:) ! (lrpar)
260 : integer, intent(inout) :: ipar(:) ! (lipar)
261 : integer, intent(out) :: ierr
262 :
263 : integer :: k
264 61 : real(dp) :: dt, kappa, delx, fnzp1
265 61 : ierr = 0
266 :
267 61 : dt = rpar(1)
268 61 : kappa = rpar(2)
269 61 : delx = rpar(3)
270 :
271 183244 : equ = 0.0
272 :
273 : ! Equation 1 : dc/dt = -div(fluxes)
274 61061 : do k = 1, nz - 1
275 61000 : equ(equ_conc, k) = &
276 122061 : x(i_conc, k) - xold(i_conc, k) + dt*(x(i_flux, k + 1) - x(i_flux, k))/delx
277 : end do
278 61 : fnzp1 = kappa*x(i_conc, nz)/delx
279 61 : equ(equ_conc, nz) = &
280 61 : x(i_conc, nz) - xold(i_conc, nz) + dt*(fnzp1 - x(i_flux, k))/delx
281 :
282 : ! Equation 2 : flux = -kappa*gradient(c)
283 61061 : do k = 2, nz
284 61061 : equ(equ_flux, k) = x(i_flux, k) + kappa*(x(i_conc, k) - x(i_conc, k - 1))/delx
285 : end do
286 61 : equ(equ_flux, 1) = x(i_flux, 1)
287 :
288 61 : end subroutine eval_equations
289 :
290 20 : subroutine eval_jacobian(ldA, A1, idiag, lrpar, rpar, lipar, ipar, ierr)
291 : integer, intent(in) :: ldA ! leading dimension of A
292 : real(dp), pointer :: A1(:) ! (ldA, nvar*nz) ! the jacobian matrix
293 : ! A(idiag+q-v, v) = partial of equation(q) wrt variable(v)
294 : integer, intent(inout) :: idiag
295 : integer, intent(in) :: lrpar, lipar
296 : real(dp), intent(inout) :: rpar(:) ! (lrpar)
297 : integer, intent(inout) :: ipar(:) ! (lipar)
298 : integer, intent(out) :: ierr
299 :
300 : integer :: k
301 20 : real(dp) :: dt, kappa, delx
302 20 : real(dp), pointer :: blk3(:, :, :, :)
303 20 : ierr = 0
304 :
305 20 : blk3(1:nvar, 1:nvar, 1:nz, 1:3) => A1(1:nvar*nvar*nz*3)
306 20 : ublk => blk3(:, :, :, 1)
307 20 : dblk => blk3(:, :, :, 2)
308 20 : lblk => blk3(:, :, :, 3)
309 :
310 20 : dt = rpar(1)
311 20 : kappa = rpar(2)
312 20 : delx = rpar(3)
313 :
314 280300 : A1 = 0.0
315 20020 : do k = 1, nz - 1 ! concentration equ
316 20000 : call e00(A1, equ_conc, i_conc, k, idiag, ldA, one)
317 20000 : call e00(A1, equ_conc, i_flux, k, idiag, ldA, -dt/delx)
318 20020 : call ep1(A1, equ_conc, i_flux, k, idiag, ldA, dt/delx)
319 : end do
320 20 : call e00(A1, equ_conc, i_conc, nz, idiag, ldA, one + kappa*dt/delx**2)
321 :
322 20020 : do k = 2, nz ! flux equ
323 20000 : call e00(A1, equ_flux, i_flux, k, idiag, ldA, one)
324 20000 : call e00(A1, equ_flux, i_conc, k, idiag, ldA, kappa/delx)
325 20020 : call em1(A1, equ_flux, i_conc, k, idiag, ldA, -kappa/delx)
326 : end do
327 20 : call e00(A1, equ_flux, i_flux, 1, idiag, ldA, one)
328 :
329 20 : end subroutine eval_jacobian
330 :
331 80040 : subroutine e00(A1, i, j, k, idiag, ldA, v) ! partial of equ(i,k) wrt var(j,k)
332 : real(dp), pointer :: A1(:)
333 : real(dp) :: v
334 80040 : real(dp), pointer :: A(:, :)
335 : integer, intent(in) :: i, j, k, idiag, ldA
336 : integer :: b, q, v00, sz, neq
337 80040 : sz = size(A1, dim=1)
338 80040 : neq = sz/ldA
339 80040 : A(1:ldA, 1:neq) => A1(1:sz)
340 80040 : if (matrix_type == square_matrix_type) then
341 0 : b = nvar*(k - 1)
342 0 : A(b + i, b + j) = A(b + i, b + j) + v
343 80040 : else if (matrix_type == banded_matrix_type) then
344 80040 : b = nvar*(k - 1)
345 80040 : q = idiag + b + i
346 80040 : v00 = b + j
347 80040 : A(q - v00, v00) = A(q - v00, v00) + v
348 0 : else if (matrix_type == block_tridiag_dble_matrix_type) then
349 0 : dblk(i, j, k) = dblk(i, j, k) + v
350 : else
351 0 : stop 'bad matrix_type'
352 : end if
353 80040 : end subroutine e00
354 :
355 20000 : subroutine em1(A1, i, j, k, idiag, ldA, v) ! partial of equ(i,k) wrt var(j,k-1)
356 : real(dp), pointer :: A1(:)
357 : real(dp) :: v
358 20000 : real(dp), pointer :: A(:, :)
359 : integer, intent(in) :: i, j, k, idiag, ldA
360 : integer :: b, q, vm1, sz, neq
361 20000 : sz = size(A1, dim=1)
362 20000 : neq = sz/ldA
363 20000 : A(1:ldA, 1:neq) => A1(1:sz)
364 20000 : if (k <= 0) stop 'bad k for em1'
365 20000 : if (matrix_type == square_matrix_type) then
366 0 : b = nvar*(k - 1)
367 0 : A(b + i, b + j - nvar) = A(b + i, b + j - nvar) + v
368 20000 : else if (matrix_type == banded_matrix_type) then
369 20000 : b = nvar*(k - 1)
370 20000 : q = idiag + b + i
371 20000 : vm1 = b + j - nvar
372 20000 : A(q - vm1, vm1) = A(q - vm1, vm1) + v
373 0 : else if (matrix_type == block_tridiag_dble_matrix_type) then
374 0 : lblk(i, j, k) = lblk(i, j, k) + v
375 : else
376 0 : stop 'bad matrix_type'
377 : end if
378 20000 : end subroutine em1
379 :
380 20000 : subroutine ep1(A1, i, j, k, idiag, ldA, v) ! partial of equ(i,k) wrt var(j,k+1)
381 : real(dp), pointer :: A1(:)
382 : real(dp) :: v
383 20000 : real(dp), pointer :: A(:, :)
384 : integer, intent(in) :: i, j, k, idiag, ldA
385 : integer :: b, q, vp1, sz, neq
386 20000 : sz = size(A1, dim=1)
387 20000 : neq = sz/ldA
388 20000 : A(1:ldA, 1:neq) => A1(1:sz)
389 20000 : if (k >= nz) stop 'bad k for ep1'
390 20000 : if (matrix_type == square_matrix_type) then
391 0 : b = nvar*(k - 1)
392 0 : A(b + i, b + j + nvar) = A(b + i, b + j + nvar) + v
393 20000 : else if (matrix_type == banded_matrix_type) then
394 20000 : b = nvar*(k - 1)
395 20000 : q = idiag + b + i
396 20000 : vp1 = b + j + nvar
397 20000 : A(q - vp1, vp1) = A(q - vp1, vp1) + v
398 0 : else if (matrix_type == block_tridiag_dble_matrix_type) then
399 0 : ublk(i, j, k) = ublk(i, j, k) + v
400 : else
401 0 : stop 'bad matrix_type'
402 : end if
403 20000 : end subroutine ep1
404 :
405 20 : subroutine enter_setmatrix( &
406 : iter, nvar, nz, neqs, x, xold, xscale, xder, need_solver_to_eval_jacobian, &
407 20 : ldA, A1, idiag, lrpar, rpar, lipar, ipar, ierr)
408 : integer, intent(in) :: iter, nvar, nz, neqs
409 : real(dp), pointer, dimension(:, :) :: x, xold, xscale, xder ! (nvar, nz)
410 : logical, intent(out) :: need_solver_to_eval_jacobian
411 : integer, intent(in) :: ldA ! leading dimension of A
412 : real(dp), pointer, dimension(:) :: A1 ! =(ldA, neqs)
413 : integer, intent(inout) :: idiag
414 : integer, intent(in) :: lrpar, lipar
415 : real(dp), intent(inout) :: rpar(:) ! (lrpar)
416 : integer, intent(inout) :: ipar(:) ! (lipar)
417 : integer, intent(out) :: ierr
418 20 : real(dp) :: epsder
419 : include 'formats'
420 : if (dbg) write (*, '(/, a)') 'enter_setmatrix'
421 20 : if (ipar(1) /= 0) then ! do numerical jacobian
422 0 : epsder = 1d-6 ! relative variation to compute numerical derivatives
423 0 : xder = epsder*(xscale + abs(xold))
424 0 : need_solver_to_eval_jacobian = .true.
425 : else
426 20 : call eval_jacobian(ldA, A1, idiag, lrpar, rpar, lipar, ipar, ierr)
427 20 : need_solver_to_eval_jacobian = .false.
428 : end if
429 20 : end subroutine enter_setmatrix
430 :
431 0 : subroutine exit_setmatrix( &
432 : iter, nvar, nz, neqs, dx, ldA, A1, idiag, xscale, lrpar, rpar, lipar, ipar, ierr)
433 : integer, intent(in) :: ldA ! leading dimension of A
434 : integer, intent(in) :: iter, nvar, nz, neqs ! number of equations, 2nd dimension of A
435 : integer, intent(inout) :: idiag ! row of A with the matrix diagonal entries
436 : real(dp), pointer, dimension(:, :) :: dx
437 : real(dp), pointer, dimension(:) :: A1
438 : real(dp), pointer, dimension(:, :) :: xscale ! (nvar, nz)
439 : integer, intent(in) :: lrpar, lipar
440 : real(dp), intent(inout) :: rpar(:) ! (lrpar)
441 : integer, intent(inout) :: ipar(:) ! (lipar)
442 : integer, intent(out) :: ierr
443 : if (dbg) write (*, '(a, /)') 'exit_setmatrix'
444 0 : ierr = 0
445 0 : end subroutine exit_setmatrix
446 :
447 0 : subroutine failed_in_setmatrix(j, lrpar, rpar, lipar, ipar, ierr)
448 : integer, intent(in) :: j
449 : integer, intent(in) :: lrpar, lipar
450 : real(dp), intent(inout) :: rpar(:) ! (lrpar)
451 : integer, intent(inout) :: ipar(:) ! (lipar)
452 : integer, intent(out) :: ierr
453 : if (dbg) write (*, '(a, /)') 'failed_in_setmatrix'
454 0 : ierr = 0
455 0 : end subroutine failed_in_setmatrix
456 :
457 : ! you might want to use a different value of xscale_min for this
458 10 : subroutine diffusion_set_xscale(nvar, nz, xold, xscale, lrpar, rpar, lipar, ipar, ierr)
459 : integer, intent(in) :: nvar, nz
460 : real(dp), pointer :: xold(:, :) ! (nvar, nz)
461 : real(dp), pointer :: xscale(:, :) ! (nvar, nz)
462 : integer, intent(in) :: lrpar, lipar
463 : real(dp), intent(inout) :: rpar(:) ! (lrpar)
464 : integer, intent(inout) :: ipar(:) ! (lipar)
465 : integer, intent(out) :: ierr
466 : ! real(dp), parameter :: xscale_min = 1d0
467 30040 : xscale = 1.d0 ! max(xscale_min, abs(xold))
468 10 : ierr = 0
469 10 : end subroutine diffusion_set_xscale
470 :
471 : end module test_newton
|
