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