Choral
ode_NL_ms_mod.f90
Go to the documentation of this file.
1 
9 
10 module ode_nl_ms_mod
11 
12  use choral_constants
13  use choral_variables
14  use real_type
15  use basic_tools
16  use algebra_set
17  use ode_def
18  use ode_problem_mod
19  use ode_solution_mod
20  !$ use OMP_LIB
21 
22  implicit none
23  private
24 
25  public :: solve_ode_nl_ms
26  public :: create_ode_nl_ms_sol
27  public :: set_solver_ode_nl_ms
28  public :: memsize_ode_nl_ms
29  public :: check_ode_method_nl_ms
30 
31  real(RP), parameter :: c1_3 = 1._rp / 3._rp
32  real(RP), parameter :: c2_3 = 2._rp / 3._rp
33  real(RP), parameter :: c4_3 = 4._rp / 3._rp
34 
35  real(RP), parameter :: c1_6 = 1._rp / 6._rp
36  real(RP), parameter :: c11_6 = 11._rp / 6._rp
37 
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 
43  real(RP), parameter :: c5_12 = 5._rp / 12._rp
44  real(RP), parameter :: c16_12 = 16._rp / 12._rp
45  real(RP), parameter :: c23_12 = 23._rp / 12._rp
46 
47  real(RP), parameter :: c9_16 = 9._rp / 16._rp
48 
49  real(RP), parameter :: c1_24 = 1._rp / 24._rp
50  real(RP), parameter :: c9_24 = 9._rp / 24._rp
51  real(RP), parameter :: c37_24 = 37._rp / 24._rp
52  real(RP), parameter :: c55_24 = 55._rp / 24._rp
53  real(RP), parameter :: c59_24 = 59._rp / 24._rp
54 
55  real(RP), parameter :: c3_25 = 3._rp / 25._rp
56  real(RP), parameter :: c12_25 = 12._rp / 25._rp
57  real(RP), parameter :: c16_25 = 16._rp / 25._rp
58  real(RP), parameter :: c36_25 = 36._rp / 25._rp
59  real(RP), parameter :: c48_25 = 48._rp / 25._rp
60  real(RP), parameter :: c72_25 = 72._rp / 25._rp
61 
62  real(RP), parameter :: c12_137 = 12._rp / 137._rp
63  real(RP), parameter :: c60_137 = 60._rp / 137._rp
64  real(RP), parameter :: c75_137 = 75._rp / 137._rp
65  real(RP), parameter :: c200_137 = 200._rp / 137._rp
66  real(RP), parameter :: c300_137 = 300._rp / 137._rp
67  real(RP), parameter :: c600_137 = 600._rp / 137._rp
68 
69 contains
70 
72  function check_ode_method_nl_ms(method) result(b)
73  logical :: b
74  integer, intent(in) :: method
75 
76  b = .false.
77 
78  select case(method)
79  case(ode_fbe, ode_erk1, ode_rl2, ode_rl3, ode_rl4, &
80  & ode_bdfsbdf2, ode_bdfsbdf3, ode_bdfsbdf4, &
81  & ode_bdfsbdf5, ode_cnab2, ode_mcnab2, &
82  & ode_eab2, ode_eab3, ode_eab4)
83 
84  b = .true.
85  end select
86 
87  end function check_ode_method_nl_ms
88 
89 
97  subroutine memsize_ode_nl_ms(n_Y, n_FY, method)
98  integer, intent(out) :: n_Y, n_FY
99  integer, intent(in) :: method
100 
101  select case(method)
102 
103  case(ode_fbe, ode_erk1)
104  n_y = 1 ; n_fy = 1
105 
106  case(ode_rl2)
107  n_y = 1 ; n_fy = 2
108 
109  case(ode_rl3)
110  n_y = 1 ; n_fy = 3
111 
112  case(ode_rl4)
113  n_y = 1 ; n_fy = 4
114 
115  case(ode_eab2, ode_bdfsbdf2, ode_cnab2, ode_mcnab2)
116  n_y = 2 ; n_fy = 2
117 
118  case(ode_bdfsbdf3)
119  n_y = 3 ; n_fy = 3
120 
121  case(ode_bdfsbdf4)
122  n_y = 4 ; n_fy = 4
123 
124  case(ode_bdfsbdf5)
125  n_y = 5 ; n_fy = 5
126 
127  case(ode_eab3)
128  n_y = 3 ; n_fy = 3
129 
130  case(ode_eab4)
131  n_y = 4 ; n_fy = 4
132 
133  case default
134  call quit( "ode_NL_ms_mod: memSize_ode_NL_ms:&
135  & uncorrect method")
136 
137  end select
138 
139  end subroutine memsize_ode_nl_ms
140 
141 
144  subroutine set_solver_ode_nl_ms(slv, method)
145  integer, intent(in) :: method
146  procedure(ode_nl_ms_solver), pointer :: slv
147 
148  select case(method)
149  case(ode_fbe)
150  slv => nl_ms_fbe
151  case(ode_erk1)
152  slv => nl_ms_erk1
153  case(ode_rl2)
154  slv => nl_ms_rl2
155  case(ode_rl3)
156  slv => nl_ms_rl3
157  case(ode_rl4)
158  slv => nl_ms_rl4
159  case(ode_eab2)
160  slv => nl_ms_eab2
161  case(ode_eab3)
162  slv => nl_ms_eab3
163  case(ode_eab4)
164  slv => nl_ms_eab4
165  case(ode_bdfsbdf2)
166  slv => nl_ms_bdfsbdf2
167  case(ode_bdfsbdf3)
168  slv => nl_ms_bdfsbdf3
169  case(ode_bdfsbdf4)
170  slv => nl_ms_bdfsbdf4
171  case(ode_bdfsbdf5)
172  slv => nl_ms_bdfsbdf5
173  case(ode_cnab2)
174  slv => nl_ms_cnab2
175  case(ode_mcnab2)
176  slv => nl_ms_mcnab2
177 
178  case default
179  call quit( "ode_NL_ms_mod: set_solver_ode_NL_ms:&
180  & uncorrect method")
181  end select
182 
183  end subroutine set_solver_ode_nl_ms
184 
187  subroutine create_ode_nl_ms_sol(sol, pb, method)
188  type(ode_solution), intent(inout) :: sol
189  type(ode_problem) , intent(in) :: pb
190  integer , intent(in) :: method
191 
192  integer :: n_Y, n_FY
193  logical :: bool
194 
195  bool = check_ode_method_nl_ms(method)
196  if (.NOT.bool) call quit(&
197  & "ode_NL_ms_mod: create_ode_NL_ms_sol: uncorrect method")
198 
199  call memsize_ode_nl_ms(n_y, n_fy, method)
200 
201  sol = ode_solution(pb, ny=n_y, nfy=n_fy)
202 
203  end subroutine create_ode_nl_ms_sol
204 
205 
208  subroutine solve_ode_nl_ms(sol, pb, t0, T, dt, method, out, &
209  & check_overflow)
210  type(ode_solution), intent(inout) :: sol
211  type(ode_problem) , intent(in) :: pb
212  real(RP) , intent(in) :: t0, T, dt
213  integer , intent(in) :: method
214  procedure(ode_output_proc) :: out
215  logical , intent(in) :: check_overflow
216 
217 
218  procedure(ode_nl_ms_solver), pointer :: slv=>null()
219  real(RP) :: tn
220  integer :: ns, ii, jj, n0, n0_F
221  logical :: exit_comp, ierr
222 
223  if (choral_verb>1) write(*,*) &
224  &"ode_NL_ms_mod : solve_ode_NL_ms"
225  if (choral_verb>2) write(*,*) &
226  &" NonLin multistep solver = ",&
227  & name_ode_method(method)
228 
229  !! set the elementary solver
230  !!
231  call set_solver_ode_nl_ms(slv, method)
232 
233 
234  !! initialise the solution indexes
235  !!
236  call ode_solution_init_indexes(sol)
237 
238  !! ode resolution
239  !!
240  exit_comp = .false.
241  sol%ierr = 0
242  ns = int( (t-t0) / dt)
243 
244  do jj=1, ns
245  tn = t0 + re(jj-1)*dt
246 
247  ! output
248  call out(tn, sol, exit_comp)
249  if (exit_comp) then
250  if (choral_verb>2) write(*,*) &
251  &"ode_NL_ms_mod : solve&
252  & time = ", tn, ': EXIT_COMPp'
253  return
254  end if
255 
256  !! perform one time-step
257  !!
258  call slv(sol, dt, pb%N, pb%Na)
259 
260  !! updates
261  !!
262  call circperm(sol%Y_i)
263  call circperm(sol%F_i)
264 
265  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY, &
266  & x=>pb%X)
267 
268  n0 = sol%Y_i(1)
269  n0_f = sol%F_i(1)
270 
271  !$OMP PARALLEL
272  !$OMP DO
273  do ii=1, sol%dof
274  call pb%AB( ay(:,n0_f, ii), by(:,n0_f, ii), &
275  & x(:,ii), tn+dt, y(:,n0, ii), pb%N, pb%Na )
276  end do
277  !$OMP END DO
278  !$OMP END PARALLEL
279 
280  !! Check overflow
281  !!
282  if (check_overflow) then
283  ierr = overflow(sol%Y(:,n0,:))
284  if (ierr) then
285  sol%ierr = 10
286  if (choral_verb>1) write(*,*) &
287  &"ode_NL_ms_mod : solve,&
288  & time = ", tn, ': OVERFLOW Y'
289  return
290  end if
291 
292  ierr = overflow(sol%BY(:,n0_f,:))
293  if (ierr) then
294  sol%ierr = 10
295  if (choral_verb>1) write(*,*) &
296  &"ode_NL_ms_mod : solve,&
297  & time = ", tn, ': OVERFLOW BY'
298  return
299  end if
300 
301  ierr = overflow(sol%AY(:,n0_f,:))
302  if (ierr) then
303  sol%ierr = 10
304  if (choral_verb>1) write(*,*) &
305  &"ode_NL_ms_mod : solve,&
306  & time = ", tn, ': OVERFLOW AY'
307  return
308  end if
309  end if
310 
311  end associate
312 
313  end do
314 
315  end subroutine solve_ode_nl_ms
316 
317 
318 
321  elemental subroutine phi_1(z)
322  real(RP), intent(inout) :: z
323 
324  if (abs(z)> real_tol) then
325  z = (exp(z) - 1._rp)/z
326  else
327  z = 1._rp
328  end if
329 
330  end subroutine phi_1
331 
332 
335  elemental subroutine phi_2(z)
336  real(RP), intent(inout) :: z
337 
338  real(RP) :: z2
339 
340  z2 = z**2
341 
342  if (abs(z2)>real_tol) then
343  z = (exp(z) - z - 1._rp)/z2
344 
345  else
346  z = 0.5_rp
347 
348  end if
349 
350  end subroutine phi_2
351 
352 
355  elemental subroutine phi_3(z)
356  real(RP), intent(inout) :: z
357 
358  real(RP) :: z2, z3
359 
360  z2 = z**2
361  z3 = z2*z
362 
363  if (abs(z3)>real_tol) then
364  z = (exp(z) - 0.5_rp*z2 - z - 1._rp)/z3
365 
366  else
367  z = c1_6
368 
369  end if
370 
371  end subroutine phi_3
372 
375  elemental subroutine phi_4(z)
376  real(RP), intent(inout) :: z
377 
378  real(RP) :: z2, z3, z4
379 
380  z2 = z *z
381  z3 = z2*z
382  z4 = z3*z
383 
384  if (abs(z4)>real_tol) then
385  z = (exp(z) - c1_6*z3 - 0.5_rp*z2 - z - 1._rp)/z4
386 
387  else
388  z = c1_24
389 
390  end if
391  end subroutine phi_4
392 
393 
397  subroutine nl_ms_fbe(sol, dt, N, P)
398  type(ode_solution), intent(inout) :: sol
399  real(RP) , intent(in) :: dt
400  integer , intent(in) :: N, P
401 
402  integer :: ii
403  integer :: n0, nl, n0_F
404 
405  n0 = sol%Y_i(1)
406  nl = sol%Y_i(sol%nY)
407  n0_f = sol%F_i(1)
408 
409  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
410 
411  !$OMP PARALLEL
412  !$OMP DO
413  do ii=1, sol%dof
414 
415  y(1:n, nl, ii) = y(1:n, n0, ii) + dt*by(1:n, n0_f, ii)
416  y(1:p, nl, ii) = y(1:p, nl, ii) / &
417  & (1._rp - dt*ay(1:p, n0_f, ii))
418  end do
419  !$OMP END DO
420  !$OMP END PARALLEL
421 
422  end associate
423 
424  end subroutine nl_ms_fbe
425 
426 
430  subroutine nl_ms_erk1(sol, dt, N, P)
431  type(ode_solution), intent(inout) :: sol
432  real(RP) , intent(in) :: dt
433  integer , intent(in) :: N, P
434 
435  real(RP), dimension(N) :: b
436  real(RP), dimension(P) :: a
437  integer :: ii
438  integer :: n0, nl, n0_F
439 
440  n0 = sol%Y_i(1)
441  nl = sol%Y_i(sol%nY)
442  n0_f = sol%F_i(1)
443 
444  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
445 
446  !$OMP PARALLEL PRIVATE(a, b)
447  !$OMP DO
448  do ii=1, sol%dof
449 
450  !! a = a_n
451  a = ay(1:p,n0_f, ii)
452 
453  !! b = b_n
454  b = by(1:n,n0_f, ii)
455 
456  !! b(1:P) := a*y_n + b_n
457  b(1:p) = a * y(1:p,n0, ii) + b(1:p)
458 
459  !! a = phi_1(a*dt)
460  a = a*dt
461  call phi_1(a)
462 
463  !! b(1:P) = a*b
464  b(1:p) = a * b(1:p)
465 
466  y(1:n,nl, ii) = y(1:n,n0, ii) + dt*b
467 
468  end do
469  !$OMP END DO
470  !$OMP END PARALLEL
471 
472  end associate
473  end subroutine nl_ms_erk1
474 
475 
479  subroutine nl_ms_rl2(sol, dt, N, P)
480  type(ode_solution), intent(inout) :: sol
481  real(RP) , intent(in) :: dt
482  integer , intent(in) :: N, P
483 
484  real(RP), dimension(N) :: b
485  real(RP), dimension(P) :: a
486  integer :: ii
487  integer :: n0, nl, n0_F, n1_F
488 
489  n0 = sol%Y_i(1)
490  nl = sol%Y_i(sol%nY)
491  n0_f = sol%F_i(1)
492  n1_f = sol%F_i(2)
493 
494  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
495 
496  !$OMP PARALLEL PRIVATE(a, b)
497  !$OMP DO
498  do ii=1, sol%dof
499 
500  !! a = a_n*1.5 - a_{n-1}*0.5
501  a = ay(1:p,n0_f, ii)*1.5_rp - ay(1:p,n1_f, ii)*0.5_rp
502 
503  !! b = b_n*1.5 - b_{n-1}*0.5
504  b = by(1:n,n0_f, ii)*1.5_rp - by(1:n,n1_f, ii)*0.5_rp
505 
506  !! b(1:P) := a*y_n + b_n
507  b(1:p) = a * y(1:p,n0, ii) + b(1:p)
508 
509  !! a = phi_1(a*dt)
510  a = a*dt
511  call phi_1(a)
512 
513  !! b(1:P) = a*b
514  b(1:p) = a * b(1:p)
515 
516  y(1:n,nl, ii) = y(1:n,n0, ii) + dt*b
517 
518  end do
519  !$OMP END DO
520  !$OMP END PARALLEL
521 
522  end associate
523  end subroutine nl_ms_rl2
524 
525 
529  subroutine nl_ms_rl3(sol, dt, N, P)
530  type(ode_solution), intent(inout) :: sol
531  real(RP) , intent(in) :: dt
532  integer , intent(in) :: N, P
533 
534  real(RP), dimension(N) :: b
535  real(RP), dimension(P) :: a
536  real(RP) :: h1
537  integer :: ii
538  integer :: n0, nl, n0_F, n1_F, n2_F
539 
540  h1 = dt / 12._rp
541 
542  n0 = sol%Y_i(1)
543  nl = sol%Y_i(sol%nY)
544  n0_f = sol%F_i(1)
545  n1_f = sol%F_i(2)
546  n2_f = sol%F_i(3)
547 
548  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
549 
550  !$OMP PARALLEL PRIVATE(a, b)
551  !$OMP DO
552  do ii=1, sol%dof
553 
554  !! a = a_n*23/15 - a_{n-1}*16/12 + a_{n-2}*5/12
555  a = ay(1:p,n0_f, ii)*c23_12 - ay(1:p,n1_f, ii)*c16_12 &
556  & + ay(1:p,n2_f, ii)*c5_12
557 
558  !! b = b_n*23/15 - b_{n-1}*16/12 + b_{n-2}*5/12
559  b = by(1:n,n0_f, ii)*c23_12 - by(1:n,n1_f, ii)*c16_12 &
560  & + by(1:n,n2_f, ii)*c5_12
561 
562  !! b(1:P) = b(1:P) + (a_n*b_{n-1} - a_{n-1}*b_n)*dt/12
563  b(1:p) = b(1:p) + h1*ay(1:p,n0_f, ii)*by(1:p,n1_f, ii)
564  b(1:p) = b(1:p) - h1*ay(1:p,n1_f, ii)*by(1:p,n0_f, ii)
565 
566  !! b(1:P) := a*y_n + b_n
567  b(1:p) = a * y(1:p,n0, ii) + b(1:p)
568 
569  !! a = phi_1(a*dt)
570  a = a*dt
571  call phi_1(a)
572 
573  !! b(1:P) = a*b
574  b(1:p) = a * b(1:p)
575 
576  y(1:n,nl, ii) = y(1:n,n0, ii) + dt*b
577 
578  end do
579  !$OMP END DO
580  !$OMP END PARALLEL
581 
582  end associate
583  end subroutine nl_ms_rl3
584 
585 
589  subroutine nl_ms_rl4(sol, dt, N, P)
590  type(ode_solution), intent(inout) :: sol
591  real(RP) , intent(in) :: dt
592  integer , intent(in) :: N, P
593 
594  real(RP), dimension(N) :: b
595  real(RP), dimension(P) :: a
596  real(RP) :: h3, h1
597  integer :: ii
598  integer :: n0, nl, n0_F, n1_F, n2_F, n3_F
599 
600  h3 = 3._rp * dt / 12._rp
601  h1 = dt / 12._rp
602 
603  n0 = sol%Y_i(1)
604  nl = sol%Y_i(sol%nY)
605  n0_f = sol%F_i(1)
606  n1_f = sol%F_i(2)
607  n2_f = sol%F_i(3)
608  n3_f = sol%F_i(4)
609 
610  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
611 
612  !$OMP PARALLEL PRIVATE(a, b)
613  !$OMP DO
614  do ii=1, sol%dof
615 
616  !! a = a_n*55/24 - a_{n-1}*559/24 + a_{n-2}*7/24
617  !! - a_{n-}*9/24
618  a = ay(1:p,n0_f, ii)*c55_24 - ay(1:p,n1_f, ii)*c59_24 &
619  & + ay(1:p,n2_f, ii)*c37_24 - ay(1:p,n3_f, ii)*c9_24
620 
621  !! b = b_n*55/24 - b_{n-1}*559/24 + b_{n-2}*7/24
622  !! - b_{n-}*9/24
623  b = by(1:n,n0_f, ii)*c55_24 - by(1:n,n1_f, ii)*c59_24 &
624  & + by(1:n,n2_f, ii)*c37_24 - by(1:n,n3_f, ii)*c9_24
625 
626  !! b(1:P) = b(1:P) + a_n*(b_{n-1}*3/12 - b_{n-2}/12)*dt
627  b(1:p) = b(1:p) + ay(1:p,n0_f, ii) * &
628  & (h3*by(1:p,n1_f, ii) - h1*by(1:p,n2_f, ii))
629 
630  !! b(1:P) = b(1:P) - b_n*(a_{n-1}*3/12 + a_{n-2}/12)*dt
631  b(1:p) = b(1:p) - by(1:p,n0_f, ii) * &
632  & (h3*ay(1:p,n1_f, ii) - h1*ay(1:p,n2_f, ii))
633 
634  !! b(1:P) := a*y_n + b
635  b(1:p) = a * y(1:p,n0, ii) + b(1:p)
636 
637  !! a = phi_1(a*dt)
638  a = a*dt
639  call phi_1(a)
640 
641  !! b(1:P) = a*b
642  b(1:p) = a * b(1:p)
643 
644  y(1:n,nl, ii) = y(1:n,n0, ii) + dt*b
645 
646  end do
647  !$OMP END DO
648  !$OMP END PARALLEL
649 
650  end associate
651  end subroutine nl_ms_rl4
652 
653 
657  subroutine nl_ms_bdfsbdf2(sol, dt, N, P)
658  type(ode_solution), intent(inout) :: sol
659  real(RP) , intent(in) :: dt
660  integer , intent(in) :: N, P
661 
662  real(RP), dimension(N) :: b
663  real(RP), dimension(P) :: a
664  real(RP) :: h2, h4
665  integer :: ii
666  integer :: n0, n1, nl, n0_F, n1_F
667 
668  h2 = dt * c2_3
669  h4 = dt * c4_3
670 
671  n0 = sol%Y_i(1)
672  n1 = sol%Y_i(2)
673  nl = sol%Y_i(sol%nY)
674  n0_f = sol%F_i(1)
675  n1_f = sol%F_i(2)
676 
677  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
678 
679  !$OMP PARALLEL PRIVATE(a, b)
680  !$OMP DO
681  do ii=1, sol%dof
682 
683  b = by(1:n, n0_f, ii)*h4 - by(1:n, n1_f, ii)*h2
684 
685  a = ay(1:p, n0_f, ii) - ay(1:p, n1_f, ii)
686  a = a * y(1:p, n1_f, ii)
687 
688  b(1:p) = b(1:p) + a * h2
689 
690  b = b + y(1:n, n0, ii)*c4_3 - y(1:n, n1, ii)*c1_3
691 
692  b(1:p) = b(1:p) / (1._rp - h2*ay(1:p, n0_f, ii))
693 
694  y(1:n, nl, ii) = b
695 
696  end do
697  !$OMP END DO
698  !$OMP END PARALLEL
699 
700  end associate
701 
702  end subroutine nl_ms_bdfsbdf2
703 
704 
705 
709  subroutine nl_ms_bdfsbdf3(sol, dt, N, P)
710  type(ode_solution), intent(inout) :: sol
711  real(RP) , intent(in) :: dt
712  integer , intent(in) :: N, P
713 
714  real(RP), dimension(N) :: b
715  real(RP), dimension(P) :: a
716  real(RP) :: h6, h18
717  integer :: ii
718  integer :: n0, n1, n2, nl, n0_F, n1_F, n2_F
719 
720  h6 = dt * c6_11
721  h18 = dt * c18_11
722 
723  n0 = sol%Y_i(1)
724  n1 = sol%Y_i(2)
725  n2 = sol%Y_i(3)
726  nl = sol%Y_i(sol%nY)
727  n0_f = sol%F_i(1)
728  n1_f = sol%F_i(2)
729  n2_f = sol%F_i(3)
730 
731  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
732 
733  !$OMP PARALLEL PRIVATE(a, b)
734  !$OMP DO
735  do ii=1, sol%dof
736 
737  b = by(1:n, n0_f, ii)*h18 - by(1:n, n1_f, ii)*h18 &
738  & + by(1:n, n2_f, ii)*h6
739 
740  a = ay(1:p, n2_f, ii) - ay(1:p, n0_f, ii)
741  a = a * y(1:p, n2_f, ii)
742  b(1:p) = b(1:p) + a * h6
743 
744  a = ay(1:p, n1_f, ii) - ay(1:p, n0_f, ii)
745  a = a * y(1:p, n1_f, ii)
746  b(1:p) = b(1:p) - a * h18
747 
748  b = b + y(1:n, n0, ii)*c18_11 - y(1:n, n1, ii)*c9_11 &
749  & + y(1:n, n2, ii)*c2_11
750 
751  b(1:p) = b(1:p) / (1._rp - h6*ay(1:p, n0_f, ii))
752 
753  y(1:n, nl, ii) = b
754 
755  end do
756  !$OMP END DO
757  !$OMP END PARALLEL
758 
759  end associate
760 
761  end subroutine nl_ms_bdfsbdf3
762 
763 
767  subroutine nl_ms_bdfsbdf4(sol, dt, N, P)
768  type(ode_solution), intent(inout) :: sol
769  real(RP) , intent(in) :: dt
770  integer , intent(in) :: N, P
771 
772  real(RP), dimension(N) :: b
773  real(RP), dimension(P) :: a
774  real(RP) :: h12, h48, h72
775  integer :: ii
776  integer :: n0, n1, n2, n3, nl, n0_F, n1_F, n2_F, n3_F
777 
778  h12 = dt * c12_25
779  h48 = dt * c48_25
780  h72 = dt * c72_25
781 
782  n0 = sol%Y_i(1)
783  n1 = sol%Y_i(2)
784  n2 = sol%Y_i(3)
785  n3 = sol%Y_i(4)
786  nl = sol%Y_i(sol%nY)
787  n0_f = sol%F_i(1)
788  n1_f = sol%F_i(2)
789  n2_f = sol%F_i(3)
790  n3_f = sol%F_i(4)
791 
792  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
793 
794  !$OMP PARALLEL PRIVATE(a, b)
795  !$OMP DO
796  do ii=1, sol%dof
797 
798  b = by(1:n, n0_f, ii)*h48 - by(1:n, n1_f, ii)*h72 &
799  & + by(1:n, n2_f, ii)*h48 - by(1:n, n3_f, ii)*h12
800 
801  a = ay(1:p, n3_f, ii) - ay(1:p, n0_f, ii)
802  a = a * y(1:p, n3_f, ii)
803  b(1:p) = b(1:p) - a * h12
804 
805  a = ay(1:p, n2_f, ii) - ay(1:p, n0_f, ii)
806  a = a * y(1:p, n2_f, ii)
807  b(1:p) = b(1:p) + a * h48
808 
809  a = ay(1:p, n1_f, ii) - ay(1:p, n0_f, ii)
810  a = a * y(1:p, n1_f, ii)
811  b(1:p) = b(1:p) - a * h72
812 
813  b = b + y(1:n, n0, ii)*c48_25 - y(1:n, n1, ii)*c36_25 &
814  & + y(1:n, n2, ii)*c16_25 - y(1:n, n3, ii)*c3_25
815 
816  b(1:p) = b(1:p) / (1._rp - h12*ay(1:p, n0_f, ii))
817 
818  y(1:n, nl, ii) = b
819 
820  end do
821  !$OMP END DO
822  !$OMP END PARALLEL
823 
824  end associate
825 
826  end subroutine nl_ms_bdfsbdf4
827 
828 
829 
833  subroutine nl_ms_bdfsbdf5(sol, dt, N, P)
834  type(ode_solution), intent(inout) :: sol
835  real(RP) , intent(in) :: dt
836  integer , intent(in) :: N, P
837 
838  real(RP), dimension(N) :: b
839  real(RP), dimension(P) :: a
840  real(RP) :: h60, h300, h600
841  integer :: ii
842  integer :: n0, n1, n2, n3, n4, nl, n0_F, n1_F, n2_F, n3_F, n4_F
843 
844  h60 = dt * c60_137
845  h300 = dt * c300_137
846  h600 = dt * c600_137
847 
848  n0 = sol%Y_i(1)
849  n1 = sol%Y_i(2)
850  n2 = sol%Y_i(3)
851  n3 = sol%Y_i(4)
852  n4 = sol%Y_i(5)
853  nl = sol%Y_i(sol%nY)
854  n0_f = sol%F_i(1)
855  n1_f = sol%F_i(2)
856  n2_f = sol%F_i(3)
857  n3_f = sol%F_i(4)
858  n4_f = sol%F_i(5)
859 
860  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
861 
862  !$OMP PARALLEL PRIVATE(a, b)
863  !$OMP DO
864  do ii=1, sol%dof
865 
866  b = by(1:n, n0_f, ii)*h300 - by(1:n, n1_f, ii)*h600 &
867  & + by(1:n, n2_f, ii)*h600 - by(1:n, n3_f, ii)*h300 &
868  & + by(1:n, n4_f, ii)*h60
869 
870  a = ay(1:p, n4_f, ii) - ay(1:p, n0_f, ii)
871  a = a * y(1:p, n4_f, ii)
872  b(1:p) = b(1:p) + a * h60
873 
874  a = ay(1:p, n3_f, ii) - ay(1:p, n0_f, ii)
875  a = a * y(1:p, n3_f, ii)
876  b(1:p) = b(1:p) - a * h300
877 
878  a = ay(1:p, n2_f, ii) - ay(1:p, n0_f, ii)
879  a = a * y(1:p, n2_f, ii)
880  b(1:p) = b(1:p) + a * h600
881 
882  a = ay(1:p, n1_f, ii) - ay(1:p, n0_f, ii)
883  a = a * y(1:p, n1_f, ii)
884  b(1:p) = b(1:p) - a * h600
885 
886  b = b + y(1:n, n0, ii)*c300_137 - y(1:n, n1, ii)*c300_137 &
887  & + y(1:n, n2, ii)*c200_137 - y(1:n, n3, ii)*c75_137 &
888  & + y(1:n, n4, ii)*c12_137
889 
890  b(1:p) = b(1:p) / (1._rp - h60*ay(1:p, n0_f, ii))
891 
892  y(1:n, nl, ii) = b
893 
894  end do
895  !$OMP END DO
896  !$OMP END PARALLEL
897 
898  end associate
899 
900  end subroutine nl_ms_bdfsbdf5
901 
902 
903 
907  subroutine nl_ms_cnab2(sol, dt, N, P)
908  type(ode_solution), intent(inout) :: sol
909  real(RP) , intent(in) :: dt
910  integer , intent(in) :: N, P
911 
912  real(RP), dimension(N) :: b
913  real(RP), dimension(P) :: a
914  real(RP) :: h1, h3
915  integer :: ii
916  integer :: n0, n1, nl, n0_F, n1_F
917 
918  h1 = dt * 0.5_rp
919  h3 = dt * 1.5_rp
920 
921  n0 = sol%Y_i(1)
922  n1 = sol%Y_i(2)
923  nl = sol%Y_i(sol%nY)
924  n0_f = sol%F_i(1)
925  n1_f = sol%F_i(2)
926 
927  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
928 
929  !$OMP PARALLEL PRIVATE(a, b)
930  !$OMP DO
931  do ii=1, sol%dof
932 
933  ! a = a_n*y_n + ( a_n - a_{n-1})*y_{n-1}
934  a = ay(1:p, n0_f, ii) - ay(1:p, n1_f, ii)
935  a = a * y(1:p, n1, ii)
936  a = a + ay(1:p, n0_f, ii)*y(1:p, n0, ii)
937 
938  ! b := y_n + (3/2)*dt*b_n - (1/2)*dt*b_{n-1}
939  b = y(1:n, n0, ii) + by(1:n, n0, ii)*h3 &
940  & - by(1:n, n1, ii)*h1
941 
942  b(1:p) = b(1:p) + a*h1
943  b(1:p) = b(1:p) / (1._rp - h1*ay(1:p, n0_f, ii))
944 
945  y(1:n, nl, ii) = b
946 
947  end do
948  !$OMP END DO
949  !$OMP END PARALLEL
950 
951  end associate
952 
953  end subroutine nl_ms_cnab2
954 
955 
959  subroutine nl_ms_mcnab2(sol, dt, N, P)
960  type(ode_solution), intent(inout) :: sol
961  real(RP) , intent(in) :: dt
962  integer , intent(in) :: N, P
963 
964  real(RP), dimension(N) :: b
965  real(RP), dimension(P) :: a
966  real(RP) :: h1, h3, CS
967  integer :: ii
968  integer :: n0, n1, nl, n0_F, n1_F
969 
970  cs = dt * c9_16
971  h1 = dt * 0.5_rp
972  h3 = dt * 1.5_rp
973 
974  n0 = sol%Y_i(1)
975  n1 = sol%Y_i(2)
976  nl = sol%Y_i(sol%nY)
977  n0_f = sol%F_i(1)
978  n1_f = sol%F_i(2)
979 
980  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
981 
982  !$OMP PARALLEL PRIVATE(a, b)
983  !$OMP DO
984  do ii=1, sol%dof
985 
986  ! a = (3/4)*an*yn + (1/8)*an*y_{n-1} + ( an-a_{n-1})*y_{n-1}
987  a = ay(1:p, n0_f, ii) - ay(1:p, n1_f, ii)
988  a = a * y(1:p, n1, ii)
989  a = a + ay(1:p, n0_f, ii)*y(1:p, n0, ii)*0.75_rp
990  a = a + ay(1:p, n0_f, ii)*y(1:p, n1, ii)*0.125_rp
991 
992 
993 
994  ! b := y_n + (3/2)*dt*b_n - (1/2)*dt*b_{n-1}
995  b = y(1:n, n0, ii) + by(1:n, n0, ii)*h3 &
996  & - by(1:n, n1, ii)*h1
997 
998  b(1:p) = b(1:p) + a*h1
999  b(1:p) = b(1:p) / (1._rp - cs*ay(1:p, n0_f, ii))
1000 
1001  y(1:n, nl, ii) = b
1002 
1003  end do
1004  !$OMP END DO
1005  !$OMP END PARALLEL
1006 
1007  end associate
1008 
1009  end subroutine nl_ms_mcnab2
1010 
1011 
1012 
1016  subroutine nl_ms_eab2(sol, dt, N, P)
1017  type(ode_solution), intent(inout) :: sol
1018  real(RP) , intent(in) :: dt
1019  integer , intent(in) :: N, P
1020 
1021  real(RP), dimension(N) :: b, w1
1022  real(RP), dimension(P) :: a
1023  integer :: ii
1024  integer :: n0, n1, nl, n0_F, n1_F
1025 
1026  n0 = sol%Y_i(1)
1027  n1 = sol%Y_i(2)
1028  nl = sol%Y_i(sol%nY)
1029  n0_f = sol%F_i(1)
1030  n1_f = sol%F_i(2)
1031 
1032  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
1033 
1034  !$OMP PARALLEL PRIVATE(a, b, w1)
1035  !$OMP DO
1036  do ii=1, sol%dof
1037 
1038  !! w1 = bn + an*yn
1039  w1 = by(1:n, n0_f, ii)
1040  w1(1:p) = w1(1:p) + ay(1:p, n0_f, ii)*y(1:p, n0, ii)
1041 
1042  ! b = gamma_1 = bn - b_{n-1} + y_{n-1}*(an - a_{n-1})
1043  b = by(1:n, n0_f, ii) - by(1:n, n1_f, ii)
1044  b(1:p) = b(1:p) + y(1:p, n1, ii) * &
1045  & ( ay(1:p, n0_f, ii) - ay(1:p, n1_f, ii) )
1046 
1047  ! a = dt*an
1048  a = dt*ay(1:p, n0_f, ii)
1049 
1050  ! b = gamma_1 + dt*an * w1
1051  b(1:p) = b(1:p) + a*w1(1:p)
1052 
1053  ! a = phi_2(an*dt)
1054  call phi_2(a)
1055 
1056  ! b = phi_2(an*dt) * (gamma_1 + dt*an*w1)
1057  b(1:p) = b(1:p) * a
1058  b(p+1:n) = b(p+1:n) * 0.5_rp
1059 
1060  ! b = w1 + phi_2(an*dt)*(gamma_1 + dt*an*w1)
1061  b = b + w1
1062 
1063  ! y_{n+1} = yn + dt * b
1064  y(1:n, nl, ii) = y(1:n, n0, ii) + dt*b
1065 
1066  end do
1067  !$OMP END DO
1068  !$OMP END PARALLEL
1069 
1070  end associate
1071 
1072  end subroutine nl_ms_eab2
1073 
1074 
1075 
1092  subroutine nl_ms_eab3(sol, dt, N, P)
1093  type(ode_solution), intent(inout) :: sol
1094  real(RP) , intent(in) :: dt
1095  integer , intent(in) :: N, P
1096 
1097  real(RP), dimension(N) :: w1, w2, b, g2, g3
1098  real(RP), dimension(P) :: a
1099  integer :: ii
1100  integer :: n0, n1, n2, nl, n0_F, n1_F, n2_F
1101 
1102  n0 = sol%Y_i(1)
1103  n1 = sol%Y_i(2)
1104  n2 = sol%Y_i(3)
1105  nl = sol%Y_i(sol%nY)
1106  n0_f = sol%F_i(1)
1107  n1_f = sol%F_i(2)
1108  n2_f = sol%F_i(3)
1109 
1110  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
1111 
1112  !$OMP PARALLEL PRIVATE(a, b, w1, w2, g2, g3)
1113  !$OMP DO
1114  do ii=1, sol%dof
1115 
1116  !! w1 = bn
1117  w1 = by(1:n, n0_f, ii)
1118 
1119  !! a = an
1120  a = ay(1:p, n0_f, ii)
1121 
1122  !! w2 = g_1 = b_{n-1} + (a_{n-1} - an) y_{n-1}
1123  w2 = by(1:n, n1_f, ii)
1124  w2(1:p) = w2(1:p) + y(1:p, n1, ii) * &
1125  & ( ay(1:p, n1_f, ii) - a)
1126 
1127  !! b = g_2 = b_{n-2} + (a_{n-2} - an) y_{n-2}
1128  b = by(1:n, n2_f, ii)
1129  b(1:p) = b(1:p) + y(1:p, n2, ii) * &
1130  & ( ay(1:p, n2_f, ii) - a)
1131 
1132  !! g2 = gamma_2 = (3/2) bn - 2 g_1 + (1/2) g_2
1133  !! g3 = gamma_3 = bn - 2 g_1 + g_2
1134  g2 = 1.5_rp*w1 -2._rp*w2 + 0.5_rp*b
1135  g3 = w1 -2._rp*w2 + b
1136 
1137  !! w1 = bn + an*yn
1138  w1(1:p) = w1(1:p) + a*y(1:p, n0, ii)
1139 
1140  !! a = an*dt
1141  a = a*dt
1142 
1143  !! w2 = gamma_2 + an h w1
1144  w2 = g2
1145  w2(1:p) = w2(1:p) + a*w1(1:p)
1146 
1147  !! b = w3 = gamma_3 + an h w2
1148  b = g3
1149  b(1:p) = b(1:p) + a*w2(1:p)
1150 
1151  !! a = phi_3(an*dt)
1152  call phi_3(a)
1153 
1154  ! b = phi_3(an*dt) * w3
1155  b(1:p) = b(1:p) * a
1156  b(p+1:n) = b(p+1:n) * c1_6
1157 
1158  !! b = w1 + w2/2 + phi_3(an*dt) * w3
1159  b = w1 + w2*0.5_rp + b
1160 
1161  ! y_{n+1} = yn + dt * b
1162  y(1:n, nl, ii) = y(1:n, n0, ii) + dt*b
1163 
1164  end do
1165  !$OMP END DO
1166  !$OMP END PARALLEL
1167 
1168  end associate
1169 
1170  end subroutine nl_ms_eab3
1171 
1190  subroutine nl_ms_eab4(sol, dt, N, P)
1191  type(ode_solution), intent(inout) :: sol
1192  real(RP) , intent(in) :: dt
1193  integer , intent(in) :: N, P
1194 
1195  real(RP), dimension(N) :: w1, w2, w3, b, g2, g3, g4
1196  real(RP), dimension(P) :: a
1197  integer :: ii
1198  integer :: n0, n1, n2, n3, nl, n0_F, n1_F, n2_F, n3_F
1199 
1200  n0 = sol%Y_i(1)
1201  n1 = sol%Y_i(2)
1202  n2 = sol%Y_i(3)
1203  n3 = sol%Y_i(4)
1204  nl = sol%Y_i(sol%nY)
1205  n0_f = sol%F_i(1)
1206  n1_f = sol%F_i(2)
1207  n2_f = sol%F_i(3)
1208  n3_f = sol%F_i(4)
1209 
1210  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY )
1211 
1212  !$OMP PARALLEL PRIVATE(a, b, w1, w2, w3, g2, g3, g4)
1213  !$OMP DO
1214  do ii=1, sol%dof
1215 
1216  !! w1 = bn
1217  w1 = by(1:n, n0_f, ii)
1218 
1219  !! a = an
1220  a = ay(1:p, n0_f, ii)
1221 
1222  !! w2 = g_1 = b_{n-1} + (a_{n-1} - an) y_{n-1}
1223  w2 = by(1:n, n1_f, ii)
1224  w2(1:p) = w2(1:p) + y(1:p, n1, ii) * &
1225  & ( ay(1:p, n1_f, ii) - a)
1226 
1227  !! w3 = g_2 = b_{n-2} + (a_{n-2} - an) y_{n-2}
1228  w3 = by(1:n, n2_f, ii)
1229  w3(1:p) = w3(1:p) + y(1:p, n2, ii) * &
1230  & ( ay(1:p, n2_f, ii) - a)
1231 
1232  !! b = g_3 = b_{n-3} + (a_{n-3} - an) y_{n-3}
1233  b = by(1:n, n3_f, ii)
1234  b(1:p) = b(1:p) + y(1:p, n3, ii) * &
1235  & ( ay(1:p, n3_f, ii) - a)
1236 
1237  ! w2 = gamma_2 = (11/6) bn - 3 g_1 + (3/2) g_2 - (1/3) g_3
1238  ! w2 = gamma_3 = 2 bn - 5 g_1 + 4 g_2 - g_3
1239  ! w2 = gamma_4 = 9 bn - 3 g_1 + 3 g_2 - g_3
1240  g2 = c11_6*w1 - 3._rp*w2 + 1.5_rp*w3 - c1_3*b
1241  g3 = 2._rp*w1 - 5._rp*w2 + 4.0_rp*w3 - b
1242  g4 = w1 - 3._rp*w2 + 3.0_rp*w3 - b
1243 
1244  !! w1 = bn + an*yn
1245  w1(1:p) = w1(1:p) + a*y(1:p, n0, ii)
1246 
1247  !! a = an*dt
1248  a = a*dt
1249 
1250  !! w2 = gamma_2 + an h w1
1251  w2 = g2
1252  w2(1:p) = w2(1:p) + a*w1(1:p)
1253 
1254  !! w3 = gamma_3 + an h w1
1255  w3 = g3
1256  w3(1:p) = w3(1:p) + a*w2(1:p)
1257 
1258  !! b = w4 = gamma_4 + an h w2
1259  b = g4
1260  b(1:p) = b(1:p) + a*w3(1:p)
1261 
1262  !! a = phi_4(an*dt)
1263  call phi_4(a)
1264 
1265  ! b = phi_4(an*dt) * w4
1266  b(1:p) = b(1:p) * a
1267  b(p+1:n) = b(p+1:n) * c1_24
1268 
1269  !! b = w1 + w2/2 + w3/6 + phi_4(an*dt) * w4
1270  b = w1 + w2*0.5_rp + w3*c1_6 + b
1271 
1272  ! y_{n+1} = yn + dt * b
1273  y(1:n, nl, ii) = y(1:n, n0, ii) + dt*b
1274 
1275  end do
1276  !$OMP END DO
1277  !$OMP END PARALLEL
1278 
1279  end associate
1280 
1281  end subroutine nl_ms_eab4
1282 
1283 end module ode_nl_ms_mod
ode_nl_ms_mod::c16_12
real(rp), parameter c16_12
Definition: ode_NL_ms_mod.f90:44
ode_nl_ms_mod::c59_24
real(rp), parameter c59_24
Definition: ode_NL_ms_mod.f90:53
ode_solution_mod::ode_output_proc
ode output management
Definition: ode_solution_mod.f90:177
ode_nl_ms_mod::c23_12
real(rp), parameter c23_12
Definition: ode_NL_ms_mod.f90:45
ode_solution_mod::ode_solution
Type ode_solution: data structure to solve ODE/PDE problems.
Definition: ode_solution_mod.f90:86
ode_nl_ms_mod::c48_25
real(rp), parameter c48_25
Definition: ode_NL_ms_mod.f90:59
ode_nl_ms_mod::c12_137
real(rp), parameter c12_137
Definition: ode_NL_ms_mod.f90:62
ode_nl_ms_mod::phi_3
elemental subroutine phi_3(z)
phi_3(z) = (exp(z) - z^2/2 - z - 1)/z^3
Definition: ode_NL_ms_mod.f90:356
ode_nl_ms_mod::c6_11
real(rp), parameter c6_11
Definition: ode_NL_ms_mod.f90:39
basic_tools
BASIC TOOLS
Definition: basic_tools.F90:8
choral_constants::ode_mcnab2
integer, parameter ode_mcnab2
ModifiedCrank Nicholson / Adamns Bashforth 2.
Definition: choral_constants.f90:206
ode_nl_ms_mod::nl_ms_bdfsbdf4
subroutine nl_ms_bdfsbdf4(sol, dt, N, P)
Lines 1 ..P = combined BDF4/SBDF4 Lines P+1..N = SBDF4.
Definition: ode_NL_ms_mod.f90:768
ode_nl_ms_mod::c2_3
real(rp), parameter c2_3
Definition: ode_NL_ms_mod.f90:32
choral_constants::ode_bdfsbdf3
integer, parameter ode_bdfsbdf3
BDF / SBDF 3.
Definition: choral_constants.f90:200
ode_nl_ms_mod::c600_137
real(rp), parameter c600_137
Definition: ode_NL_ms_mod.f90:67
ode_nl_ms_mod::c300_137
real(rp), parameter c300_137
Definition: ode_NL_ms_mod.f90:66
ode_nl_ms_mod::c1_3
real(rp), parameter c1_3
Definition: ode_NL_ms_mod.f90:31
ode_problem_mod
DERIVED TYPE ode_problem: definition of ODE/PDE problems
Definition: ode_problem_mod.f90:95
ode_nl_ms_mod::c18_11
real(rp), parameter c18_11
Definition: ode_NL_ms_mod.f90:41
ode_nl_ms_mod::nl_ms_rl2
subroutine nl_ms_rl2(sol, dt, N, P)
Lines 1 ..P = Rush-Larsen 2 Lines P+1..N = AB2.
Definition: ode_NL_ms_mod.f90:480
ode_nl_ms_mod::nl_ms_mcnab2
subroutine nl_ms_mcnab2(sol, dt, N, P)
Lines 1 ..P = combined Modified-CN/AB2 Lines P+1..N = AB2.
Definition: ode_NL_ms_mod.f90:960
basic_tools::quit
subroutine, public quit(message)
Stop program execution, display an error messahe.
Definition: basic_tools.F90:106
ode_nl_ms_mod::c55_24
real(rp), parameter c55_24
Definition: ode_NL_ms_mod.f90:52
ode_nl_ms_mod::check_ode_method_nl_ms
logical function, public check_ode_method_nl_ms(method)
is 'method' a multi-step non-linear ODE solver ?
Definition: ode_NL_ms_mod.f90:73
ode_nl_ms_mod::c4_3
real(rp), parameter c4_3
Definition: ode_NL_ms_mod.f90:33
ode_nl_ms_mod::c9_24
real(rp), parameter c9_24
Definition: ode_NL_ms_mod.f90:50
ode_nl_ms_mod::c72_25
real(rp), parameter c72_25
Definition: ode_NL_ms_mod.f90:60
choral_constants::ode_eab4
integer, parameter ode_eab4
Exponential Adamns Bashforth 4.
Definition: choral_constants.f90:198
ode_nl_ms_mod::c75_137
real(rp), parameter c75_137
Definition: ode_NL_ms_mod.f90:64
ode_def::overflow
logical function, public overflow(yy)
Detects overflow.
Definition: ode_def.f90:254
ode_nl_ms_mod::nl_ms_eab3
subroutine nl_ms_eab3(sol, dt, N, P)
Lines 1 ..P = Exponential AB3 Lines P+1..N = AB3.
Definition: ode_NL_ms_mod.f90:1093
ode_def::name_ode_method
character(len=15) function, public name_ode_method(method)
Get ODE method name.
Definition: ode_def.f90:68
ode_nl_ms_mod::c37_24
real(rp), parameter c37_24
Definition: ode_NL_ms_mod.f90:51
real_type::re
conversion integers or rational to real
Definition: real_type.F90:153
real_type
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
choral_constants::ode_eab3
integer, parameter ode_eab3
Exponential Adamns Bashforth 3.
Definition: choral_constants.f90:196
ode_solution_mod::ode_solution_init_indexes
subroutine, public ode_solution_init_indexes(sol)
initialise the ode_indexes
Definition: ode_solution_mod.f90:316
ode_nl_ms_mod::c1_24
real(rp), parameter c1_24
Definition: ode_NL_ms_mod.f90:49
choral_constants::ode_fbe
integer, parameter ode_fbe
Forward / backward Euler.
Definition: choral_constants.f90:186
ode_nl_ms_mod::c60_137
real(rp), parameter c60_137
Definition: ode_NL_ms_mod.f90:63
ode_solution_mod
DERIVED TYPE ode_solution: data straucture to solve ODEs
Definition: ode_solution_mod.f90:57
ode_nl_ms_mod::phi_2
elemental subroutine phi_2(z)
phi_2 = (exp(z) - z - 1._RP)/z^2
Definition: ode_NL_ms_mod.f90:336
ode_nl_ms_mod::phi_4
elemental subroutine phi_4(z)
phi_4(z) = (exp(z) - z^3/6 - z^2/2 - z - 1)/z^4
Definition: ode_NL_ms_mod.f90:376
ode_nl_ms_mod::c1_6
real(rp), parameter c1_6
Definition: ode_NL_ms_mod.f90:35
choral_constants
CHORAL CONSTANTS
Definition: choral_constants.f90:13
ode_nl_ms_mod::memsize_ode_nl_ms
subroutine, public memsize_ode_nl_ms(n_Y, n_FY, method)
required sizes to allocate memory
Definition: ode_NL_ms_mod.f90:98
ode_nl_ms_mod::c2_11
real(rp), parameter c2_11
Definition: ode_NL_ms_mod.f90:38
algebra_set
ALGEBRAIC OPERATIONS ON SETS
Definition: algebra_set.F90:14
choral_constants::ode_erk1
integer, parameter ode_erk1
Exponential Euler.
Definition: choral_constants.f90:189
ode_nl_ms_mod::c5_12
real(rp), parameter c5_12
Definition: ode_NL_ms_mod.f90:43
ode_nl_ms_mod::nl_ms_rl4
subroutine nl_ms_rl4(sol, dt, N, P)
Lines 1 ..P = Rush-Larsen 4 Lines P+1..N = AB4.
Definition: ode_NL_ms_mod.f90:590
ode_nl_ms_mod::nl_ms_erk1
subroutine nl_ms_erk1(sol, dt, N, P)
Lines 1 ..P = Rush-Larsen 1 Lines P+1..N = Forward Euler.
Definition: ode_NL_ms_mod.f90:431
ode_nl_ms_mod::c36_25
real(rp), parameter c36_25
Definition: ode_NL_ms_mod.f90:58
ode_nl_ms_mod::solve_ode_nl_ms
subroutine, public solve_ode_nl_ms(sol, pb, t0, T, dt, method, out, check_overflow)
solve with constant time-step
Definition: ode_NL_ms_mod.f90:210
ode_nl_ms_mod::nl_ms_bdfsbdf5
subroutine nl_ms_bdfsbdf5(sol, dt, N, P)
Lines 1 ..P = combined BDF5/SBDF5 Lines P+1..N = SBDF5.
Definition: ode_NL_ms_mod.f90:834
ode_nl_ms_mod::nl_ms_cnab2
subroutine nl_ms_cnab2(sol, dt, N, P)
Lines 1 ..P = combined CN/AB2 Lines P+1..N = AB2.
Definition: ode_NL_ms_mod.f90:908
ode_nl_ms_mod::nl_ms_bdfsbdf2
subroutine nl_ms_bdfsbdf2(sol, dt, N, P)
Lines 1 ..P = combined BDF2/SBDF2 Lines P+1..N = SBDF2.
Definition: ode_NL_ms_mod.f90:658
choral_constants::ode_bdfsbdf5
integer, parameter ode_bdfsbdf5
BDF / SBDF 5.
Definition: choral_constants.f90:202
ode_def
BOTTOM LEVEL MODULE FOR ODEs
Definition: ode_def.f90:12
ode_nl_ms_mod::nl_ms_fbe
subroutine nl_ms_fbe(sol, dt, N, P)
Lines 1 ..P = combined FORWARD/BacWard Euler Lines P+1..N = FORWARD Euler.
Definition: ode_NL_ms_mod.f90:398
choral_variables::choral_verb
integer choral_verb
Verbosity level.
Definition: choral_variables.f90:22
ode_nl_ms_mod::c200_137
real(rp), parameter c200_137
Definition: ode_NL_ms_mod.f90:65
ode_nl_ms_mod::nl_ms_bdfsbdf3
subroutine nl_ms_bdfsbdf3(sol, dt, N, P)
Lines 1 ..P = combined BDF3/SBDF3 Lines P+1..N = SBDF3.
Definition: ode_NL_ms_mod.f90:710
ode_nl_ms_mod::c9_11
real(rp), parameter c9_11
Definition: ode_NL_ms_mod.f90:40
choral_constants::ode_rl3
integer, parameter ode_rl3
Rush Larsen 3.
Definition: choral_constants.f90:191
real_type::real_tol
real(rp), parameter, public real_tol
Definition: real_type.F90:91
choral_constants::ode_bdfsbdf4
integer, parameter ode_bdfsbdf4
BDF / SBDF 4.
Definition: choral_constants.f90:201
algebra_set::circperm
subroutine, public circperm(E)
Circular permutation of an array of integer.
Definition: algebra_set.F90:271
choral_variables
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
Definition: choral_variables.f90:13
ode_nl_ms_mod::set_solver_ode_nl_ms
subroutine, public set_solver_ode_nl_ms(slv, method)
set the resolution solver
Definition: ode_NL_ms_mod.f90:145
ode_nl_ms_mod::phi_1
elemental subroutine phi_1(z)
phi_1 = ( exp(z) - 1 ) / z
Definition: ode_NL_ms_mod.f90:322
ode_nl_ms_mod::nl_ms_rl3
subroutine nl_ms_rl3(sol, dt, N, P)
Lines 1 ..P = Rush-Larsen 3 Lines P+1..N = AB3.
Definition: ode_NL_ms_mod.f90:530
ode_nl_ms_mod
MULTISTEP SOLVERS FOR NON LINEAR ODEs
Definition: ode_NL_ms_mod.f90:10
ode_nl_ms_mod::nl_ms_eab2
subroutine nl_ms_eab2(sol, dt, N, P)
Lines 1 ..P = Exponential AB2 Lines P+1..N = AB2.
Definition: ode_NL_ms_mod.f90:1017
choral_constants::ode_cnab2
integer, parameter ode_cnab2
Crank Nicholson / Adamns Bashforth 2.
Definition: choral_constants.f90:204
ode_nl_ms_mod::c11_6
real(rp), parameter c11_6
Definition: ode_NL_ms_mod.f90:36
choral_constants::ode_bdfsbdf2
integer, parameter ode_bdfsbdf2
BDF / SBDF 2.
Definition: choral_constants.f90:199
ode_nl_ms_mod::nl_ms_eab4
subroutine nl_ms_eab4(sol, dt, N, P)
Lines 1 ..P = Exponential AB4 Lines P+1..N = AB4.
Definition: ode_NL_ms_mod.f90:1191
choral_constants::ode_eab2
integer, parameter ode_eab2
Exponential Adamns Bashforth 2.
Definition: choral_constants.f90:194
ode_problem_mod::ode_problem
Type ode_problem: definition of ODE/PDE problems.
Definition: ode_problem_mod.f90:122
ode_nl_ms_mod::c9_16
real(rp), parameter c9_16
Definition: ode_NL_ms_mod.f90:47
ode_nl_ms_mod::c16_25
real(rp), parameter c16_25
Definition: ode_NL_ms_mod.f90:57
choral_constants::ode_rl2
integer, parameter ode_rl2
Rush Larsen 2.
Definition: choral_constants.f90:190
ode_nl_ms_mod::c12_25
real(rp), parameter c12_25
Definition: ode_NL_ms_mod.f90:56
ode_nl_ms_mod::c3_25
real(rp), parameter c3_25
Definition: ode_NL_ms_mod.f90:55
ode_nl_ms_mod::create_ode_nl_ms_sol
subroutine, public create_ode_nl_ms_sol(sol, pb, method)
create memory for the ode_solution structure 'sol'
Definition: ode_NL_ms_mod.f90:188
choral_constants::ode_rl4
integer, parameter ode_rl4
Rush Larsen 4.
Definition: choral_constants.f90:192
Generated on Mon May 17 2021 15:14:13 for Choral by   doxygen 1.8.13