SICOPOLIS V5-dev  Revision 1368
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  function bilinint(x1, x2, y1, y2, z11, z12, z21, z22, x, y)
255  implicit none
256 
257  real(dp), intent(in) :: x1, x2, y1, y2, z11, z12, z21, z22, x, y
258 
259  real(dp) :: t, u
260  real(dp) :: bilinint
261 
262  real(dp), parameter :: I = 1.0_dp
263 
264  t = (x-x1)/(x2-x1)
265  u = (y-y1)/(y2-y1)
266 
267  bilinint = (i-t)*(i-u)*z11 + (i-t)*u*z12 + t*(i-u)*z21 + t*u*z22
268 
269  end function bilinint
270 
271 !-------------------------------------------------------------------------------
272 !> Computation of the complementary error function erfc(x) = 1-erf(x)
273 !! with a fractional error everywhere less than 1.2 x 10^(-7)
274 !! (formula by Press et al., 'Numerical Recipes in Fortran 77').
275 !<------------------------------------------------------------------------------
276 #if (defined(ALLOW_OPENAD) || defined(ALLOW_GRDCHK))
277  subroutine my_erfc(x, retval)
278 
279  implicit none
280 
281  real(dp), intent(in) :: x
282  real(dp), intent(out) :: retval
283 
284  real(dp) :: t, z
285 
286  z = abs(x)
287  t = 1.0_dp/(1.0_dp+0.5_dp*z)
288 
289  retval = t * exp( -z*z -1.26551223_dp &
290  + t * ( 1.00002368_dp &
291  + t * ( 0.37409196_dp &
292  + t * ( 0.09678418_dp &
293  + t * ( -0.18628806_dp &
294  + t * ( 0.27886807_dp &
295  + t * ( -1.13520398_dp &
296  + t * ( 1.48851587_dp &
297  + t * ( -0.82215223_dp &
298  + t * 0.17087277_dp ) ) ) ) ) ) ) ) )
299 
300  if (x < 0.0_dp) retval = 2.0_dp-retval
301 
302  end subroutine my_erfc
303 #endif
304 
305 !-------------------------------------------------------------------------------
306 !> OpenAD needs a template to help differentiate through the LIS solver when
307 !! it is used. This is substituted in for adjoint modes in Antarctica with
308 !! ice shelves.
309 !<------------------------------------------------------------------------------
310 #ifdef BUILD_LIS
311 #ifdef ALLOW_OPENAD
312 #include "lisf.h"
313 subroutine sico_lis_solver_local(nmax, nnz, &
314 #else
315 subroutine sico_lis_solver(nmax, nnz, &
316 #endif
317  lgs_a_ptr, lgs_a_index, &
318  lgs_a_value, lgs_b_value, lgs_x_value)
319 
320 implicit none
321 
322 integer(i4b) :: ierr
323 integer(i4b) :: iter
324 integer(i4b) :: nc, nr
325 integer(i4b), intent(in) :: nmax
326 integer(i4b), intent(in) :: nnz
327 integer(i4b), dimension(nmax+1), intent(in) :: lgs_a_ptr
328 integer(i4b), dimension(nnz), intent(in) :: lgs_a_index
329 
330 lis_matrix :: lgs_a
331 lis_vector :: lgs_b
332 lis_vector :: lgs_x
333 lis_solver :: solver
334 
335 real(dp), dimension(nnz), intent(in) :: lgs_a_value
336 real(dp), dimension(nmax), intent(in) :: lgs_b_value
337 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD))
338 real(dp), dimension(nmax), intent(in) :: lgs_x_value
339 #else
340 real(dp), dimension(nmax), intent(inout) :: lgs_x_value
341 #endif
342 
343 character(len=256) :: ch_solver_set_option
344 
345 ! ------ Settings for Lis
346 
347 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
348 call lis_vector_create(lis_comm_world, lgs_b, ierr)
349 call lis_vector_create(lis_comm_world, lgs_x, ierr)
350 
351 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
352 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
353 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
354 
355 do nr=1, nmax
356 
357  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
358  call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
359  lgs_a_value(nc), lgs_a, ierr)
360  end do
361 
362  call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
363  call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
364 
365 end do
366 
367 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
368 call lis_matrix_assemble(lgs_a, ierr)
369 
370 ! ------ Solution with Lis
371 
372 call lis_solver_create(solver, ierr)
373 
374 #if (defined(LIS_OPTS))
375  ch_solver_set_option = trim(lis_opts)
376 #else
377  ch_solver_set_option = '-i bicgsafe -p jacobi '// &
378  '-maxiter 2000 -tol 1.0e-18 -initx_zeros false'
379 #endif
380 
381 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
382 call chkerr(ierr)
383 
384 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
385 call chkerr(ierr)
386 
387 call lis_solver_get_iter(solver, iter, ierr)
388 write(6,'(10x,a,i0)') 'calc_thk_sia_impl: iter = ', iter
389 
390 lgs_x_value = 0.0_dp
391 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
392 call lis_matrix_destroy(lgs_a, ierr)
393 call lis_vector_destroy(lgs_b, ierr)
394 call lis_vector_destroy(lgs_x, ierr)
395 call lis_solver_destroy(solver, ierr)
396 
397 #ifdef ALLOW_OPENAD
398 end subroutine sico_lis_solver_local
399 #else
400 end subroutine sico_lis_solver
401 #endif
402 #endif
403 !-------------------------------------------------------------------------------
404 
405 end module sico_maths_m
406 !
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.