Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011 Bill Paxton & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module test_block_tri_dble
21 :
22 : use mtx_lib
23 : use mtx_def
24 : use const_def, only: dp
25 : use utils_lib, only: is_bad, mesa_error
26 :
27 : implicit none
28 :
29 : integer, parameter :: fltp = dp
30 : integer, parameter :: caller_id = 0
31 :
32 : contains
33 :
34 1 : subroutine do_test_block_tri_dble
35 1 : call test_block(bcyclic_dble, .true.)
36 1 : call test_block(block_thomas_dble, .true.)
37 1 : end subroutine do_test_block_tri_dble
38 :
39 2 : subroutine test_block(which_decsol_option, for_release)
40 :
41 : integer, intent(in) :: which_decsol_option
42 : logical, intent(in) :: for_release
43 :
44 2 : real(fltp), pointer :: lblk1(:), dblk1(:), ublk1(:) ! (nvar,nvar,nz)
45 2 : real(fltp), pointer :: lblk(:, :, :), dblk(:, :, :), ublk(:, :, :) ! (nvar,nvar,nz)
46 2 : real(fltp), pointer :: x(:, :), xcorrect(:, :), brhs(:, :), work(:, :) ! (nvar,nz)
47 2 : real(fltp), pointer :: x1(:) ! =(nvar,nz)
48 2 : integer, pointer :: ipiv1(:) ! =(nvar,nz)
49 2 : integer, pointer :: ipiv(:, :) ! (nvar,nz)
50 :
51 2 : real(dp), pointer :: rpar_decsol(:) ! (lrd)
52 2 : integer, pointer :: ipar_decsol(:) ! (lid)
53 2 : real(fltp) :: time_factor, time_solve, time_refine, time_dealloc
54 :
55 : integer :: ierr, lid, lrd, nvar, nz
56 : character(len=255) :: fname, which_decsol_str
57 :
58 : include 'formats'
59 :
60 2 : ierr = 0
61 :
62 2 : call decsol_option_str(which_decsol_option, which_decsol_str, ierr)
63 2 : if (ierr /= 0) return
64 :
65 1 : if (for_release) then
66 1 : fname = 'block_tri.data'
67 : else
68 0 : fname = 'block_tri_12.data'
69 : end if
70 :
71 1 : time_factor = 0; time_solve = 0; time_refine = 0; time_dealloc = 0
72 :
73 1 : call read_testfile(fname)
74 :
75 1 : if (which_decsol_option == bcyclic_dble) then
76 1 : write (*, *) 'bcyclic_dble'
77 1 : call bcyclic_dble_work_sizes(nvar, nz, lrd, lid)
78 : else
79 0 : write (*, *) 'bad value for which_decsol_option in test_block'
80 0 : call mesa_error(__FILE__, __LINE__)
81 : end if
82 :
83 : allocate ( &
84 : rpar_decsol(lrd), ipar_decsol(lid), x1(nvar*nz), xcorrect(nvar, nz), &
85 1 : brhs(nvar, nz), ipiv1(nvar*nz), work(nvar, nz), stat=ierr)
86 1 : if (ierr /= 0) then
87 0 : write (*, *) 'failed in alloc'
88 0 : call mesa_error(__FILE__, __LINE__)
89 : end if
90 1 : ipiv(1:nvar, 1:nz) => ipiv1(1:nvar*nz)
91 1 : x(1:nvar, 1:nz) => x1(1:nvar*nz)
92 :
93 1 : call set_xcorrect
94 1 : call set_brhs(lblk, dblk, ublk)
95 :
96 1 : if (which_decsol_option == bcyclic_dble) then
97 1 : call solve_blocks(bcyclic_dble_decsolblk)
98 : else
99 0 : write (*, *) 'missing case for which_decsol_option', which_decsol_option
100 0 : call mesa_error(__FILE__, __LINE__)
101 : end if
102 :
103 1 : call check_x
104 :
105 1 : if (which_decsol_option == bcyclic_dble) then
106 1 : write (*, *) 'done bcyclic_dble'
107 0 : else if (which_decsol_option == block_thomas_dble) then
108 0 : write (*, *) 'done block_thomas_dble'
109 : end if
110 :
111 1 : write (*, *)
112 :
113 0 : deallocate (rpar_decsol, ipar_decsol, x1, xcorrect, work, &
114 3 : brhs, ipiv1, lblk1, dblk1, ublk1)
115 :
116 : contains
117 :
118 1 : subroutine solve_blocks(decsolblk)
119 : interface
120 : include 'mtx_decsolblk_dble.dek'
121 : end interface
122 :
123 : integer :: iop, rep
124 : integer :: j, k
125 :
126 : include 'formats'
127 :
128 1 : iop = 0 ! factor A
129 : call decsolblk( &
130 1 : iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr)
131 1 : if (ierr /= 0) then
132 0 : write (*, *) 'decsolblk failed for factor'
133 0 : call mesa_error(__FILE__, __LINE__)
134 : end if
135 :
136 2 : do rep = 1, 1
137 :
138 1 : iop = 1 ! solve A*x = b
139 :
140 21 : do k = 1, nz
141 521 : do j = 1, nvar
142 520 : x(j, k) = brhs(j, k)
143 : end do
144 : end do
145 : call decsolblk( &
146 1 : iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr)
147 2 : if (ierr /= 0) then
148 0 : write (*, *) 'decsolblk failed for solve'
149 0 : call mesa_error(__FILE__, __LINE__)
150 : end if
151 :
152 : end do
153 :
154 1 : iop = 2 ! deallocate
155 : call decsolblk( &
156 1 : iop, caller_id, nvar, nz, lblk1, dblk1, ublk1, x1, ipiv1, lrd, rpar_decsol, lid, ipar_decsol, ierr)
157 1 : if (ierr /= 0) then
158 0 : write (*, *) 'decsolblk failed for deallocate'
159 0 : call mesa_error(__FILE__, __LINE__)
160 : end if
161 :
162 1 : end subroutine solve_blocks
163 :
164 2 : subroutine read_testfile(fname)
165 : character(len=*), intent(in) :: fname
166 : integer :: iounit, ierr
167 : !write(*,*) 'reading ' // trim(fname)
168 1 : iounit = 33; ierr = 0
169 1 : open (unit=iounit, file=trim(fname), status='old', action='read', iostat=ierr)
170 1 : if (ierr /= 0) then
171 0 : write (*, *) 'failed to open '//trim(fname)
172 0 : call mesa_error(__FILE__, __LINE__)
173 : end if
174 :
175 1 : call mtx_read_block_tridiagonal(iounit, nvar, nz, lblk1, dblk1, ublk1, ierr)
176 :
177 1 : if (ierr /= 0) then
178 0 : write (*, *) 'failed to read '//trim(fname)
179 0 : call mesa_error(__FILE__, __LINE__)
180 : end if
181 1 : close (iounit)
182 :
183 1 : lblk(1:nvar, 1:nvar, 1:nz) => lblk1(1:nvar*nvar*nz)
184 1 : dblk(1:nvar, 1:nvar, 1:nz) => dblk1(1:nvar*nvar*nz)
185 1 : ublk(1:nvar, 1:nvar, 1:nz) => ublk1(1:nvar*nvar*nz)
186 :
187 : !return
188 1 : nz = 20
189 1 : write (*, *) 'testing with nvar,nz', nvar, nz
190 :
191 1 : end subroutine read_testfile
192 :
193 1 : subroutine set_brhs(lblk, dblk, ublk)
194 : real(fltp), pointer, dimension(:, :, :) :: lblk, dblk, ublk
195 : integer :: k, j
196 : include 'formats'
197 : ! set brhs = A*xcorrect
198 :
199 1 : call block_dble_mv(nvar, nz, lblk, dblk, ublk, xcorrect, brhs)
200 :
201 1 : return
202 : do k = 1, 2 !nz
203 : do j = 1, nvar
204 : if (brhs(j, k) /= 0) write (*, 3) 'brhs xcorrect', j, k, brhs(j, k), xcorrect(j, k)
205 : end do
206 : end do
207 : write (*, *) 'end set_brhs'
208 : stop
209 : end subroutine set_brhs
210 :
211 1 : subroutine check_x
212 1 : real(fltp) :: max_err, atol, rtol, avg_err
213 : integer :: i_max, j_max, i, j
214 : include 'formats'
215 1 : atol = 1d-4
216 1 : rtol = 1d-4
217 1 : call check1_x(avg_err, max_err, atol, rtol, i_max, j_max)
218 1 : i = i_max; j = j_max
219 1 : if (max_err > 1) then
220 0 : write (*, 3) 'BAD: err, x, xcorrect', i, j, max_err, x(i, j), xcorrect(i, j)
221 : !write(*,3) 'BAD: avg err, max err, x, xcorrect', i, j, avg_err, max_err, x(i,j), xcorrect(i,j)
222 : end if
223 1 : end subroutine check_x
224 :
225 1 : subroutine check1_x(avg_err, max_err, atol, rtol, i_max, j_max)
226 : real(fltp), intent(out) :: avg_err, max_err
227 : real(fltp), intent(in) :: atol, rtol
228 : integer, intent(out) :: i_max, j_max
229 : integer :: i, j
230 1 : real(fltp) :: err_sum
231 1 : real(fltp) :: err
232 : include 'formats'
233 1 : max_err = 0; i_max = 0; j_max = 0; err_sum = 0
234 21 : do j = 1, nz
235 521 : do i = 1, nvar
236 500 : if (is_bad(x(i, j))) then
237 0 : write (*, 3) 'x xcorrect', i, j, x(i, j), xcorrect(i, j)
238 0 : stop 'check1_x'
239 : end if
240 500 : err = abs(x(i, j) - xcorrect(i, j))/(atol + rtol*max(abs(x(i, j)), abs(xcorrect(i, j))))
241 500 : err_sum = err_sum + err
242 520 : if (err > max_err) then
243 4 : max_err = err; i_max = i; j_max = j
244 : end if
245 : end do
246 : end do
247 1 : avg_err = err_sum/(nz*nvar)
248 : !write(*,1) 'avg_err', avg_err
249 1 : end subroutine check1_x
250 :
251 1 : subroutine set_xcorrect
252 1 : real(fltp) :: cnt
253 : integer :: k, j
254 1 : cnt = 1d0
255 21 : do k = 1, nz
256 521 : do j = 1, nvar
257 500 : cnt = cnt + 1d-3
258 520 : xcorrect(j, k) = cnt
259 : end do
260 : end do
261 1 : end subroutine set_xcorrect
262 :
263 : end subroutine test_block
264 : end module test_block_tri_dble
|
