solver_internal.f90 8.89 KB
Newer Older
1
! solver_interal.f90 --- internal state and logic of the solver, helper functions 
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
!
! Copyright  (C)  2010  Florian Wilhelm <Florian.Wilhelm@kit.edu>
!
! Version: 1.0
! Keywords: scales ppm solver internal
! Author: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! Maintainer: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! URL: https://www.dkrz.de/redmine/projects/show/scales-ppm
!
! Redistribution and use in source and binary forms, with or without
! modification, are  permitted provided that the following conditions are
! met:
!
! Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
!
! Neither the name of the DKRZ GmbH nor the names of its contributors
! may be used to endorse or promote products derived from this software
! without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! Commentary:
Florian Wilhelm's avatar
Florian Wilhelm committed
39
!
40
41
42
43
44
45
! Avoid including/importing this module as it is supposed to be only used
! from other solver modules.
!
! Code:
!

46
MODULE solver_internal
47
    USE mo_kind, ONLY: wp, dp
48
    USE ppm_extents, ONLY: extent_type => extent !, extent_start, extent_end, extent_size
49
50

    IMPLICIT NONE
Florian Wilhelm's avatar
Florian Wilhelm committed
51

52
    PRIVATE
Florian Wilhelm's avatar
Florian Wilhelm committed
53

54
55
56
57
58
59
60
#ifdef USE_MPI
    INCLUDE 'mpif.h'
    INTEGER :: mpi_real_dp
#else
    INTEGER, PARAMETER :: mpi_comm_world = 0
#endif

Florian Wilhelm's avatar
Florian Wilhelm committed
61
    TYPE stencil_type
62
        REAL(wp), DIMENSION(:,:), POINTER :: zonal => null(), meridional => null(), central => null()
63
64
        TYPE(extent_type), DIMENSION(:), POINTER :: extent => null()
        REAL(wp) :: shift = 0.0_wp ! For shifting when determing least dominant eigenvalue
65
    END TYPE stencil_type
Florian Wilhelm's avatar
Florian Wilhelm committed
66

67
    TYPE solver_config_type
68
69
        INTEGER :: solver              ! which solver to use
        INTEGER :: preconditioner      ! which preconditioner to use
70
71
72
73
74
        LOGICAL :: do_exchange         ! perform boundary exchanges or not
        REAL(wp):: tol                 ! tolerance for solver
        INTEGER :: maxiter             ! maximum iteration
        REAL(wp):: schwarz_tol         ! tolerance for solver
        INTEGER :: schwarz_maxiter     ! maximum iteration
75
76
        REAL(wp):: lambda_min          ! minimum absolut eigenvalue
        REAL(wp):: lambda_max          ! maximum absolut eigenvalue
77
78
        REAL(wp):: lambda_tol          ! tolerance for calculation of lambda_min, lambda_max
        INTEGER :: lambda_maxiter      ! maxiter for calculation of lambda_min, lambda_max
79
        REAL(wp):: ssor_param          ! SSOR parameter for preconditioner
80
        INTEGER :: icc_param           ! level of fill-in for ICC preconditioner
Florian Wilhelm's avatar
Florian Wilhelm committed
81
        INTEGER :: cheby_miniter       ! minimum number of iterations
82
        INTEGER :: cheby_checkrate     ! rate of checking for break condition
Florian Wilhelm's avatar
Florian Wilhelm committed
83
        INTEGER :: schwarz_miniter     ! minimum number of iterations
84
        INTEGER :: schwarz_checkrate   ! rate of checking for break condition
85
86
    END TYPE solver_config_type

Florian Wilhelm's avatar
Florian Wilhelm committed
87
    !-------------------------------------------------------------------
88
    ! Internal State Variables
Florian Wilhelm's avatar
Florian Wilhelm committed
89
    !-------------------------------------------------------------------
90
91
    TYPE(stencil_type), SAVE :: stencil
    TYPE(solver_config_type), SAVE :: config
Florian Wilhelm's avatar
Florian Wilhelm committed
92

93
94
95
96
97
98
99
    INTERFACE abort_unless_normal
        MODULE PROCEDURE abort_unless_normal0
        MODULE PROCEDURE abort_unless_normal1
        MODULE PROCEDURE abort_unless_normal2
        MODULE PROCEDURE abort_unless_normal3
    END INTERFACE

100
    PUBLIC :: stencil_type, solver_config_type, stencil, config, &
101
         int2str, abort_unless_normal, clear_halos
102
103
104
#ifdef USE_MPI
    PUBLIC :: mpi_real_dp, mpi_comm_world, mpi_sum
#endif
Florian Wilhelm's avatar
Florian Wilhelm committed
105
106
    CONTAINS

107
    ! Returns string representation of given integer i
108
109
110
111
    PURE FUNCTION int2str(i) RESULT(cout)
        INTEGER, INTENT(IN) :: i
        CHARACTER(len=32) :: cin
        CHARACTER(len=floor(log10(REAL(max(1,abs(i)))))+1-min(sign(1,i),0)) :: cout
Florian Wilhelm's avatar
Florian Wilhelm committed
112

113
114
        WRITE(cin,"(I32)") i
        cout = ADJUSTL(cin)
Florian Wilhelm's avatar
Florian Wilhelm committed
115

116
    END FUNCTION int2str
117

118
119
    

120
121
    ! Aborts if abnormal (NaN, Inf etc.) entries are present
    SUBROUTINE abort_unless_normal0(x, x_str)
122
        USE ppm_base, ONLY: abort_ppm
123
124
125
126
        USE ieee_arithmetic, ONLY: is_normal => ieee_is_normal

        REAL(wp), INTENT(IN) :: x
        CHARACTER(len=*), INTENT(IN) :: x_str ! string representation of x
Florian Wilhelm's avatar
Florian Wilhelm committed
127
        INTEGER :: error, rank
128

Florian Wilhelm's avatar
Florian Wilhelm committed
129
        CALL MPI_COMM_RANK(mpi_comm_world, rank, error)
130
131

        IF ( .NOT. is_normal(x) ) THEN
Florian Wilhelm's avatar
Florian Wilhelm committed
132
133
134
135
            CALL abort_ppm("Not normal number found in " // x_str // &
                " on rank " // int2str(rank) // ".", &
                __FILE__, &
                __LINE__)
136
137
138
139
140
141
        ENDIF

    END SUBROUTINE abort_unless_normal0

    ! Aborts if abnormal (NaN, Inf etc.) entries are present
    SUBROUTINE abort_unless_normal1(x, x_str)
142
        USE ppm_base, ONLY: abort_ppm
143
144
145
146
        USE ieee_arithmetic, ONLY: is_normal => ieee_is_normal

        REAL(wp), INTENT(IN) :: x(:)
        CHARACTER(len=*), INTENT(IN) :: x_str ! string representation of x
Florian Wilhelm's avatar
Florian Wilhelm committed
147
        INTEGER :: i, ie, error, rank
148
149
150

        ie = SIZE(x,1)

Florian Wilhelm's avatar
Florian Wilhelm committed
151
        CALL MPI_COMM_RANK(mpi_comm_world, rank, error)
152
153
154

        DO i=1,ie
            IF ( .NOT. is_normal(x(i)) ) THEN
Florian Wilhelm's avatar
Florian Wilhelm committed
155
156
157
158
159
                CALL abort_ppm("Not normal number found in " // x_str // &
                    ", i=" // int2str(i) // &
                    " on rank " // int2str(rank) // ".", &
                    __FILE__, &
                    __LINE__)
160
161
162
163
164
165
166
            ENDIF
        ENDDO

    END SUBROUTINE abort_unless_normal1

    ! Aborts if abnormal (NaN, Inf etc.) entries are present
    SUBROUTINE abort_unless_normal2(x, x_str)
167
        USE ppm_base, ONLY: abort_ppm
168
169
170
171
        USE ieee_arithmetic, ONLY: is_normal => ieee_is_normal

        REAL(wp), INTENT(IN) :: x(:,:)
        CHARACTER(len=*), INTENT(IN) :: x_str ! string representation of x
Florian Wilhelm's avatar
Florian Wilhelm committed
172
        INTEGER :: i, j, ie, je, error, rank
173
174
175
176

        ie = SIZE(x,1)
        je = SIZE(x,2)

Florian Wilhelm's avatar
Florian Wilhelm committed
177
        CALL MPI_COMM_RANK(mpi_comm_world, rank, error)
178
179
180
181

        DO j=1,je
            DO i=1,ie
                IF ( .NOT. is_normal(x(i,j)) ) THEN
Florian Wilhelm's avatar
Florian Wilhelm committed
182
183
184
185
186
                    CALL abort_ppm("Not normal number found in " // x_str // &
                        ", i=" // int2str(i) // ", j=" // int2str(j) // &
                        " on rank " // int2str(rank) // ".", &
                        __FILE__, &
                        __LINE__)
187
188
189
190
191
192
193
194
                ENDIF
            ENDDO
        ENDDO

    END SUBROUTINE abort_unless_normal2

    ! Aborts if abnormal (NaN, Inf etc.) entries are present
    SUBROUTINE abort_unless_normal3(x, x_str)
195
        USE ppm_base, ONLY: abort_ppm
196
197
198
199
        USE ieee_arithmetic, ONLY: is_normal => ieee_is_normal

        REAL(wp), INTENT(IN) :: x(:,:,:)
        CHARACTER(len=*), INTENT(IN) :: x_str ! string representation of x
Florian Wilhelm's avatar
Florian Wilhelm committed
200
        INTEGER :: i, j, k, ie, je, ke, error, rank
201
202
203
204
205

        ie = SIZE(x,1)
        je = SIZE(x,2)
        ke = SIZE(x,3)

Florian Wilhelm's avatar
Florian Wilhelm committed
206
        CALL MPI_COMM_RANK(mpi_comm_world, rank, error)
207
208
209
210
211

        DO k=1,ke
            DO j=1,je
                DO i=1,ie
                    IF ( .NOT. is_normal(x(i,j,k)) ) THEN
Florian Wilhelm's avatar
Florian Wilhelm committed
212
213
214
215
216
217
                        CALL abort_ppm("Not normal number found in " // x_str &
                            // ", i=" // int2str(i) // ", j=" // int2str(j) // &
                            ", k=" // int2str(k) // " on rank " &
                            // int2str(rank) // ".", &
                            __FILE__, &
                            __LINE__)
218
219
220
221
222
223
224
                    ENDIF
                ENDDO
            ENDDO
        ENDDO

    END SUBROUTINE abort_unless_normal3

225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    ! sets halos to zero
    SUBROUTINE clear_halos(x, ext_x)
        USE partition_descriptors, ONLY: extent_start, extent_end

        REAL(wp), INTENT(INOUT) :: x(:,:)
        TYPE(extent_type), INTENT(IN) :: ext_x(:)  
        INTEGER :: jb, ib, il, jl, ie, je

        ie = SIZE(x,1)
        je = SIZE(x,2)
        ib = extent_start(ext_x(1))
        jb = extent_start(ext_x(2))
        il = extent_end(ext_x(1))
        jl = extent_end(ext_x(2))

        x(1:ib-1,:) = 0._wp
        x(il+1:ie,:) = 0._wp
        x(ib:il,1:jb-1) = 0._wp
        x(ib:il,jl+1:je) = 0._wp

    END SUBROUTINE clear_halos

247
END MODULE solver_internal