SPEED
GET_CFL.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
38! 'not' deltat will be changed if is extremely small for the time integration
45
46 subroutine get_cfl(time_step, nn_loc, loc_n_num, nm, tm, pm, sdeg, &
47 xx_loc, yy_loc, zz_loc, cs_nnz_loc, cs_loc, &
48 time_step_cfl, fmax, deltat_fixed, &
49 mpi_comm, mpi_np, mpi_id, b_failCFL)
50
51 use speed_exit_codes
52
53 implicit none
54
55 include 'SPEED.MPI'
56
57 character*3 :: deltat_fixed
58
59 integer*4 :: nm, im, nn_loc, cs_nnz_loc, mpi_comm, mpi_np, mpi_err, mpi_id
60 integer*4 :: nel_loc, ic1, ic2, n1, n2, istart, ifin
61 integer*4 :: ie, i, j, nn, mcode, smcode, sdeg_deltat_cfl, sdeg_npoints
62
63 integer*4, dimension(1) :: pos, pos2
64 integer*4, dimension(nm) :: tm, sdeg
65 integer*4, dimension(nn_loc) :: loc_n_num
66
67 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
68
69 real*8 :: time_step, time_step_cfl, fmax
70 real*8 :: length_min,length,vs_length_min,vs_length,length_vp_min,length_vp
71 real*8 :: percent_deltat
72 real*8 :: num_of_points,vs_npoints,vp_deltat_cfl
73 real*8 :: rho,lambda,mu,vs,vp
74 real*8 :: x1_length_vp_min,y1_length_vp_min,z1_length_vp_min
75 real*8 :: x2_length_vp_min,y2_length_vp_min,z2_length_vp_min
76
77 real*8, dimension(:), allocatable :: ct,ww
78 real*8, dimension(:), allocatable :: time_step_glo, vs_length_glo
79 real*8, dimension(nn_loc) :: xx_loc, yy_loc, zz_loc
80
81 real*8, dimension(:,:), allocatable :: dd
82 real*8, dimension(nm,4) :: pm
83
84 logical :: b_failCFL, b_CFL_failure = .false.
85
86 nel_loc = cs_loc(0) - 1
87
88 length_vp_min = 1.d10
89 vs_length_min = 1.d10
90
91
92
93
94 do ie = 1, nel_loc
95
96 im = cs_loc(cs_loc(ie -1) +0)
97 nn = sdeg(im) +1
98
99 smcode=sdeg(im)
100 rho = pm(im,1)
101 lambda = pm(im,2)
102 mu = pm(im,3)
103 vs = dsqrt(mu/rho)
104 vp = dsqrt((lambda+2*mu)/rho)
105
106 n1 = nn*nn*(1 -1) +nn*(1 -1) +1
107 n2 = nn*nn*(nn -1) +nn*(nn -1) +nn
108 !ic1 = cs_loc(cs_loc(ie -1) +n1)
109 !ic2 = cs_loc(cs_loc(ie -1) +n2)
110
111 istart = cs_loc(ie -1) +n1
112 ifin = cs_loc(ie -1) +n2
113
114
115! write(*,*) ic1, ic2
116! write(*,*) xx_loc(ic1),yy_loc(ic1),zz_loc(ic1)
117! write(*,*) xx_loc(ic2),yy_loc(ic2),zz_loc(ic2)
118! read(*,*)
119! write(*,*) cs_loc(ie -1) +n2, cs_loc(cs_loc(ie -1) +n2), cs_loc
120! read(*,*)
121
122 do ic1 = istart, ifin
123 do ic2 = ic1 + 1, ifin
124
125 if(cs_loc(ic1) .ne. 0 .and. cs_loc(ic2) .ne. 0) then
126
127 length=dsqrt((xx_loc(cs_loc(ic1))-xx_loc(cs_loc(ic2)))**2 + &
128 (yy_loc(cs_loc(ic1))-yy_loc(cs_loc(ic2)))**2 + &
129 (zz_loc(cs_loc(ic1))-zz_loc(cs_loc(ic2)))**2)
130
131 length_vp = length/vp
132 !write(*,*) cs_loc(ic1),cs_loc(ic2), length, vp
133
134
135 if (length_vp_min .gt. length_vp) then
136 length_vp_min = length_vp
137 x1_length_vp_min = xx_loc(cs_loc(ic1))
138 y1_length_vp_min = yy_loc(cs_loc(ic1))
139 z1_length_vp_min = zz_loc(cs_loc(ic1))
140 x2_length_vp_min = xx_loc(cs_loc(ic2))
141 y2_length_vp_min = yy_loc(cs_loc(ic2))
142 z2_length_vp_min = zz_loc(cs_loc(ic2))
143 vp_deltat_cfl = vp
144 sdeg_deltat_cfl = smcode
145 length_min = length
146 endif
147
148
149 vs_length = vs/length
150
151 if (vs_length_min .gt. vs_length) then
152 vs_length_min = vs_length
153 vs_npoints = vs
154 sdeg_npoints = smcode
155 endif
156
157
158
159 endif
160
161 enddo
162 enddo
163
164 !read(*,*)
165 enddo
166
167 nn = sdeg_deltat_cfl + 1
168
169 !write(*,*) length_vp_min,length_min,vp
170 !read(*,*)
171
172 time_step_cfl = length_vp_min !/((nn-1)**2)
173
174 allocate(time_step_glo(mpi_np),vs_length_glo(mpi_np))
175
176 call mpi_barrier(mpi_comm, mpi_err)
177 call mpi_allgather(time_step_cfl, 1, speed_double, time_step_glo, 1, speed_double, mpi_comm, mpi_err)
178 call mpi_allgather(vs_length_min, 1, speed_double, vs_length_glo, 1, speed_double, mpi_comm, mpi_err)
179
180 pos = minloc(time_step_glo)
181 pos2 = minloc(vs_length_glo)
182
183 deallocate(time_step_glo, vs_length_glo)
184
185 if((mpi_id + 1) .eq. pos(1)) then
186
187 write(*,'(A)')' '
188 write(*,'(A)') '--------Stability for time advancing algorithm --------'
189 write(*,'(A,E12.4)') 'Mininum distance beteween grid nodes :', length_vp_min*vp_deltat_cfl
190 write(*,'(A,E12.4)') 'Vp for CFL condition :', vp_deltat_cfl
191 write(*,'(A,E12.4,E12.4,E12.4)')'Min. node 1 : ',x1_length_vp_min,y1_length_vp_min,z1_length_vp_min
192 write(*,'(A,E12.4,E12.4,E12.4)')'Min. node 2 : ',x2_length_vp_min,y2_length_vp_min,z2_length_vp_min
193 write(*,'(A)') '-------------------------------------------------------'
194
195 if (time_step.le.time_step_cfl) then
196 percent_deltat=time_step/time_step_cfl*100
197 if (percent_deltat.le.1) then
198 write(*,'(A,E12.4)')'WARNING!!! deltat CFL =',&
199 time_step_cfl
200 write(*,'(A,F6.2,A)')'You are using ',percent_deltat,&
201 '% of critical time step.'
202 write(*,'(A)')'This time step is excedingly lower the deltat CFL'
203 if (deltat_fixed.eq.'not') then
204 write(*,'(A)')'deltat chosen will be substituted with 10% of deltat CFL'
205 time_step=time_step_cfl*0.1
206 write(*,'(A,E12.4)')'deltat chosen :',time_step
207 endif
208 elseif (percent_deltat.le.25) then
209 write(*,'(A,E12.4)')'OK!!! deltat CFL =',&
210 time_step_cfl
211 write(*,'(A,F6.2,A)')'You are using ',percent_deltat,&
212 '% of critical time step.'
213 else
214 write(*,'(A,E12.4)')'WARNING!!! deltat CFL =',&
215 time_step_cfl
216 write(*,'(A,F6.2,A)')'You are using ',percent_deltat,&
217 '% of critical time step.'
218 write(*,'(A)')'This could be not enough for the correct working of ABC'
219 if (deltat_fixed.eq.'not') then
220 write(*,'(A)')'deltat chosen will be substituted with 20% of deltat CFL'
221 time_step=time_step_cfl*0.2
222 write(*,'(A,E12.4)')'deltat chosen :',time_step
223 endif
224 endif
225
226 elseif (time_step .gt. time_step_cfl) then
227 write(*,'(2X,A,E12.4)')'ERROR!!! deltat CFL = ',&
228 time_step_cfl,' < deltat = ',time_step
229 write(*,'(A)')'The time advancing scheme will be unstable!'
230 if (deltat_fixed.eq.'not') then
231 write(*,'(A)')'deltat chosen will be substituted with 20% of deltat CFL'
232 time_step=time_step_cfl*0.2
233 write(*,'(A,E12.4)')'deltat chosen :',time_step
234 else
235 b_cfl_failure = .true.
236 endif
237 endif
238
239 write(*,'(A)')'-------------------------------------------------------'
240 write(*,'(A)')' '
241
242 if (b_failcfl .and. b_cfl_failure) then
243 write(*,*) 'CFL does NOT hold, asked to quit. Program finished.'
244 call mpi_finalize(mpi_err)
245 call exit(exit_cfl)
246 endif
247
248 !Writing of the results for the maximum number of points per wave length
249
250 endif
251
252 call mpi_barrier(mpi_comm, mpi_err)
253
254 allocate(ct(nn),ww(nn),dd(nn,nn))
255 call make_lgl_nw(nn,ct,ww,dd)
256
257 num_of_points = vs_length_min*(ct(1)-ct(nn))/(ct(int(nn/2))-ct(int(nn/2)+1))/fmax
258
259 !write(*,*) mpi_id, 'id', fmax, 'ppppp', vs_length_min, num_of_points
260 !call MPI_BARRIER(mpi_comm, mpi_err)
261
262 if((mpi_id + 1) .eq. pos2(1)) then
263 write(*,'(A)')'-----------Number of points per wave length------------'
264 !write(*,'(A,E12.4)')'Max. el. length :',1/vs_length_min*vs_npoints !length
265 write(*,'(A,E12.4)')'Minum vs wrt max frequency :',vs_npoints
266 write(*,'(A,E12.4)')'Points per wavelength :',num_of_points
267 write(*,'(A)') ' '
268 endif
269
270 deallocate(ct,ww,dd)
271
272 call mpi_barrier(mpi_comm, mpi_err)
273
274 return
275
276 end subroutine get_cfl
277
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29
speed_exit_codes
SPEED exit codes.
Definition MODULES.f90:25
speed_exit_codes::exit_cfl
integer, parameter exit_cfl
Definition MODULES.f90:30