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
60 integer,
intent(in) :: method
82 integer,
intent(out) :: n_V, n_FY
83 integer,
intent(in) :: method
106 call quit(
"ode_SL_ms_mod: memSize_ode_SL_ms: error")
116 integer,
intent(in) :: method
117 procedure(ode_lin_solver),
pointer :: slv
138 call quit(
"ode_SL_ms_mod:& 139 & set_solver_ode_SL_ms: uncorrect method")
150 integer ,
intent(in) :: method
156 if (.NOT.bool)
call quit(&
157 &
"ode_SL_ms_mod: create_ode_SL_ms_sol: uncorrect method")
159 if (pb%N/=1)
call quit(&
160 &
"ode_SL_ms_mod: create_ode_SL_ms_sol: wrong problem size")
172 & out, check_overflow, Kinv, kry)
175 real(RP) ,
intent(in) :: t0, T, dt
176 integer ,
intent(in) :: method
178 logical ,
intent(in) :: check_overflow
179 procedure(linsystem_solver),
optional :: Kinv
180 type(
krylov) ,
optional :: kry
182 procedure(ode_lin_solver) ,
pointer :: slv=>null()
183 procedure(linsystem_solver),
pointer :: KInv_1
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
192 &
"ode_SL_ms_mod : solve_ode_SL_ms" 194 &
" SemiLin multistep solver = ",&
197 &
" K_inv provided = ",&
199 if (.NOT.
present(kinv))
then 201 &
" Krylov settings provided = ",&
211 if (
present(kinv))
then 215 if (
present(kry)) kry_1 = kry
234 ns = int( (t-t0) / dt)
239 tn = t0 +
re(jj-1)*dt
243 call out(tn, sol, exit_comp)
246 &
"ode_SL_ms_mod : solve& 247 & time = ", tn,
': EXIT_COMPp' 250 t_out = t_out +
clock() - cpu
255 call slv(sol, ierr, dt, pb, kinv_1)
256 t_sl = t_sl +
clock() - cpu
263 &
"ode_SL_ms_mod : solve,& 264 & time = ", tn,
': linear system solver error' 271 associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY, &
290 y(1, n0, ii) = v(ii, n0_v)
292 call pb%AB(ay0, by(:, n0_f,ii), &
293 & x(:,ii), tn+dt, y(:,n0,ii), 1, 0 )
301 y(1, n0, ii) = v(ii, n0_v)
303 call pb%AB( ay(:,n0_f,ii), by(:,n0_f,ii), &
304 & x(:,ii), tn+dt, y(:,n0,ii), 1, 1 )
306 by(1,n0_f, ii) = by(1,n0_f, ii) &
307 & + ay(1,n0_f, ii)*v(ii, n0_v)
315 if (check_overflow)
then 320 &
"ode_NL_ms_mod : solve,& 321 & time = ", tn,
': OVERFLOW Y' 329 &
"ode_NL_ms_mod : solve,& 330 & time = ", tn,
': OVERFLOW BY' 339 &
"ode_NL_ms_mod : solve,& 340 & time = ", tn,
': OVERFLOW AY' 347 t_reac = t_reac +
clock() - cpu
352 write(*,*)
"ode_SL_NL_mod : solve, end = timing" 353 write(*,*)
" Semilin. eq. integration =", &
355 write(*,*)
" Reaction term evaluation =", &
357 write(*,*)
" Output =", &
359 t_tot =
clock() - t_tot
360 write(*,*)
" Total CPU =", &
370 real(RP),
dimension(:),
intent(inout) :: x
371 logical ,
intent(out) :: bool
372 real(RP),
dimension(:),
intent(in) :: b
382 real(RP),
dimension(:),
intent(out) :: y
383 real(RP),
dimension(:),
intent(in) :: x
386 call pb%S(sol%aux, x)
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
403 integer :: n0_V, nl_V, n0_F
406 nl_v = sol%V_i(sol%nV)
409 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
411 call xpay(aux, v(:, n0_v), dt, by(pb%N, n0_f, :))
415 if (n0_v /= nl_v)
then 416 call copy(v(:,nl_v), v(:,n0_v))
420 call kinv(v(:,nl_v), ierr, rhs)
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
436 integer :: n0_V, nl_V, n0_F
439 nl_v = sol%V_i(sol%nV)
442 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
444 call xpay(aux, v(:, n0_v), dt, by(pb%N, n0_f, :))
446 call pb%S(aux, v(:, n0_v))
447 call xpay(rhs, -dt, aux)
451 if (n0_v /= nl_v)
then 452 call copy(v(:,nl_v), v(:,n0_v))
456 call kinv(v(:,nl_v), ierr, rhs)
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
472 integer :: ii, N, n0_V, nl_V, n0_F, n1_F
473 real(RP) :: h3_2, h1_2
480 nl_v = sol%V_i(sol%nV)
484 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
490 aux(ii) = v(ii,n0_v) &
491 & + h3_2 * by(n,n0_f,ii) - h1_2 * by(n,n1_f,ii)
498 call pb%S(aux, v(:,n0_v))
499 call xpay(rhs, -h1_2, aux)
502 if (sol%NV>1)
call copy( v(:,nl_v), v(:,n0_v))
505 call kinv(v(:,nl_v), ierr, rhs)
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
521 integer :: ii, N, n0_V, n1_V, nl_V, n0_F, n1_F
522 real(RP) :: h3_2, h1_2
530 nl_v = sol%V_i(sol%nV)
534 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
540 aux(ii) = v(ii,n0_v) &
541 & + h3_2 * by(n,n0_f,ii) - h1_2 * by(n,n1_f,ii)
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)
554 call copy( v(:,nl_v), v(:,n0_v))
557 call kinv(v(:,nl_v), ierr, rhs)
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
591 integer :: ii, N, n0_V, n1_V, nl_V, n0_F, n1_F
592 real(RP) :: h4_3, h2_3
600 nl_v = sol%V_i(sol%nV)
604 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
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)
620 call copy( v(:,nl_v), v(:,n0_v))
623 call kinv(v(:,nl_v), ierr, rhs)
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
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
649 nl_v = sol%V_i(sol%nV)
654 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
661 & +
c2_11 * v(ii,n2_v)
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)
675 call copy( v(:,nl_v), v(:,n0_v))
678 call kinv(v(:,nl_v), ierr, rhs)
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
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
708 nl_v = sol%V_i(sol%nV)
714 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
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)
735 call copy( v(:,nl_v), v(:,n0_v))
738 call kinv(v(:,nl_v), ierr, rhs)
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
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
768 nl_v = sol%V_i(sol%nV)
775 associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)
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)
798 call copy( v(:,nl_v), v(:,n0_v))
801 call kinv(v(:,nl_v), ierr, rhs)
Type ode_solution: data structure to solve ODE/PDE problems.
integer, parameter ode_mcnab2
ModifiedCrank Nicholson / Adamns Bashforth 2.
integer, parameter ode_bdfsbdf3
BDF / SBDF 3.
subroutine k_default(y, x)
DERIVED TYPE ode_problem: definition of ODE/PDE problems
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...
real(rp), parameter c12_25
The type krylov defines the settings of a linear solver.
subroutine, public set_solver_ode_sl_ms(slv, method)
set the resolution solver
real(rp), parameter c300_137
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
subroutine sl_ms_bdfsbdf5(sol, ierr, dt, pb, KInv)
BDF5/SBDF5.
logical function, public overflow(yy)
Detects overflow.
character(len=15) function, public name_ode_method(method)
Get ODE method name.
real(rp), parameter c6_11
conversion integers or rational to real
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
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.
subroutine kinv_default(x, bool, b)
Default solver for K*x = b.
subroutine sl_ms_mcnab2(sol, ierr, dt, pb, KInv)
Modified-CN/AB2.
real(rp), parameter c48_25
real(rp), parameter c3_25
ALGEBRAIC OPERATIONS ON SETS
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.
integer, parameter ode_bdfsbdf5
BDF / SBDF 5.
BOTTOM LEVEL MODULE FOR ODEs
integer choral_verb
Verbosity level.
real(rp), parameter c12_137
x = x + b*y // OR // z = x + b*y
real(rp), parameter c2_11
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
real(rp), parameter c75_137
real(rp), parameter c600_137