Choral
ode_SL_ms_mod.f90
Go to the documentation of this file.
1 
9 
10 
12 
13  !$ use OMP_LIB
16  use real_type
17  use algebra_set
18  use basic_tools
19  use r1d_mod, only: xpay, copy, axpby
20  use ode_def
21  use ode_problem_mod
23  use krylov_mod
24 
25  implicit none
26  private
27 
28  public :: solve_ode_sl_ms
29  public :: create_ode_sl_ms_sol
30  public :: set_solver_ode_sl_ms
31  public :: memsize_ode_sl_ms
32  public :: check_ode_method_sl_ms
33 
34 
35  real(RP), parameter :: c4_3 = 4._rp / 3._rp
36  real(RP), parameter :: c2_3 = 2._rp / 3._rp
37  real(RP), parameter :: c_3 = 1._rp / 3._rp
38  real(RP), parameter :: c2_11 = 2._rp / 11._rp
39  real(RP), parameter :: c6_11 = 6._rp / 11._rp
40  real(RP), parameter :: c9_11 = 9._rp / 11._rp
41  real(RP), parameter :: c18_11 = 18._rp / 11._rp
42  real(RP), parameter :: c3_25 = 3._rp / 25._rp
43  real(RP), parameter :: c12_25 = 12._rp / 25._rp
44  real(RP), parameter :: c16_25 = 16._rp / 25._rp
45  real(RP), parameter :: c36_25 = 36._rp / 25._rp
46  real(RP), parameter :: c48_25 = 48._rp / 25._rp
47  real(RP), parameter :: c72_25 = 72._rp / 25._rp
48  real(RP), parameter :: c12_137 = 12._rp /137._rp
49  real(RP), parameter :: c60_137 = 60._rp /137._rp
50  real(RP), parameter :: c75_137 = 75._rp /137._rp
51  real(RP), parameter :: c200_137= 200._rp /137._rp
52  real(RP), parameter :: c300_137= 300._rp /137._rp
53  real(RP), parameter :: c600_137= 600._rp /137._rp
54 
55 contains
56 
58  function check_ode_method_sl_ms(method) result(b)
59  logical :: b
60  integer, intent(in) :: method
61 
62  b = .false.
63 
64  select case(method)
67 
68  b = .true.
69  end select
70 
71  end function check_ode_method_sl_ms
72 
73 
81  subroutine memsize_ode_sl_ms(n_V, n_FY, method)
82  integer, intent(out) :: n_V, n_FY
83  integer, intent(in) :: method
84 
85  select case(method)
86 
87  case(ode_fbe, ode_fe)
88  n_v = 1 ; n_fy = 1
89 
90  case(ode_cnab2)
91  n_v = 1 ; n_fy = 2
92 
94  n_v = 2 ; n_fy = 2
95 
96  case(ode_bdfsbdf3)
97  n_v = 3 ; n_fy = 3
98 
99  case(ode_bdfsbdf4)
100  n_v = 4 ; n_fy = 4
101 
102  case(ode_bdfsbdf5)
103  n_v = 5 ; n_fy = 5
104 
105  case default
106  call quit( "ode_SL_ms_mod: memSize_ode_SL_ms: error")
107 
108  end select
109 
110  end subroutine memsize_ode_sl_ms
111 
112 
115  subroutine set_solver_ode_sl_ms(slv, method)
116  integer, intent(in) :: method
117  procedure(ode_lin_solver), pointer :: slv
118 
119  select case(method)
120  case(ode_fbe)
121  slv => sl_ms_fbe
122  case(ode_fe)
123  slv => sl_ms_fe
124  case(ode_bdfsbdf2)
125  slv => sl_ms_bdfsbdf2
126  case(ode_bdfsbdf3)
127  slv => sl_ms_bdfsbdf3
128  case(ode_bdfsbdf4)
129  slv => sl_ms_bdfsbdf4
130  case(ode_bdfsbdf5)
131  slv => sl_ms_bdfsbdf5
132  case(ode_cnab2)
133  slv => sl_ms_cnab2
134  case(ode_mcnab2)
135  slv => sl_ms_mcnab2
136 
137  case default
138  call quit( "ode_SL_ms_mod:&
139  & set_solver_ode_SL_ms: uncorrect method")
140  end select
141 
142  end subroutine set_solver_ode_sl_ms
143 
144 
147  subroutine create_ode_sl_ms_sol(sol, pb, method)
148  type(ode_solution), intent(inout) :: sol
149  type(ode_problem) , intent(in) :: pb
150  integer , intent(in) :: method
151 
152  integer :: n_V, n_FY
153  logical :: bool
154 
155  bool = check_ode_method_sl_ms(method)
156  if (.NOT.bool) call quit(&
157  & "ode_SL_ms_mod: create_ode_SL_ms_sol: uncorrect method")
158 
159  if (pb%N/=1) call quit(&
160  & "ode_SL_ms_mod: create_ode_SL_ms_sol: wrong problem size")
161 
162  call memsize_ode_sl_ms(n_v, n_fy, method)
163  sol = ode_solution(pb, nv=n_v, nfy=n_fy, ny=1)
164 
165  end subroutine create_ode_sl_ms_sol
166 
167 
171  subroutine solve_ode_sl_ms(sol, pb, t0, T, dt, method, &
172  & out, check_overflow, Kinv, kry)
173  type(ode_solution) , intent(inout) :: sol
174  type(ode_problem) , intent(in) :: pb
175  real(RP) , intent(in) :: t0, T, dt
176  integer , intent(in) :: method
177  procedure(ode_output_proc) :: out
178  logical , intent(in) :: check_overflow
179  procedure(linsystem_solver), optional :: Kinv
180  type(krylov) , optional :: kry
181 
182  procedure(ode_lin_solver) , pointer :: slv=>null()
183  procedure(linsystem_solver), pointer :: KInv_1
184  type(krylov) :: kry_1
185  real(RP) :: tn, CS
186  real(DP) :: t_SL, t_out, t_reac, t_tot, cpu
187  logical :: ierr, exit_comp
188  integer :: ns, ii, jj, n0, n0_F, n0_V
189  real(RP), dimension(0) :: ay0
190 
191  if (choral_verb>1) write(*,*) &
192  &"ode_SL_ms_mod : solve_ode_SL_ms"
193  if (choral_verb>2) write(*,*) &
194  &" SemiLin multistep solver = ",&
195  & name_ode_method(method)
196  if (choral_verb>2) write(*,*) &
197  &" K_inv provided = ",&
198  & present(kinv)
199  if (.NOT.present(kinv)) then
200  if (choral_verb>2) write(*,*) &
201  &" Krylov settings provided = ",&
202  & present(kry)
203  end if
204 
205  !! Set CS to define K = M + CS*S
206  !!
207  cs = s_prefactor(method, dt)
208 
209  !! set the linear system solver for the system K*x = rhs
210  !!
211  if (present(kinv)) then
212  kinv_1 => kinv
213  else
214  kinv_1 => kinv_default
215  if (present(kry)) kry_1 = kry
216  end if
217 
218  !! set the elementary solver
219  !!
220  call set_solver_ode_sl_ms(slv, method)
221 
222  !! initialise the solution indexes
223  !!
224  call ode_solution_init_indexes(sol)
225 
226  !! timer initialisation
227  t_sl = 0.0_dp
228  t_out = 0.0_dp
229  t_reac = .0_dp
230  t_tot = clock()
231 
232  !! ode resolution
233  !!
234  ns = int( (t-t0) / dt)
235  exit_comp = .false.
236  sol%ierr = 0
237 
238  do jj=1, ns
239  tn = t0 + re(jj-1)*dt
240 
241  ! output
242  cpu = clock()
243  call out(tn, sol, exit_comp)
244  if (exit_comp) then
245  if (choral_verb>2) write(*,*) &
246  &"ode_SL_ms_mod : solve&
247  & time = ", tn, ': EXIT_COMPp'
248  return
249  end if
250  t_out = t_out + clock() - cpu
251 
252  !! perform one time-step
253  !!
254  cpu = clock()
255  call slv(sol, ierr, dt, pb, kinv_1)
256  t_sl = t_sl + clock() - cpu
257 
258  !! check linear system resolution
259  !!
260  if (ierr) then
261  sol%ierr = 1
262  if (choral_verb>1) write(*,*) &
263  &"ode_SL_ms_mod : solve,&
264  & time = ", tn, ': linear system solver error'
265  return
266  end if
267 
268  !! updates
269  !!
270  cpu = clock()
271  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY, &
272  & x=>pb%X, v=>sol%V)
273 
274  call circperm(sol%Y_i)
275  call circperm(sol%F_i)
276  call circperm(sol%V_i)
277 
278  n0 = sol%Y_i(1)
279  n0_f = sol%F_i(1)
280  n0_v = sol%V_i(1)
281 
282  ! 1. recast V into Y (note that pb%N = 1 here)
283  !
284  ! 2. compute reaction terms
285  !
286  if (pb%Na==0) then
287  !$OMP PARALLEL PRIVATE(ay0)
288  !$OMP DO
289  do ii=1, sol%dof
290  y(1, n0, ii) = v(ii, n0_v)
291 
292  call pb%AB(ay0, by(:, n0_f,ii), &
293  & x(:,ii), tn+dt, y(:,n0,ii), 1, 0 )
294  end do
295  !$OMP END DO
296  !$OMP END PARALLEL
297  else
298  !$OMP PARALLEL
299  !$OMP DO
300  do ii=1, sol%dof
301  y(1, n0, ii) = v(ii, n0_v)
302 
303  call pb%AB( ay(:,n0_f,ii), by(:,n0_f,ii), &
304  & x(:,ii), tn+dt, y(:,n0,ii), 1, 1 )
305 
306  by(1,n0_f, ii) = by(1,n0_f, ii) &
307  & + ay(1,n0_f, ii)*v(ii, n0_v)
308  end do
309  !$OMP END DO
310  !$OMP END PARALLEL
311  end if
312 
313  !! Check overflow
314  !!
315  if (check_overflow) then
316  ierr = overflow(sol%Y(:,n0,:))
317  if (ierr) then
318  sol%ierr = 10
319  if (choral_verb>2) write(*,*) &
320  &"ode_NL_ms_mod : solve,&
321  & time = ", tn, ': OVERFLOW Y'
322  return
323  end if
324 
325  ierr = overflow(sol%BY(:,n0_f,:))
326  if (ierr) then
327  sol%ierr = 10
328  if (choral_verb>2) write(*,*) &
329  &"ode_NL_ms_mod : solve,&
330  & time = ", tn, ': OVERFLOW BY'
331  return
332  end if
333 
334  if (pb%Na>0) then
335  ierr = overflow(sol%AY(:,n0_f,:))
336  if (ierr) then
337  sol%ierr = 10
338  if (choral_verb>2) write(*,*) &
339  &"ode_NL_ms_mod : solve,&
340  & time = ", tn, ': OVERFLOW AY'
341  return
342  end if
343  end if
344  end if
345 
346  end associate
347  t_reac = t_reac + clock() - cpu
348 
349  end do
350 
351  if (choral_verb>1) then
352  write(*,*)"ode_SL_NL_mod : solve, end = timing"
353  write(*,*)" Semilin. eq. integration =", &
354  & real(t_sl, sp)
355  write(*,*)" Reaction term evaluation =", &
356  & real(t_reac, sp)
357  write(*,*)" Output =", &
358  & real(t_out, sp)
359  t_tot = clock() - t_tot
360  write(*,*)" Total CPU =", &
361  & real(t_tot, sp)
362 
363  end if
364 
365  contains
366 
369  subroutine kinv_default(x, bool, b)
370  real(RP), dimension(:), intent(inout) :: x
371  logical , intent(out) :: bool
372  real(RP), dimension(:), intent(in) :: b
373 
374  call solve(x, kry_1, b, k_default)
375  bool = kry_1%ierr
376 
377  end subroutine kinv_default
378 
379  !! matrix-vector product x --> K*x, K = M + CS*S
380  !!
381  subroutine k_default(y, x)
382  real(RP), dimension(:), intent(out) :: y
383  real(RP), dimension(:), intent(in) :: x
384 
385  call pb%M(y , x)
386  call pb%S(sol%aux, x)
387  call xpay(y, cs, sol%aux)
388 
389  end subroutine k_default
390 
391  end subroutine solve_ode_sl_ms
392 
393 
396  subroutine sl_ms_fbe(sol, ierr, dt, pb, KInv)
397  type(ode_solution) , intent(inout) :: sol
398  logical , intent(out) :: ierr
399  real(RP) , intent(in) :: dt
400  type(ode_problem) , intent(in) :: pb
401  procedure(linsystem_solver) :: KInv
402 
403  integer :: n0_V, nl_V, n0_F
404 
405  n0_v = sol%V_i(1)
406  nl_v = sol%V_i(sol%nV)
407  n0_f = sol%F_i(1)
408 
409  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
410 
411  call xpay(aux, v(:, n0_v), dt, by(pb%N, n0_f, :))
412  call pb%M(rhs, aux)
413 
414  !! initial guess
415  if (n0_v /= nl_v) then
416  call copy(v(:,nl_v), v(:,n0_v))
417  end if
418 
419  !! solve the linear system
420  call kinv(v(:,nl_v), ierr, rhs)
421 
422  end associate
423 
424  end subroutine sl_ms_fbe
425 
426 
429  subroutine sl_ms_fe(sol, ierr, dt, pb, KInv)
430  type(ode_solution) , intent(inout) :: sol
431  logical , intent(out) :: ierr
432  real(RP) , intent(in) :: dt
433  type(ode_problem) , intent(in) :: pb
434  procedure(linsystem_solver) :: KInv
435 
436  integer :: n0_V, nl_V, n0_F
437 
438  n0_v = sol%V_i(1)
439  nl_v = sol%V_i(sol%nV)
440  n0_f = sol%F_i(1)
441 
442  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
443 
444  call xpay(aux, v(:, n0_v), dt, by(pb%N, n0_f, :))
445  call pb%M(rhs, aux)
446  call pb%S(aux, v(:, n0_v))
447  call xpay(rhs, -dt, aux)
448 
449 
450  !! initial guess
451  if (n0_v /= nl_v) then
452  call copy(v(:,nl_v), v(:,n0_v))
453  end if
454 
455  !! solve the linear system
456  call kinv(v(:,nl_v), ierr, rhs)
457 
458  end associate
459 
460  end subroutine sl_ms_fe
461 
462 
465  subroutine sl_ms_cnab2(sol, ierr, dt, pb, KInv)
466  type(ode_solution) , intent(inout) :: sol
467  logical , intent(out) :: ierr
468  real(RP) , intent(in) :: dt
469  type(ode_problem) , intent(in) :: pb
470  procedure(linsystem_solver) :: KInv
471 
472  integer :: ii, N, n0_V, nl_V, n0_F, n1_F
473  real(RP) :: h3_2, h1_2
474 
475  h3_2 = dt * 1.5_rp
476  h1_2 = dt * 0.5_rp
477 
478  n = pb%N
479  n0_v = sol%V_i(1)
480  nl_v = sol%V_i(sol%nV)
481  n0_f = sol%F_i(1)
482  n1_f = sol%F_i(2)
483 
484  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
485 
486  !$OMP PARALLEL
487  !$OMP DO
488  do ii=1, sol%dof
489 
490  aux(ii) = v(ii,n0_v) &
491  & + h3_2 * by(n,n0_f,ii) - h1_2 * by(n,n1_f,ii)
492 
493  end do
494  !$OMP END DO
495  !$OMP END PARALLEL
496 
497  call pb%M(rhs, aux)
498  call pb%S(aux, v(:,n0_v))
499  call xpay(rhs, -h1_2, aux)
500 
501  !! initial guess
502  if (sol%NV>1) call copy( v(:,nl_v), v(:,n0_v))
503 
504  !! solve the linear system
505  call kinv(v(:,nl_v), ierr, rhs)
506 
507  end associate
508 
509  end subroutine sl_ms_cnab2
510 
511 
514  subroutine sl_ms_mcnab2(sol, ierr, dt, pb, KInv)
515  type(ode_solution) , intent(inout) :: sol
516  logical , intent(out) :: ierr
517  real(RP) , intent(in) :: dt
518  type(ode_problem) , intent(in) :: pb
519  procedure(linsystem_solver) :: KInv
520 
521  integer :: ii, N, n0_V, n1_V, nl_V, n0_F, n1_F
522  real(RP) :: h3_2, h1_2
523 
524  h3_2 = dt * 1.5_rp
525  h1_2 = dt * 0.5_rp
526 
527  n = pb%N
528  n0_v = sol%V_i(1)
529  n1_v = sol%V_i(2)
530  nl_v = sol%V_i(sol%nV)
531  n0_f = sol%F_i(1)
532  n1_f = sol%F_i(2)
533 
534  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
535 
536  !$OMP PARALLEL
537  !$OMP DO
538  do ii=1, sol%dof
539 
540  aux(ii) = v(ii,n0_v) &
541  & + h3_2 * by(n,n0_f,ii) - h1_2 * by(n,n1_f,ii)
542 
543  end do
544  !$OMP END DO
545  !$OMP END PARALLEL
546 
547  call pb%M(rhs, aux)
548  call pb%S(aux, v(:,n0_v))
549  call xpay(rhs, -re(3,8)*dt, aux)
550  call pb%S(aux, v(:,n1_v))
551  call xpay(rhs, -re(1,16)*dt, aux)
552 
553  !! initial guess
554  call copy( v(:,nl_v), v(:,n0_v))
555 
556  !! solve the linear system
557  call kinv(v(:,nl_v), ierr, rhs)
558 
559  end associate
560 
561  end subroutine sl_ms_mcnab2
562 
563 
564  ! call copy(aux, V(:,n0_V))
565  ! call add_aux_AY_p_B(1, dt*1.5_RP)
566  ! call add_aux_AY_p_B(2,-dt/2._RP)
567 
568  ! call M(rhs, aux)
569  ! call S(aux, V(:,n0_V))
570  ! call xpay(rhs, -re(3,8)*dt, aux)
571  ! call S(aux, V(:,n1_V))
572  ! call xpay(rhs, -re(1,16)*dt, aux)
573 
574  ! !! initial guess
575  ! if (size(V)>1) call copy(V(NV), V(:,n0_V))
576 
577  ! !! solve the linear system
578  ! call KInv(V(NV), ierr, rhs)
579 
580 
581 
584  subroutine sl_ms_bdfsbdf2(sol, ierr, dt, pb, KInv)
585  type(ode_solution) , intent(inout) :: sol
586  logical , intent(out) :: ierr
587  real(RP) , intent(in) :: dt
588  type(ode_problem) , intent(in) :: pb
589  procedure(linsystem_solver) :: KInv
590 
591  integer :: ii, N, n0_V, n1_V, nl_V, n0_F, n1_F
592  real(RP) :: h4_3, h2_3
593 
594  h4_3 = dt * c4_3
595  h2_3 = dt * c2_3
596 
597  n = pb%N
598  n0_v = sol%V_i(1)
599  n1_v = sol%V_i(2)
600  nl_v = sol%V_i(sol%nV)
601  n0_f = sol%F_i(1)
602  n1_f = sol%F_i(2)
603 
604  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
605 
606  !$OMP PARALLEL
607  !$OMP DO
608  do ii=1, sol%dof
609 
610  aux(ii) = c4_3 * v(ii,n0_v) - c_3 * v(ii,n1_v) &
611  & + h4_3 * by(n,n0_f,ii) - h2_3 * by(n,n1_f,ii)
612 
613  end do
614  !$OMP END DO
615  !$OMP END PARALLEL
616 
617  call pb%M(rhs, aux)
618 
619  !! initial guess
620  call copy( v(:,nl_v), v(:,n0_v))
621 
622  !! solve the linear system
623  call kinv(v(:,nl_v), ierr, rhs)
624 
625  end associate
626 
627  end subroutine sl_ms_bdfsbdf2
628 
629 
632  subroutine sl_ms_bdfsbdf3(sol, ierr, dt, pb, KInv)
633  type(ode_solution) , intent(inout) :: sol
634  logical , intent(out) :: ierr
635  real(RP) , intent(in) :: dt
636  type(ode_problem) , intent(in) :: pb
637  procedure(linsystem_solver) :: KInv
638 
639  integer :: ii, N, n0_V, n1_V, n2_V, nl_V, n0_F, n1_F, n2_F
640  real(RP) :: a, h18_11, h6_11
641 
642  h18_11 = dt * c18_11
643  h6_11 = dt * c6_11
644 
645  n = pb%N
646  n0_v = sol%V_i(1)
647  n1_v = sol%V_i(2)
648  n2_v = sol%V_i(3)
649  nl_v = sol%V_i(sol%nV)
650  n0_f = sol%F_i(1)
651  n1_f = sol%F_i(2)
652  n2_f = sol%F_i(3)
653 
654  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
655 
656  !$OMP PARALLEL PRIVATE(a)
657  !$OMP DO
658  do ii=1, sol%dof
659 
660  a = c18_11 * v(ii,n0_v) - c9_11 * v(ii,n1_v) &
661  & + c2_11 * v(ii,n2_v)
662 
663  a = a + h18_11 * by(n,n0_f,ii) - h18_11 * by(n,n1_f,ii) &
664  & + h6_11 * by(n,n2_f,ii)
665 
666  aux(ii) = a
667 
668  end do
669  !$OMP END DO
670  !$OMP END PARALLEL
671 
672  call pb%M(rhs, aux)
673 
674  !! initial guess
675  call copy( v(:,nl_v), v(:,n0_v))
676 
677  !! solve the linear system
678  call kinv(v(:,nl_v), ierr, rhs)
679 
680  end associate
681 
682  end subroutine sl_ms_bdfsbdf3
683 
684 
685 
688  subroutine sl_ms_bdfsbdf4(sol, ierr, dt, pb, KInv)
689  type(ode_solution) , intent(inout) :: sol
690  logical , intent(out) :: ierr
691  real(RP) , intent(in) :: dt
692  type(ode_problem) , intent(in) :: pb
693  procedure(linsystem_solver) :: KInv
694 
695  integer :: ii, N, n0_V, n1_V, n2_V, nl_V, n0_F, n1_F, n2_F
696  integer :: n3_V, n3_F
697  real(RP) :: a, h12_25, h48_25, h72_25
698 
699  h12_25 = dt * c12_25
700  h48_25 = dt * c48_25
701  h72_25 = dt * c72_25
702 
703  n = pb%N
704  n0_v = sol%V_i(1)
705  n1_v = sol%V_i(2)
706  n2_v = sol%V_i(3)
707  n3_v = sol%V_i(4)
708  nl_v = sol%V_i(sol%nV)
709  n0_f = sol%F_i(1)
710  n1_f = sol%F_i(2)
711  n2_f = sol%F_i(3)
712  n3_f = sol%F_i(4)
713 
714  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
715 
716  !$OMP PARALLEL PRIVATE(a)
717  !$OMP DO
718  do ii=1, sol%dof
719 
720  a = c48_25 * v(ii,n0_v) - c36_25 * v(ii,n1_v) &
721  & + c16_25 * v(ii,n2_v) - c3_25 * v(ii,n3_v)
722 
723  a = a + h48_25 * by(n,n0_f,ii) - h72_25 * by(n,n1_f,ii) &
724  & + h48_25 * by(n,n2_f,ii) - h12_25 * by(n,n3_f,ii)
725 
726  aux(ii) = a
727 
728  end do
729  !$OMP END DO
730  !$OMP END PARALLEL
731 
732  call pb%M(rhs, aux)
733 
734  !! initial guess
735  call copy( v(:,nl_v), v(:,n0_v))
736 
737  !! solve the linear system
738  call kinv(v(:,nl_v), ierr, rhs)
739 
740  end associate
741 
742  end subroutine sl_ms_bdfsbdf4
743 
744 
747  subroutine sl_ms_bdfsbdf5(sol, ierr, dt, pb, KInv)
748  type(ode_solution) , intent(inout) :: sol
749  logical , intent(out) :: ierr
750  real(RP) , intent(in) :: dt
751  type(ode_problem) , intent(in) :: pb
752  procedure(linsystem_solver) :: KInv
753 
754  integer :: ii, N, n0_V, n1_V, n2_V, nl_V, n0_F, n1_F, n2_F
755  integer :: n3_V, n3_F, n4_V, n4_F
756  real(RP) :: a, h60_137, h300_137, h600_137
757 
758  h60_137 = dt * c60_137
759  h300_137 = dt * c300_137
760  h600_137 = dt * c600_137
761 
762  n = pb%N
763  n0_v = sol%V_i(1)
764  n1_v = sol%V_i(2)
765  n2_v = sol%V_i(3)
766  n3_v = sol%V_i(4)
767  n4_v = sol%V_i(5)
768  nl_v = sol%V_i(sol%nV)
769  n0_f = sol%F_i(1)
770  n1_f = sol%F_i(2)
771  n2_f = sol%F_i(3)
772  n3_f = sol%F_i(4)
773  n4_f = sol%F_i(5)
774 
775  associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
776 
777  !$OMP PARALLEL PRIVATE(a)
778  !$OMP DO
779  do ii=1, sol%dof
780 
781  a = c300_137 * v(ii,n0_v) - c300_137 * v(ii,n1_v) &
782  & + c200_137 * v(ii,n2_v) - c75_137 * v(ii,n3_v) &
783  & + c12_137 * v(ii,n4_v)
784 
785  a = a + h300_137* by(n,n0_f,ii) - h600_137* by(n,n1_f,ii) &
786  & + h600_137* by(n,n2_f,ii) - h300_137* by(n,n3_f,ii) &
787  & + h60_137 * by(n,n4_f,ii)
788 
789  aux(ii) = a
790 
791  end do
792  !$OMP END DO
793  !$OMP END PARALLEL
794 
795  call pb%M(rhs, aux)
796 
797  !! initial guess
798  call copy( v(:,nl_v), v(:,n0_v))
799 
800  !! solve the linear system
801  call kinv(v(:,nl_v), ierr, rhs)
802 
803  end associate
804 
805  end subroutine sl_ms_bdfsbdf5
806 
807 
808 end module ode_sl_ms_mod
Type ode_solution: data structure to solve ODE/PDE problems.
BASIC TOOLS
Definition: basic_tools.F90:8
integer, parameter ode_mcnab2
ModifiedCrank Nicholson / Adamns Bashforth 2.
double precision function, public clock()
Returns the internal clock time
Definition: basic_tools.F90:94
integer, parameter ode_bdfsbdf3
BDF / SBDF 3.
subroutine k_default(y, x)
DERIVED TYPE ode_problem: definition of ODE/PDE problems
real(rp), parameter c2_3
subroutine sl_ms_bdfsbdf2(sol, ierr, dt, pb, KInv)
BDF2/SBDF2.
real(rp), parameter c36_25
real(rp) function, public s_prefactor(method, dt)
When discretising with two matrices, this function returns the prefactor Cs for the matrix S...
Definition: ode_def.f90:205
real(rp), parameter c4_3
real(rp), parameter c12_25
The type krylov defines the settings of a linear solver.
Definition: krylov_mod.f90:63
subroutine, public set_solver_ode_sl_ms(slv, method)
set the resolution solver
real(rp), parameter c300_137
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine, public create_ode_sl_ms_sol(sol, pb, method)
create memory for the ode_solution structure 'sol'
integer, parameter, public sp
simple precision for real numbers
Definition: real_type.F90:60
subroutine sl_ms_bdfsbdf5(sol, ierr, dt, pb, KInv)
BDF5/SBDF5.
logical function, public overflow(yy)
Detects overflow.
Definition: ode_def.f90:254
character(len=15) function, public name_ode_method(method)
Get ODE method name.
Definition: ode_def.f90:68
real(rp), parameter c6_11
conversion integers or rational to real
Definition: real_type.F90:153
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
real(rp) function cs(x)
subroutine, public ode_solution_init_indexes(sol)
initialise the ode_indexes
subroutine, public memsize_ode_sl_ms(n_V, n_FY, method)
required sizes to allocate memory
integer, parameter ode_fbe
Forward / backward Euler.
real(rp), parameter c72_25
DERIVED TYPE ode_solution: data straucture to solve ODEs
subroutine sl_ms_cnab2(sol, ierr, dt, pb, KInv)
CN/AB2.
real(rp), parameter c9_11
subroutine sl_ms_bdfsbdf4(sol, ierr, dt, pb, KInv)
BDF4/SBDF4.
y = x
Definition: R1d_mod.F90:44
subroutine kinv_default(x, bool, b)
Default solver for K*x = b.
subroutine sl_ms_mcnab2(sol, ierr, dt, pb, KInv)
Modified-CN/AB2.
CHORAL CONSTANTS
real(rp), parameter c48_25
real(rp), parameter c3_25
ALGEBRAIC OPERATIONS ON SETS
Definition: algebra_set.F90:14
real(rp), parameter c16_25
logical function, public check_ode_method_sl_ms(method)
is 'method' a multi-step semilinear ODE solver ?
real(rp), parameter c200_137
real(rp), parameter c18_11
subroutine sl_ms_fbe(sol, ierr, dt, pb, KInv)
FORWARD/BacWard Euler.
real(rp), parameter c_3
integer, parameter ode_bdfsbdf5
BDF / SBDF 5.
BOTTOM LEVEL MODULE FOR ODEs
Definition: ode_def.f90:12
integer choral_verb
Verbosity level.
real(rp), parameter c12_137
x = x + b*y // OR // z = x + b*y
Definition: R1d_mod.F90:60
real(rp), parameter c2_11
z = a*x + b*y
Definition: R1d_mod.F90:66
subroutine sl_ms_fe(sol, ierr, dt, pb, KInv)
FORWARD Euler.
integer, parameter ode_bdfsbdf4
BDF / SBDF 4.
subroutine, public circperm(E)
Circular permutation of an array of integer.
MULTISTEP SOLVERS FOR SEMILINEAR ODEs
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
subroutine, public solve_ode_sl_ms(sol, pb, t0, T, dt, method, out, check_overflow, Kinv, kry)
solve with constant time-step Case where pbN = 1
integer, parameter ode_cnab2
Crank Nicholson / Adamns Bashforth 2.
subroutine sl_ms_bdfsbdf3(sol, ierr, dt, pb, KInv)
BDF3/SBDF3.
integer, parameter ode_fe
Forward Euler.
integer, parameter ode_bdfsbdf2
BDF / SBDF 2.
Type ode_problem: definition of ODE/PDE problems.
real(rp), parameter c60_137
DERIVED TYPE krylov: for the resolution of linear systems
Definition: krylov_mod.f90:25
real(rp), parameter c75_137
real(rp), parameter c600_137