FE-Project
Loading...
Searching...
No Matches
scale_timeint_rk_butcher_tab.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
27#include "scaleFElib.h"
28module scale_timeint_rk_butcher_tab
29 !-----------------------------------------------------------------------------
30 !
31 !++ used modules
32 !
33 use scale_precision
34 use scale_io
35 use scale_prc
36 !-----------------------------------------------------------------------------
37 implicit none
38 private
39
40 !-----------------------------------------------------------------------------
41 !
42 !++ Public type & procedure
43 !
44 public :: timeint_rk_butcher_tab_get_info
45 public :: timeint_rk_butcher_tab_get
46
47 !-----------------------------------------------------------------------------
48 !
49 !++ Private procedure
50 !
51
52 !-----------------------------------------------------------------------------
53 !
54 !++ Private parameters & variables
55 !
56
57 !-----------------------------------------------------------------------------
58contains
59!OCL SERIAL
60 subroutine timeint_rk_butcher_tab_get_info( rk_scheme_name, &
61 nstage, tend_buf_size, low_storage_flag, imex_flag )
62 implicit none
63 character(len=*), intent(in) :: rk_scheme_name
64 integer, intent(out) :: nstage
65 integer, intent(out) :: tend_buf_size
66 logical, intent(out) :: low_storage_flag
67 logical, intent(out) :: imex_flag
68 !-----------------------------------------------------------
69
70 low_storage_flag = .false.
71
72 select case(rk_scheme_name)
73 case( 'ERK_1s1o', 'ERK_Euler' )
74 nstage = 1
75 tend_buf_size = 1
76 imex_flag = .false.
77 case( 'ERK_4s4o', 'ERK_RK4' )
78 nstage = 4
79 tend_buf_size = 1
80 imex_flag = .false.
81 case( 'ERK_SSP_2s2o' )
82 nstage = 2
83 tend_buf_size = 1
84 low_storage_flag = .true.
85 imex_flag = .false.
86 case( 'ERK_SSP_3s3o' ) ! SSP coefficient: 1, effective SSP coefficient Ceff: 1/3
87 nstage = 3
88 tend_buf_size = 1
89 low_storage_flag = .true.
90 imex_flag = .false.
91 case( 'ERK_SSP_4s3o' )
92 nstage = 4
93 tend_buf_size = 1
94 low_storage_flag = .true.
95 imex_flag = .false.
96 case( 'ERK_SSP_5s3o_2N2*' ) ! Higueras and Roldan, 2018: New third order low-storage SSP explicit Runge–Kutta methods
97 nstage = 5
98 tend_buf_size = 1
99 low_storage_flag = .true.
100 imex_flag = .false.
101 case ( 'ERK_SSP_10s4o_2N' ) ! Ketcheson, 2008: HIGHLY EFFICIENT STRONG STABILITY-PRESERVING RUNGE–KUTTA METHODS WITH LOW-STORAGE IMPLEMENTATIONS, SIAM J. SCI. COMPUT.
102 ! SSP coefficient: 6, effective SSP coefficient Ceff: 0.6
103 nstage = 10
104 tend_buf_size = 1
105 low_storage_flag = .true.
106 imex_flag = .false.
107 case( 'IMEX_ARK232' ) ! Giraldo et al. (2013): Implicit-Explicit Formulations of a Three-Dimensional Nonhydrostatic Unified Model
108 ! of the Atmosphere (NUMA), SIAM J. Sci. Comp.
109 nstage = 3
110 tend_buf_size = 3
111 imex_flag = .true.
112 case( 'IMEX_ARK324' ) ! Kennedy and Carpenter, 2003: Additive Runge-Kutta schemes for convection-diffusion-reaction equations, Appl. Numer. Math.
113 nstage = 4
114 tend_buf_size = 4
115 imex_flag = .true.
116 case default
117 log_error("timeint_rk_butcher_tab_get_info",*) trim(rk_scheme_name)//' is not supported. Check!'
118 call prc_abort
119 end select
120
121 return
122 end subroutine timeint_rk_butcher_tab_get_info
123
124!OCL SERIAL
125 subroutine timeint_rk_butcher_tab_get( &
126 rk_scheme_name, nstage, imex_flag, & ! (in)
127 coef_a_ex, coef_b_ex, coef_c_ex, & ! (in)
128 coef_sig_ex, coef_gam_ex, & ! (in)
129 coef_a_im, coef_b_im, coef_c_im, & ! (in)
130 tend_buf_indmap ) ! (out)
131
132 implicit none
133 character(len=*), intent(in) :: rk_scheme_name
134 integer, intent(in) :: nstage
135 logical, intent(in) :: imex_flag
136 real(rp), intent(out) :: coef_a_ex(nstage,nstage)
137 real(rp), intent(out) :: coef_b_ex(nstage)
138 real(rp), intent(out) :: coef_c_ex(nstage)
139 real(rp), intent(out) :: coef_sig_ex(nstage+1,nstage)
140 real(rp), intent(out) :: coef_gam_ex(nstage+1,nstage)
141 real(rp), intent(out) :: coef_a_im(nstage,nstage)
142 real(rp), intent(out) :: coef_b_im(nstage)
143 real(rp), intent(out) :: coef_c_im(nstage)
144 integer, intent(out) :: tend_buf_indmap(nstage)
145
146 integer :: n
147 real(rp) :: alp, gam, del
148 logical :: call_shuosher2butcher
149
150 !-----------------------------------------------------------
151
152 coef_a_ex(:,:) = 0.0_rp
153 coef_b_ex(:) = 0.0_rp
154 coef_c_ex(:) = 0.0_rp
155
156 coef_sig_ex(:,:) = 0.0_rp
157 coef_gam_ex(:,:) = 0.0_rp
158
159 if (imex_flag) then
160 coef_a_im(:,:) = 0.0_rp
161 coef_b_im(:) = 0.0_rp
162 coef_c_im(:) = 0.0_rp
163 end if
164
165 call_shuosher2butcher = .false.
166
167 !---
168
169 select case(rk_scheme_name)
170 case( 'ERK_1s1o', 'ERK_Euler' )
171
172 coef_a_ex(1,1) = 0.0_rp
173 coef_b_ex(:) = (/ 1.0_rp /)
174
175 coef_sig_ex(2,1 ) = 1.0_rp
176 coef_gam_ex(2,1 ) = 1.0_rp
177
178 tend_buf_indmap(:) = 1
179
180 case( 'ERK_4s4o', 'ERK_RK4' )
181
182 coef_a_ex(2,1) = 0.5_rp
183 coef_a_ex(3,2) = 0.5_rp
184 coef_a_ex(4,3) = 1.0_rp
185 coef_b_ex(:) = (/ 1.0_rp, 2.0_rp, 2.0_rp, 1.0_rp /)/6.0_rp
186
187 tend_buf_indmap(:) = 1
188
189 case( 'ERK_SSP_2s2o' )
190
191 coef_sig_ex(2,1 ) = 1.0_rp
192 coef_sig_ex(3,1:2) = (/ 1.0_rp, 1.0_rp /) / 2.0_rp
193 coef_gam_ex(2,1 ) = 1.0_rp
194 coef_gam_ex(3,2 ) = 1.0_rp / 2.0_rp
195
196 tend_buf_indmap(:) = 1
197 call_shuosher2butcher = .true.
198
199 case( 'ERK_SSP_3s3o' )
200
201 coef_sig_ex(2,1 ) = 1.0_rp
202 coef_sig_ex(3,1:2) = (/ 3.0_rp, 1.0_rp /) / 4.0_rp
203 coef_sig_ex(4,1:3) = (/ 1.0_rp, 0.0_rp, 2.0_rp /) / 3.0_rp
204 coef_gam_ex(2,1 ) = 1.0_rp
205 coef_gam_ex(3,2 ) = 1.0_rp / 4.0_rp
206 coef_gam_ex(4,3 ) = 2.0_rp / 3.0_rp
207
208 tend_buf_indmap(:) = 1
209 call_shuosher2butcher = .true.
210
211 case( 'ERK_SSP_4s3o' )
212
213 coef_sig_ex(2,1 ) = 1.0_rp
214 coef_sig_ex(3,1:2) = (/ 0.0_rp, 1.0_rp /)
215 coef_sig_ex(4,1:3) = (/ 2.0_rp, 0.0_rp, 1.0_rp /) / 3.0_rp
216 coef_sig_ex(5,1:4) = (/ 0.0_rp, 0.0_rp, 0.0_rp, 1.0_rp /)
217 coef_gam_ex(2,1 ) = 0.5_rp
218 coef_gam_ex(3,2 ) = 0.5_rp
219 coef_gam_ex(4,3 ) = 1.0_rp / 6.0_rp
220 coef_gam_ex(5,4 ) = 0.5_rp
221
222 tend_buf_indmap(:) = 1
223 call_shuosher2butcher = .true.
224
225 case( 'ERK_SSP_5s3o_2N2*' )
226
227 coef_sig_ex(2,1 ) = 1.0_rp
228 coef_sig_ex(3,1:2) = (/ 0.0_rp, 1.0_rp /)
229 coef_sig_ex(4,1:3) = (/ 0.682342861037239_rp, 0.0_rp, 0.317657138962761_rp /)
230 coef_sig_ex(5,1:4) = (/ 0.0_rp, 0.0_rp, 0.0_rp, 1.0_rp /)
231 coef_sig_ex(6,1:5) = (/ 0.045230974482400_rp, 0.0_rp, 0.0_rp, 0.0_rp, 0.954769025517600_rp /)
232 coef_gam_ex(2,1 ) = 0.465388589249323_rp
233 coef_gam_ex(3,2 ) = 0.465388589249323_rp
234 coef_gam_ex(4,3 ) = 0.124745797313998_rp
235 coef_gam_ex(5,4 ) = 0.465388589249323_rp
236 coef_gam_ex(6,5 ) = 0.154263303748666_rp
237
238 tend_buf_indmap(:) = 1
239 call_shuosher2butcher = .true.
240
241 case( 'ERK_SSP_10s4o_2N' )
242
243 do n=1, 4
244 coef_sig_ex(n+1,n) = 1.0_rp
245 coef_gam_ex(n+1,n) = 1.0_rp / 6.0_rp
246 end do
247 coef_sig_ex(6,1:5) = (/ 3.0_rp, 0.0_rp, 0.0_rp, 0.0_rp, 2.0_rp /) / 5.0_rp
248 coef_gam_ex(6,5 ) = 1.0_rp / 15.0_rp
249 do n=6, 9
250 coef_sig_ex(n+1,n) = 1.0_rp
251 coef_gam_ex(n+1,n) = 1.0_rp / 6.0_rp
252 end do
253 coef_sig_ex(11,1:10) = (/ 0.2_rp, 0.0_rp, 0.0_rp, 0.0_rp, 1.8_rp, 0.0_rp, 0.0_rp, 0.0_rp, 0.0_rp, 3.0_rp /) * 0.2_rp
254 coef_gam_ex(11,1:10) = (/ 0.0_rp, 0.0_rp, 0.0_rp, 0.0_rp, 1.8_rp, 0.0_rp, 0.0_rp, 0.0_rp, 0.0_rp, 3.0_rp /) / 30.0_rp
255
256 tend_buf_indmap(:) = 1
257 call_shuosher2butcher = .true.
258
259 case('IMEX_ARK232')
260
261 alp = (3.0_rp + 2.0_rp*sqrt(2.0_rp))/6.0_rp
262 gam = 1.0_rp - 1.0_rp/sqrt(2.0_rp)
263 del = 1.0_rp/(2.0_rp*sqrt(2.0_rp))
264
265 coef_a_ex(2,1) = 2.0_rp*gam
266 coef_a_ex(3,:) = (/ 1.0_rp - alp, alp, 0.0_rp /)
267 coef_b_ex(:) = (/ del, del, gam /)
268
269 coef_a_im(2,:) = (/ gam, gam, 0.0_rp /)
270 coef_a_im(3,:) = (/ del, del, gam /)
271 coef_b_im(:) = coef_b_ex(:)
272
273 tend_buf_indmap(:) = (/ 1, 2, 3 /)
274
275 case('IMEX_ARK324')
276
277 coef_a_ex(2,1) = 1767732205903.0_rp/2027836641118.0_rp
278 coef_a_ex(3,1) = 5535828885825.0_rp/10492691773637.0_rp
279 coef_a_ex(3,2) = 788022342437.0_rp/10882634858940.0_rp
280 coef_a_ex(4,1) = 6485989280629.0_rp/16251701735622.0_rp
281 coef_a_ex(4,2) = -4246266847089.0_rp/9704473918619.0_rp
282 coef_a_ex(4,3) = 10755448449292.0_rp/10357097424841.0_rp
283 coef_b_ex(1) = 1471266399579.0_rp/7840856788654.0_rp
284 coef_b_ex(2) = -4482444167858.0_rp/7529755066697.0_rp
285 coef_b_ex(3) = 11266239266428.0_rp/11593286722821.0_rp
286 coef_b_ex(4) = 1767732205903.0_rp/4055673282236.0_rp
287
288 coef_a_im(2,1:2) = 1767732205903.0_rp/4055673282236.0_rp
289 coef_a_im(3,1) = 2746238789719.0_rp/10658868560708.0_rp
290 coef_a_im(3,2) = -640167445237.0_rp/6845629431997.0_rp
291 coef_a_im(3,3) = 1767732205903.0_rp/4055673282236.0_rp
292 coef_a_im(4,:) = coef_b_ex(:)
293 coef_b_im(:) = coef_b_ex(:)
294
295 tend_buf_indmap(:) = (/ 1, 2, 3, 4 /)
296
297 case default
298 log_error("timeint_rk_butcher_tab_get",*) trim(rk_scheme_name)//' is not supported. Check!'
299 call prc_abort
300
301 end select
302
303 !--
304 if ( call_shuosher2butcher ) then
305 call shuosher2butcher( &
306 coef_sig_ex, coef_gam_ex, nstage, & ! (in)
307 coef_a_ex, coef_b_ex ) ! (out)
308 end if
309
310 !--
311 coef_c_ex(:) = 0.0_rp
312 do n=1, nstage
313 coef_c_ex(n) = sum( coef_a_ex(n,1:n) )
314 end do
315 if ( imex_flag ) then
316 coef_c_ex(:) = 0.0_rp
317 do n=1, nstage
318 coef_c_ex(n) = sum( coef_a_ex(n,1:n) )
319 end do
320 end if
321
322 return
323 end subroutine timeint_rk_butcher_tab_get
324
325 !- private --------------------------------------------
326
329!OCL SERIAL
330 subroutine shuosher2butcher( &
331 coef_sig, coef_gam, nstage, &
332 coef_a, coef_b )
333
334 use scale_linalgebra, only: &
335 linalgebra_solvelineq
336 implicit none
337
338 integer, intent(in) :: nstage
339 real(rp), intent(in) :: coef_sig(nstage+1,nstage)
340 real(rp), intent(in) :: coef_gam(nstage+1,nstage)
341 real(rp), intent(out) :: coef_a(nstage,nstage)
342 real(rp), intent(out) :: coef_b(nstage)
343
344 integer :: n
345 real(rp) :: l(nstage,nstage)
346 !-----------------------------------------------------------
347
348 ! L=I-SIG
349 l(:,:) = - coef_sig(1:nstage,:)
350 do n=1, nstage
351 l(n,n) = 1.0_rp + l(n,n)
352 end do
353 ! A=L^-1 GAM
354 call linalgebra_solvelineq( l(:,:), coef_gam(1:nstage,:), coef_a(:,:) )
355
356 ! b = GAM_N,: + SIG_N,: * A
357 do n=1, nstage
358 coef_b(n) = coef_gam(nstage+1,n) + sum( coef_sig(nstage+1,:) * coef_a(:,n) )
359 end do
360
361 return
362 end subroutine shuosher2butcher
363
364end module scale_timeint_rk_butcher_tab
scale_linalgebra::linalgebra_solvelineq
Definition scale_linalgebra.F90:32
scale_linalgebra
module common / Linear algebra
Definition scale_linalgebra.F90:11
scale_timeint_rk_butcher_tab
module common / Runge-Kutta scheme
Definition scale_timeint_rk_butcher_tab.F90:28
scale_timeint_rk_butcher_tab::timeint_rk_butcher_tab_get_info
subroutine, public timeint_rk_butcher_tab_get_info(rk_scheme_name, nstage, tend_buf_size, low_storage_flag, imex_flag)
Definition scale_timeint_rk_butcher_tab.F90:62
scale_timeint_rk_butcher_tab::timeint_rk_butcher_tab_get
subroutine, public timeint_rk_butcher_tab_get(rk_scheme_name, nstage, imex_flag, coef_a_ex, coef_b_ex, coef_c_ex, coef_sig_ex, coef_gam_ex, coef_a_im, coef_b_im, coef_c_im, tend_buf_indmap)
Definition scale_timeint_rk_butcher_tab.F90:131
Generated by doxygen 1.10.0