SICOPOLIS V5-dev  Revision 1288
sico_maths_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s i c o _ m a t h s _ m
4 !
5 !> @file
6 !!
7 !! Several mathematical tools used by SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2018 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Several mathematical tools used by SICOPOLIS.
34 !<------------------------------------------------------------------------------
36 
37  use sico_types_m
38 
39  implicit none
40 
41 #ifndef ALLOW_OPENAD
42  public
43 #else
44 
45  public :: sor_sprs, tri_sle, bilinint, my_erfc
46 
47 #ifdef BUILD_LIS
48  public :: sico_lis_solver
49 #endif
50 
51  interface sor_sprs
52  module procedure sor_sprs_local
53  end interface
54 
55  interface tri_sle
56  module procedure tri_sle_local
57  end interface
58 
59 #ifdef BUILD_LIS
60  interface sico_lis_solver
61  module procedure sico_lis_solver_local
62  end interface
63 #endif
64 
65 #endif
66 
67 contains
68 
69 !-------------------------------------------------------------------------------
70 !> SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b
71 !! [matrix storage: compressed sparse row CSR,
72 !! represented by arrays lgs_a_value(values), lgs_a_index (indices)
73 !! and lgs_a_ptr (pointers)].
74 !<------------------------------------------------------------------------------
75 #ifndef ALLOW_OPENAD
76  subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, &
77 #else
78  subroutine sor_sprs_local(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, &
79 #endif
80  lgs_b_value, &
81  nnz, nmax, &
82 #ifdef ALLOW_OPENAD
83  n_sprs, &
84 #endif
85  omega, eps_sor, lgs_x_value, ierr)
86 
87  implicit none
88 
89 #ifdef ALLOW_OPENAD
90  integer(i4b), intent(in) :: n_sprs
91 #endif
92  integer(i4b), intent(in) :: nnz, nmax
93  real(dp), intent(in) :: omega, eps_sor
94  integer(i4b), dimension(nmax+1), intent(in) :: lgs_a_ptr
95  integer(i4b), dimension(nnz), intent(in) :: lgs_a_index
96  integer(i4b), dimension(nmax), intent(in) :: lgs_a_diag_index
97  real(dp), dimension(nnz), intent(in) :: lgs_a_value
98  real(dp), dimension(nmax), intent(in) :: lgs_b_value
99 
100  integer(i4b), intent(out) :: ierr
101  real(dp), dimension(nmax), intent(inout) :: lgs_x_value
102 
103  integer(i4b) :: iter
104  integer(i4b) :: iter_max
105  integer(i4b) :: nr, k
106 #ifndef ALLOW_OPENAD
107  real(dp), allocatable, dimension(:) :: lgs_x_value_prev
108 #else
109  real(dp), dimension(nmax) :: lgs_x_value_prev
110  real(dp) :: temp1, temp2
111  logical :: isnanflag1, isnanflag2, isnanflag3
112 #endif
113  real(dp) :: b_nr
114  logical :: flag_convergence
115 
116 #if (ITER_MAX_SOR > 0)
117  iter_max = iter_max_sor
118 #else
119  iter_max = 1000 ! default value
120 #endif
121 
122 #ifndef ALLOW_OPENAD
123  allocate(lgs_x_value_prev(nmax))
124 #endif
125 
126  iter_loop : do iter=1, iter_max
127 
128  lgs_x_value_prev = lgs_x_value
129 
130  do nr=1, nmax
131 
132  b_nr = 0.0_dp
133 
134  do k=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
135  b_nr = b_nr + lgs_a_value(k)*lgs_x_value(lgs_a_index(k))
136  end do
137 
138  lgs_x_value(nr) = lgs_x_value(nr) &
139  -omega*(b_nr-lgs_b_value(nr)) &
140  /lgs_a_value(lgs_a_diag_index(nr))
141 
142  end do
143 
144  flag_convergence = .true.
145  do nr=1, nmax
146  if (abs(lgs_x_value(nr)-lgs_x_value_prev(nr)) > eps_sor) then
147  flag_convergence = .false.
148  exit
149  end if
150  end do
151 
152  if (flag_convergence) then
153  write(6,'(10x,a,i0)') 'sor_sprs: iter = ', iter
154  ierr = 0 ! convergence criterion fulfilled
155 #ifndef ALLOW_OPENAD
156  deallocate(lgs_x_value_prev)
157 #endif
158  return
159  end if
160 
161  end do iter_loop
162 
163  write(6,'(10x,a,i0)') 'sor_sprs: iter = ', iter
164  ierr = -1 ! convergence criterion not fulfilled
165 #ifndef ALLOW_OPENAD
166  deallocate(lgs_x_value_prev)
167 #endif
168 
169 #ifndef ALLOW_OPENAD
170  end subroutine sor_sprs
171 #else
172  end subroutine sor_sprs_local
173 #endif
174 
175 !-------------------------------------------------------------------------------
176 !> Solution of a system of linear equations Ax=b with tridiagonal matrix A.
177 !! @param[in] a0 a0(j) is element A_(j,j-1) of Matrix A
178 !! @param[in] a1 a1(j) is element A_(j,j) of Matrix A
179 !! @param[in] a2 a2(j) is element A_(j,j+1) of Matrix A
180 !! @param[in] b inhomogeneity vector
181 !! @param[in] nrows size of matrix A (indices run from 0 (!!!) to nrows)
182 !! @param[out] x Solution vector.
183 !<------------------------------------------------------------------------------
184 #ifndef ALLOW_OPENAD
185  subroutine tri_sle(a0, a1, a2, x, b, nrows)
186 #else
187  subroutine tri_sle_local(a0, a1, a2, x, b, nrows)
188 #endif
189 
190  implicit none
191 
192  integer(i4b), intent(in) :: nrows
193  real(dp), dimension(0:*), intent(in) :: a0, a2
194  real(dp), dimension(0:*), intent(inout) :: a1, b
195 
196  real(dp), dimension(0:*), intent(out) :: x
197 
198  real(dp), allocatable, dimension(:) :: help_x
199  integer(i4b) :: n
200 
201 !-------- Generate an upper triangular matrix
202 ! ('obere Dreiecksmatrix') --------
203 
204  do n=1, nrows
205  a1(n) = a1(n) - a0(n)/a1(n-1)*a2(n-1)
206  end do
207 
208  do n=1, nrows
209  b(n) = b(n) - a0(n)/a1(n-1)*b(n-1)
210  ! a0(n) = 0.0_dp , not needed in the following, therefore
211  ! not set
212  end do
213 
214 !-------- Iterative solution of the new system --------
215 
216  ! x(nrows) = b(nrows)/a1(nrows)
217 
218  ! do n=nrows-1, 0, -1
219  ! x(n) = (b(n)-a2(n)*x(n+1))/a1(n)
220  ! end do
221 
222  allocate(help_x(0:nrows))
223 
224  help_x(0) = b(nrows)/a1(nrows)
225 
226  do n=1, nrows
227  help_x(n) = b(nrows-n)/a1(nrows-n) &
228  -a2(nrows-n)/a1(nrows-n)*help_x(n-1)
229  end do
230 
231  do n=0, nrows
232  x(n) = help_x(nrows-n)
233  end do
234 
235  ! (The trick with the help_x was introduced in order to avoid
236  ! the negative step in the original, blanked-out loop.)
237 
238  deallocate(help_x)
239 
240  ! WARNING: Subroutine does not check for elements of the main
241  ! diagonal becoming zero. In this case it crashes even
242  ! though the system may be solvable. Otherwise ok.
243 
244 #ifndef ALLOW_OPENAD
245  end subroutine tri_sle
246 #else
247  end subroutine tri_sle_local
248 #endif
249 
250 !-------------------------------------------------------------------------------
251 !> Bilinear interpolation.
252 !<------------------------------------------------------------------------------
253 #ifndef ALLOW_OPENAD
254  function bilinint(x1, x2, y1, y2, z11, z12, z21, z22, x, y)
256  implicit none
257 
258  real(dp), intent(in) :: x1, x2, y1, y2, z11, z12, z21, z22, x, y
259 
260  real(dp) :: t, u
261  real(dp) :: bilinint
262 
263  real(dp), parameter :: I = 1.0_dp
264 
265  t = (x-x1)/(x2-x1)
266  u = (y-y1)/(y2-y1)
267 
268  bilinint = (i-t)*(i-u)*z11 + (i-t)*u*z12 + t*(i-u)*z21 + t*u*z22
269 
270  end function bilinint
271 #else
272  subroutine bilinint(x1, x2, y1, y2, z11, z12, z21, z22, x, y, retval)
273 
274  implicit none
275 
276  real(dp), intent(in) :: x1, x2, y1, y2, z11, z12, z21, z22, x, y
277 
278  real(dp) :: t, u
279  real(dp), intent(out) :: retval
280 
281  real(dp), parameter :: I = 1.0_dp
282 
283  t = (x-x1)/(x2-x1)
284  u = (y-y1)/(y2-y1)
285 
286  retval = (i-t)*(i-u)*z11 + (i-t)*u*z12 + t*(i-u)*z21 + t*u*z22
287 
288  end subroutine bilinint
289 #endif
290 
291 !-------------------------------------------------------------------------------
292 !> Computation of the complementary error function erfc(x) = 1-erf(x)
293 !! with a fractional error everywhere less than 1.2 x 10^(-7)
294 !! (formula by Press et al., 'Numerical Recipes in Fortran 77').
295 !<------------------------------------------------------------------------------
296 #if (defined(ALLOW_OPENAD) || defined(ALLOW_GRDCHK))
297  subroutine my_erfc(x, retval)
298 
299  implicit none
300 
301  real(dp), intent(in) :: x
302  real(dp), intent(out) :: retval
303 
304  real(dp) :: t, z
305 
306  z = abs(x)
307  t = 1.0_dp/(1.0_dp+0.5_dp*z)
308 
309  retval = t * exp( -z*z -1.26551223_dp &
310  + t * ( 1.00002368_dp &
311  + t * ( 0.37409196_dp &
312  + t * ( 0.09678418_dp &
313  + t * ( -0.18628806_dp &
314  + t * ( 0.27886807_dp &
315  + t * ( -1.13520398_dp &
316  + t * ( 1.48851587_dp &
317  + t * ( -0.82215223_dp &
318  + t * 0.17087277_dp ) ) ) ) ) ) ) ) )
319 
320  if (x < 0.0_dp) retval = 2.0_dp-retval
321 
322  end subroutine my_erfc
323 #endif
324 
325 #ifdef BUILD_LIS
326 #ifdef ALLOW_OPENAD
327 #include "lisf.h"
328 subroutine sico_lis_solver_local(nmax, nnz, &
329 #else
330 subroutine sico_lis_solver(nmax, nnz, &
331 #endif
332  lgs_a_ptr, lgs_a_index, &
333  lgs_a_value, lgs_b_value, lgs_x_value)
334 
335 use sico_types_m
336 
337 implicit none
338 
339 integer(i4b) :: ierr
340 integer(i4b) :: iter
341 integer(i4b) :: nc, nr
342 integer(i4b), intent(in) :: nmax
343 integer(i4b), intent(in) :: nnz
344 integer(i4b), dimension(nmax+1), intent(in) :: lgs_a_ptr
345 integer(i4b), dimension(nnz), intent(in) :: lgs_a_index
346 
347 lis_matrix :: lgs_a
348 lis_vector :: lgs_b
349 lis_vector :: lgs_x
350 lis_solver :: solver
351 
352 real(dp), dimension(nnz), intent(in) :: lgs_a_value
353 real(dp), dimension(nmax), intent(in) :: lgs_b_value
354 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD))
355 real(dp), dimension(nmax), intent(in) :: lgs_x_value
356 #else
357 real(dp), dimension(nmax), intent(inout) :: lgs_x_value
358 #endif
359 
360 character(len=256) :: ch_solver_set_option
361 
362 !LIS_INTEGER :: ierr
363 !LIS_INTEGER :: iter
364 !LIS_INTEGER :: nc, nr
365 !LIS_INTEGER, parameter :: nmax = (IMAX+1)*(JMAX+1)
366 !LIS_INTEGER, parameter :: n_sprs = 10*(IMAX+1)*(JMAX+1)
367 !LIS_INTEGER, allocatable, dimension(:) :: lgs_a_ptr, lgs_a_index
368 !LIS_INTEGER, allocatable, dimension(:) :: lgs_a_diag_index
369 !LIS_MATRIX :: lgs_a
370 !LIS_VECTOR :: lgs_b, lgs_x
371 !LIS_SCALAR, allocatable, dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
372 !LIS_SOLVER :: solver
373 !character(len=256) :: ch_solver_set_option
374 
375 ! ------ Settings for Lis
376 
377 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
378 call lis_vector_create(lis_comm_world, lgs_b, ierr)
379 call lis_vector_create(lis_comm_world, lgs_x, ierr)
380 
381 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
382 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
383 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
384 
385 do nr=1, nmax
386 
387  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
388  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
389  lgs_a_value(nc), lgs_a, ierr)
390  end do
391 
392  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
393  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
394 
395 end do
396 
397 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
398 call lis_matrix_assemble(lgs_a, ierr)
399 
400 ! ------ Solution with Lis
401 
402 call lis_solver_create(solver, ierr)
403 
404 #if (defined(LIS_OPTS))
405  ch_solver_set_option = trim(lis_opts)
406 #else
407  ch_solver_set_option = '-i bicgsafe -p jacobi '// &
408  '-maxiter 2000 -tol 1.0e-18 -initx_zeros false'
409 #endif
410 
411 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
412 call chkerr(ierr)
413 
414 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
415 call chkerr(ierr)
416 
417 call lis_solver_get_iter(solver, iter, ierr)
418 write(6,'(10x,a,i0)') 'calc_thk_sia_impl: iter = ', iter
419 
420 lgs_x_value = 0.0_dp
421 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
422 call lis_matrix_destroy(lgs_a, ierr)
423 call lis_vector_destroy(lgs_b, ierr)
424 call lis_vector_destroy(lgs_x, ierr)
425 call lis_solver_destroy(solver, ierr)
426 
427 #ifdef ALLOW_OPENAD
428 end subroutine sico_lis_solver_local
429 #else
430 end subroutine sico_lis_solver
431 #endif
432 #endif
433 !-------------------------------------------------------------------------------
434 
435 end module sico_maths_m
436 !
Several mathematical tools used by SICOPOLIS.
subroutine tri_sle(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr,
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
Declarations of kind types for SICOPOLIS.
real(dp) function bilinint(x1, x2, y1, y2, z11, z12, z21, z22, x, y)
Bilinear interpolation.