Line data Source code
1 : module test_beam
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_beam, only: beam_feval, beam_jeval, beam_init, beam_solut
8 :
9 : implicit none
10 :
11 : contains
12 :
13 0 : subroutine beam_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 : integer, intent(out) :: ierr
21 0 : ierr = 0
22 0 : ipar(i_nfcn) = ipar(i_nfcn) + 1
23 0 : call beam_feval(n, x, vars, dvars_dx, ierr, rpar, ipar)
24 0 : end subroutine beam_derivs
25 :
26 0 : subroutine beam_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) :: dfdy(:, :) !(ld_dfdy,n)
31 : real(dp), intent(inout) :: f(:) !(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 beam_jeval(ld_dfdy, n, x, y, yprime, dfdy, ierr, rpar, ipar)
39 0 : if (ierr == 0) call beam_derivs(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
40 0 : end subroutine beam_jacob
41 :
42 0 : subroutine beam_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 beam_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 beam_sjac
68 :
69 0 : subroutine beam_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 : irtrn = 0
97 0 : end subroutine beam_solout
98 :
99 0 : subroutine do_test_beam(which_solver, which_decsol, numerical_jacobian, show_all, quiet)
100 : use test_support, only: show_results, show_statistics, check_results
101 : use test_int_support, only: do_test_stiff_int
102 : integer, intent(in) :: which_solver, which_decsol
103 : logical, intent(in) :: numerical_jacobian, show_all, quiet
104 :
105 : integer, parameter :: n = 80 ! the number of variables in the "beam" system of ODEs
106 0 : real(dp), target :: y_ary(n), yprime(n), yexact(n)
107 : real(dp), pointer :: y(:)
108 : integer, parameter :: lrpar = 1, lipar = 3, iout = 1
109 : logical :: consis
110 : integer, parameter :: ndisc = 0, n_soln = 9
111 0 : real(dp) :: result(n_soln), soln(n_soln)
112 0 : real(dp) :: h0, t(0:ndisc + 1), atol(1), rtol(1)
113 : integer :: i, mujac, mljac, matrix_type_spec, ierr, indsol(n_soln), imas, mlmas, mumas, m1, m2, itol, nstep
114 0 : real(dp), target :: rpar_ary(lrpar)
115 : integer, target :: ipar_ary(lipar)
116 : real(dp), pointer :: rpar(:)
117 : integer, pointer :: ipar(:)
118 : integer :: caller_id, nvar, nz
119 0 : real(dp), dimension(:), pointer :: lblk, dblk, ublk ! =(nvar,nvar,nz)
120 0 : real(dp), dimension(:), pointer :: uf_lblk, uf_dblk, uf_ublk ! =(nvar,nvar,nz)
121 :
122 0 : nullify (lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk)
123 0 : caller_id = 0
124 0 : nvar = 0
125 0 : nz = 0
126 :
127 0 : rpar => rpar_ary
128 0 : ipar => ipar_ary
129 0 : y => y_ary
130 :
131 0 : if (.not. quiet) write (*, *) 'beam'
132 :
133 0 : t(0) = 0
134 0 : t(1) = 5d0
135 :
136 0 : itol = 0 ! scalar tolerances
137 0 : rtol(1) = 1d-3
138 0 : atol(1) = 1d-3
139 0 : h0 = 1d-4 ! initial step size
140 :
141 0 : m1 = n/2
142 0 : m2 = 0
143 :
144 0 : mljac = n
145 0 : mujac = n
146 0 : matrix_type_spec = square_matrix_type
147 :
148 0 : imas = 0
149 0 : mlmas = 0
150 0 : mumas = 0
151 :
152 0 : if (.not. numerical_jacobian) then
153 0 : write (*, *) 'beam test only supports numerical jacobian'
154 0 : return
155 : end if
156 :
157 0 : call beam_init(n, y, yprime, consis)
158 0 : nstep = 0
159 : call do_test_stiff_int(which_solver, which_decsol, numerical_jacobian, &
160 : beam_derivs, beam_jacob, beam_sjac, beam_solout, iout, &
161 : null_fcn_blk_dble, null_jac_blk_dble, &
162 : caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
163 : n, ndisc, mljac, mujac, matrix_type_spec, null_mas, imas, mlmas, mumas, m1, m2, &
164 0 : t, rtol, atol, itol, h0, y, nstep, lrpar, rpar, lipar, ipar, quiet, ierr)
165 0 : if (ierr /= 0) then
166 0 : write (*, *) 'test_beam ierr', ierr
167 0 : call mesa_error(__FILE__, __LINE__)
168 : end if
169 :
170 0 : call beam_solut(n, 0d0, yexact)
171 0 : indsol(1:n_soln) = [1, 10, 20, 30, 40, 50, 60, 70, 80]
172 0 : do i = 1, n_soln
173 0 : result(i) = y(indsol(i))
174 0 : soln(i) = yexact(indsol(i))
175 : end do
176 :
177 0 : call check_results(n, y, yexact, rtol(1)*1d2, ierr)
178 0 : if (ierr /= 0) then
179 0 : write (*, *) 'check results ierr', ierr
180 0 : call mesa_error(__FILE__, __LINE__) ! do_test_vdpol
181 : end if
182 :
183 0 : if (quiet) return
184 :
185 0 : call show_results(n_soln, result, soln, show_all)
186 0 : call show_statistics(ipar(i_nfcn), ipar(i_njac), nstep, show_all)
187 :
188 0 : end subroutine do_test_beam
189 :
190 : end module test_beam
|
