SICOPOLIS V5-dev  Revision 1173
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-2017 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  public
42 
43 contains
44 
45 !-------------------------------------------------------------------------------
46 !> SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b
47 !! [matrix storage: compressed sparse row CSR,
48 !! represented by arrays lgs_a_value(values), lgs_a_index (indices)
49 !! and lgs_a_ptr (pointers)].
50 !<------------------------------------------------------------------------------
51  subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, &
52  lgs_b_value, &
53  nnz, nmax, omega, eps_sor, lgs_x_value, ierr)
54 
55  implicit none
56 
57  integer(i4b), intent(in) :: nnz, nmax
58  real(dp), intent(in) :: omega, eps_sor
59  integer(i4b), dimension(nmax+1), intent(in) :: lgs_a_ptr
60  integer(i4b), dimension(nnz), intent(in) :: lgs_a_index
61  integer(i4b), dimension(nmax), intent(in) :: lgs_a_diag_index
62  real(dp), dimension(nnz), intent(in) :: lgs_a_value
63  real(dp), dimension(nmax), intent(in) :: lgs_b_value
64 
65  integer(i4b), intent(out) :: ierr
66  real(dp), dimension(nmax), intent(inout) :: lgs_x_value
67 
68  integer(i4b) :: iter
69  integer(i4b) :: iter_max
70  integer(i4b) :: nr, k
71  real(dp), allocatable, dimension(:) :: lgs_x_value_prev
72  real(dp) :: b_nr
73  logical :: flag_convergence
74 
75 #if (ITER_MAX_SOR > 0)
76  iter_max = iter_max_sor
77 #else
78  iter_max = 1000 ! default value
79 #endif
80 
81  allocate(lgs_x_value_prev(nmax))
82 
83  iter_loop : do iter=1, iter_max
84 
85  lgs_x_value_prev = lgs_x_value
86 
87  do nr=1, nmax
88 
89  b_nr = 0.0_dp
90 
91  do k=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
92  b_nr = b_nr + lgs_a_value(k)*lgs_x_value(lgs_a_index(k))
93  end do
94 
95  lgs_x_value(nr) = lgs_x_value(nr) &
96  -omega*(b_nr-lgs_b_value(nr)) &
97  /lgs_a_value(lgs_a_diag_index(nr))
98 
99  end do
100 
101  flag_convergence = .true.
102  do nr=1, nmax
103  if (abs(lgs_x_value(nr)-lgs_x_value_prev(nr)) > eps_sor) then
104  flag_convergence = .false.
105  exit
106  end if
107  end do
108 
109  if (flag_convergence) then
110  write(6,'(10x,a,i0)') 'sor_sprs: iter = ', iter
111  ierr = 0 ! convergence criterion fulfilled
112  deallocate(lgs_x_value_prev)
113  return
114  end if
115 
116  end do iter_loop
117 
118  write(6,'(10x,a,i0)') 'sor_sprs: iter = ', iter
119  ierr = -1 ! convergence criterion not fulfilled
120  deallocate(lgs_x_value_prev)
121 
122  end subroutine sor_sprs
123 
124 !-------------------------------------------------------------------------------
125 !> Solution of a system of linear equations Ax=b with tridiagonal matrix A.
126 !! @param[in] a0 a0(j) is element A_(j,j-1) of Matrix A
127 !! @param[in] a1 a1(j) is element A_(j,j) of Matrix A
128 !! @param[in] a2 a2(j) is element A_(j,j+1) of Matrix A
129 !! @param[in] b inhomogeneity vector
130 !! @param[in] nrows size of matrix A (indices run from 0 (!!!) to nrows)
131 !! @param[out] x Solution vector.
132 !<------------------------------------------------------------------------------
133  subroutine tri_sle(a0, a1, a2, x, b, nrows)
135  implicit none
136 
137  integer(i4b), intent(in) :: nrows
138  real(dp), dimension(0:*), intent(in) :: a0, a2
139  real(dp), dimension(0:*), intent(inout) :: a1, b
140 
141  real(dp), dimension(0:*), intent(out) :: x
142 
143  real(dp), allocatable, dimension(:) :: help_x
144  integer(i4b) :: n
145 
146 !-------- Generate an upper triangular matrix
147 ! ('obere Dreiecksmatrix') --------
148 
149  do n=1, nrows
150  a1(n) = a1(n) - a0(n)/a1(n-1)*a2(n-1)
151  end do
152 
153  do n=1, nrows
154  b(n) = b(n) - a0(n)/a1(n-1)*b(n-1)
155  ! a0(n) = 0.0_dp , not needed in the following, therefore
156  ! not set
157  end do
158 
159 !-------- Iterative solution of the new system --------
160 
161  ! x(nrows) = b(nrows)/a1(nrows)
162 
163  ! do n=nrows-1, 0, -1
164  ! x(n) = (b(n)-a2(n)*x(n+1))/a1(n)
165  ! end do
166 
167  allocate(help_x(0:nrows))
168 
169  help_x(0) = b(nrows)/a1(nrows)
170 
171  do n=1, nrows
172  help_x(n) = b(nrows-n)/a1(nrows-n) &
173  -a2(nrows-n)/a1(nrows-n)*help_x(n-1)
174  end do
175 
176  do n=0, nrows
177  x(n) = help_x(nrows-n)
178  end do
179 
180  ! (The trick with the help_x was introduced in order to avoid
181  ! the negative step in the original, blanked-out loop.)
182 
183  deallocate(help_x)
184 
185  ! WARNING: Subroutine does not check for elements of the main
186  ! diagonal becoming zero. In this case it crashes even
187  ! though the system may be solvable. Otherwise ok.
188 
189  end subroutine tri_sle
190 
191 !-------------------------------------------------------------------------------
192 !> Bilinear interpolation.
193 !<------------------------------------------------------------------------------
194  function bilinint(x1, x2, y1, y2, z11, z12, z21, z22, x, y)
196  implicit none
197 
198  real(dp), intent(in) :: x1, x2, y1, y2, z11, z12, z21, z22, x, y
199 
200  real(dp) :: t, u
201  real(dp) :: bilinint
202 
203  real(dp), parameter :: I = 1.0_dp
204 
205  t = (x-x1)/(x2-x1)
206  u = (y-y1)/(y2-y1)
207 
208  bilinint = (i-t)*(i-u)*z11 + (i-t)*u*z12 + t*(i-u)*z21 + t*u*z22
209 
210  end function bilinint
211 
212 !-------------------------------------------------------------------------------
213 
214 end module sico_maths_m
215 !
subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, lgs_b_value, nnz, nmax, omega, eps_sor, lgs_x_value, ierr)
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
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.
Declarations of kind types for SICOPOLIS.
real(dp) function bilinint(x1, x2, y1, y2, z11, z12, z21, z22, x, y)
Bilinear interpolation.