Choral
ode_opSplt_mod.f90
Go to the documentation of this file.
1 
28 
29 
31 
34  use real_type
35  use basic_tools
36  use r1d_mod, only: copy
37  use krylov_mod
38  use r1d_mod, only: xpay
39  use ode_def
40  use ode_problem_mod
42  use ode_nl_1s_mod
43  use ode_lin_1s_mod
44 
45  implicit none
46  private
47 
48  ! %----------------------------------------%
49  ! | |
50  ! | PUBLIC DATA |
51  ! | |
52  ! %----------------------------------------%
53  !! DERIVED TYPE
54  public :: ode_opsplt
55 
56  !! GENERIC SUBROUTINES
57  public :: clear
58  public :: valid
59  public :: print
60 
61  !! NON GENERIC SUBROUTINES
62  public :: solve_ode_opsplt
63  public :: name_ode_method_opsplt
64  public :: order_ode_method_opsplt
65  public :: create_ode_opsplt_sol
66 
67 
68 
69  ! %----------------------------------------%
70  ! | |
71  ! | DERIVED TYPE |
72  ! | |
73  ! %----------------------------------------%
77  type :: ode_opsplt
78 
80  integer :: method = -1
81 
83  integer :: ns = -1
84 
86  integer :: l_meth = -1
87 
89  integer :: nl_meth = -1
90 
92  real(RP), dimension(:), allocatable :: a_nl
93 
95  real(RP), dimension(:), allocatable :: a_l
96 
99  logical , dimension(:), allocatable :: nl_first
100 
102  procedure(ode_lin_solver ), nopass, pointer :: slv_l =>null()
103 
105  procedure(ode_nl_1s_solver), nopass, pointer :: slv_nl=>null()
106 
108  logical :: check_overflow = .false.
109 
110  contains
111 
114 
115  end type ode_opsplt
116 
117 
118  ! %----------------------------------------%
119  ! | |
120  ! | GENERIc SUBROUTINES |
121  ! | |
122  ! %----------------------------------------%
123 
125  interface clear
126  module procedure ode_opsplt_clear
127  end interface clear
128 
129  interface ode_opsplt
130  module procedure ode_opsplt_create
131  end interface ode_opsplt
132 
134  interface valid
135  module procedure ode_opsplt_valid
136  end interface valid
137 
139  interface print
140  module procedure ode_opsplt_print
141  end interface print
142 
143 
144 contains
145 
148  subroutine ode_opsplt_clear(os)
149  type(ode_opsplt), intent(inout):: os
150 
151  os%method = -1
152  os%ns = -1
153  os%L_meth = -1
154  os%NL_meth = -1
155  os%check_overflow = .false.
156 
157  call freemem(os%a_NL)
158  call freemem(os%a_L)
159  call freemem(os%NL_first)
160 
161  os%slv_L =>null()
162  os%slv_NL=>null()
163 
164 
165  end subroutine ode_opsplt_clear
166 
167 
181  function ode_opsplt_create(method, L_meth, NL_meth, &
182  & check_overflow) result(os)
183  type(ode_opsplt) :: os
184  integer , intent(in) :: method, L_meth, NL_meth
185  logical, optional, intent(in) :: check_overflow
186 
187  logical :: bool
188 
189  if (choral_verb>0) write(*,*) &
190  &"ode_opSplt_mod : ode_opSplt_create"
191  call clear(os)
192 
193  os%method = method
194  os%L_meth = l_meth
195  os%NL_meth = nl_meth
196 
197  bool = valid(os)
198  if (.NOT.bool) call quit("ode_opSplt_mod: ode_opSplt_set")
199 
200  call set_solver_ode_lin_1s(os%slv_L, l_meth)
201  call set_solver_ode_nl_1s(os%slv_NL, nl_meth)
202  call def_splitting(os)
203 
204  !! check overflow when solving
205  if (present(check_overflow)) then
206  os%check_overflow = check_overflow
207  end if
208 
209  if (.NOT.valid(os)) call quit(&
210  & "ode_opSplt_mod: ode_opSplt_create: not valid")
211 
212  end function ode_opsplt_create
213 
214 
217  subroutine def_splitting(os)
218  type(ode_opsplt), intent(inout) :: os
219 
220  real(RP) :: theta
221 
222  !! number of stages
223  !!
224  select case(os%method)
225  case(ode_os_lie)
226  os%ns = 1
227  case(ode_os_strang)
228  os%ns = 2
229  case(ode_os_ruth, ode_os_aks3)
230  os%ns = 3
231  case(ode_os_yoshida)
232  os%ns = 4
233  case default
234  call quit( 'rdEq_mod: def_splitting: 1' )
235  end select
236 
237  call allocmem(os%a_NL, os%ns)
238  call allocmem(os%a_L , os%ns)
239  call allocmem(os%NL_first, os%ns)
240 
241  !! Definition of a_L(:), a_NL(:) and NL_first(:)
242  !!
243  os%NL_first = .true. ! default = reaction first
244  select case(os%method)
245  case(ode_os_lie)
246  os%a_NL=(/1._rp/)
247  os%a_L=(/1._rp/)
248 
249  case(ode_os_strang)
250  os%a_NL=(/.5_rp, .5_rp/)
251  os%a_L=(/.5_rp, .5_rp/)
252  os%NL_first(2) = .false.
253 
254  case(ode_os_ruth)
255  os%a_L=(/ re(2,3) ,-re(2,3), 1._rp/)
256  os%a_NL=(/ re(7,24), re(3,4),-re(1,24)/)
257  os%NL_first = .true.
258 
259  case(ode_os_aks3)
260  os%a_L(1)= 0.91966152301739985_rp
261  os%a_L(2)=-0.18799161879915978_rp
262  os%a_L(3)= 0.26833009578175992_rp
263 
264  os%a_NL(1)=os%a_L(3)
265  os%a_NL(2)=os%a_L(2)
266  os%a_NL(3)=os%a_L(1)
267 
268  os%NL_first = .true.
269 
270  case(ode_os_yoshida)
271  theta = 1._rp/(2._rp - 2._rp ** ( 1._rp / 3._rp ) )
272 
273  os%a_NL=(/ theta/ 2._rp, (1._rp-theta)/2._rp, &
274  & (1._rp-theta)/2._rp, theta/2._rp /)
275 
276  os%a_L=(/ theta, 1._rp-2._rp*theta, theta, 0._rp/)
277 
278  end select
279 
280  end subroutine def_splitting
281 
282 
283  !! check ode_opSplt
284  !!
285  function ode_opsplt_valid(os) result(bool)
286  type(ode_opsplt), intent(in):: os
287 
288  logical :: bool
289 
290  bool = (os%method>0) .AND. (os%method <= ode_os_tot_nb)
291  bool = bool .AND. check_ode_method_lin_1s(os%L_meth)
292  bool = bool .AND. check_ode_method_nl_1s(os%NL_meth)
293 
294  end function ode_opsplt_valid
295 
296 
297  !! print ode_opSplt
298  !!
299  subroutine ode_opsplt_print(os)
300  type(ode_opsplt), intent(in):: os
301 
302  write(*,*)"ode_opSplt_mod : print"
303 
304  write(*,*)" operator spltiting Method = ",&
305  & trim( name_ode_method_opsplt(os%method) )
306 
307  write(*,*)" one step Non-linear ODE solver = ",&
308  & trim(name_ode_method(os%NL_meth))
309 
310  write(*,*)" one step linear ODE solver = ",&
311  & trim(name_ode_method(os%L_meth))
312 
313  write(*,*)" check overfow =", os%check_overflow
314 
315  end subroutine ode_opsplt_print
316 
317 
320  function name_ode_method_opsplt(method) result(name)
321  integer, intent(in) :: method
322  character(len=15) :: name
323 
324  select case(method)
325  case(ode_os_lie)
326  name="Lie OS"
327 
328  case(ode_os_strang)
329  name="Strang OS"
330 
331  case(ode_os_ruth)
332  name="Ruth OS"
333 
334  case(ode_os_aks3)
335  name="AKS3 OS"
336 
337  case(ode_os_yoshida)
338  name="Yoshida OS"
339 
340  case default
341  name="invalid"
342 
343  end select
344 
345  end function name_ode_method_opsplt
346 
347 
350  function order_ode_method_opsplt(method) result(o)
351  integer, intent(in) :: method
352  integer :: o
353 
354  select case(method)
355  case(ode_os_lie)
356  o = 1
357 
358  case(ode_os_strang)
359  o = 2
360 
361  case(ode_os_ruth, ode_os_aks3)
362  o = 3
363 
364  case(ode_os_yoshida)
365  o = 4
366 
367  case default
368  o= -1
369 
370  end select
371 
372  end function order_ode_method_opsplt
373 
374 
375 
376  subroutine create_ode_opsplt_sol(sol, pb, os)
377  type(ode_solution), intent(inout) :: sol
378  type(ode_problem) , intent(in) :: pb
379  type(ode_opsplt) , intent(in) :: os
380 
381  integer :: n_V, n_Y, n_FY
382  logical :: bool
383 
384  bool = valid(os)
385  if (.NOT.bool) call quit(&
386  & "ode_opSplt_mod: create_ode_opSplt_sol: 1")
387 
388  call memsize_ode_nl_1s(n_y, n_fy, os%NL_meth)
389  call memsize_ode_lin_1s(n_v, os%L_meth)
390 
391  sol = ode_solution(pb, nv=n_v, ny=n_y, nfy=n_fy)
392 
393  end subroutine create_ode_opsplt_sol
394 
395 
398  subroutine solve_ode_opsplt(sol, t0, T, dt, os, pb, kry, out)
399  type(ode_solution), intent(inout) :: sol
400  real(RP) , intent(in) :: dt, t0, T
401  type(ode_opsplt) , intent(in) :: os
402  type(ode_problem) , intent(in) :: pb
403  type(krylov) , intent(inout) :: kry
404  procedure(ode_output_proc) :: out
405 
406  real(RP) :: CS, tn, ts, h
407  integer :: ii, jj, nStep
408  logical :: exit_comp, ierr
409 
410  if (choral_verb>1) write(*,*) &
411  &"ode_opSplt_mod : solve"
412  if (choral_verb>2) write(*,*) &
413  &" Operator splitting method = ",&
414  & name_ode_method_opsplt(os%method)
415  if (choral_verb>2) write(*,*) &
416  &" NonLin onestep solver = ",&
417  & name_ode_method(os%NL_meth)
418  if (choral_verb>2) write(*,*) &
419  &" Linear onestep solver = ",&
420  & name_ode_method(os%L_meth)
421 
422  sol%ierr = 0
423  nstep = int( (t-t0) / dt)
424  exit_comp = .false.
425 
426  !! Time loop
427  !!
428  time_loop: do ii=1, nstep
429  tn = t0 + re(ii-1)*dt
430 
431  ! output
432  call out(tn, sol, exit_comp)
433  if (exit_comp) then
434  if (choral_verb>2) &
435  &write(*,*)"ode_opSplt_mod :&
436  & solve = ", tn, ': EXIT_COMPp'
437  return
438  end if
439 
440  !! Loop on operator splitting stages
441  !!
442  ts = tn
443  stage_loop: do jj=1, os%ns
444 
445  if (os%NL_first(jj)) then
446  !! Stage starting with the non-linear solver
447 
448  !! non-linear equation
449  h = dt * os%a_NL(jj)
450  call os%slv_NL(sol, h, ts, pb)
451  call copy(sol%V(:,1), sol%Y(sol%N, 1, :))
452  if (os%check_overflow) then
453  ierr = overflow(sol%Y(:,1,:))
454  if (ierr) then
455  sol%ierr = 10
456  exit time_loop
457  end if
458  end if
459 
460  !! linear equation
461  h = dt * os%a_L(jj)
462  cs = s_prefactor(os%L_meth, h)
463  call os%slv_L(sol, ierr, h, pb, kinv)
464  call copy(sol%Y(sol%N,1,:), sol%V(:,1))
465  if (ierr) then
466  sol%ierr = 1
467  exit time_loop
468  end if
469 
470  else
471  !! Stage starting with the linear equationN
472 
473  !! linear equation
474  h = dt * os%a_L(jj)
475  cs = s_prefactor(os%L_meth, h)
476  call os%slv_L(sol, ierr, h, pb, kinv)
477  call copy(sol%Y(sol%N,1,:), sol%V(:,1))
478 
479  if (ierr) then
480  sol%ierr = 1
481  exit time_loop
482  end if
483 
484  ! non-linear equation
485  h = dt * os%a_NL(jj)
486  call os%slv_NL(sol, h, ts, pb)
487  call copy(sol%V(:,1), sol%Y(sol%N,1,:))
488 
489  if (os%check_overflow) then
490  ierr = overflow(sol%Y(:,1,:))
491  if (ierr) then
492  sol%ierr = 10
493  exit time_loop
494  end if
495  end if
496 
497  end if
498 
499  ts = ts + os%a_NL(jj)*dt
500  end do stage_loop
501 
502  end do time_loop
503 
504  if (sol%ierr>0) then
505 
506  if (choral_verb>1) then
507  select case(sol%ierr)
508  case(10)
509  write(*,*) "ode_opSplt_mod : solve,&
510  & time =", real(tn, SP), ", stage =", int(jj,1), &
511  & ': OVERFLOW Y'
512 
513  case(1)
514  write(*,*) "ode_opSplt_mod : solve,&
515  & time =", real(tn, SP), ", stage =", int(jj,1), &
516  & ': linear pb inversion error'
517 
518  case default
519  write(*,*) "ode_opSplt_mod : solve,&
520  & time =", real(tn, SP), ", stage =", int(jj,1), &
521  & ': ERROR'
522 
523  end select
524  end if
525  end if
526 
527  contains
528 
531  subroutine kinv(x, bool, b)
532  real(RP), dimension(:), intent(inout) :: x
533  logical , intent(out) :: bool
534  real(RP), dimension(:), intent(in) :: b
535 
536  call solve(x, kry, b, k)
537  bool = kry%ierr
538 
539  end subroutine kinv
540 
541  !! matrix-vector product x --> K*x, K = M + CS*S
542  !!
543  subroutine k(y, x)
544  real(RP), dimension(:), intent(out) :: y
545  real(RP), dimension(:), intent(in) :: x
546 
547  call pb%M(y , x)
548  call pb%S(sol%aux, x)
549  call xpay(y, cs, sol%aux)
550 
551  end subroutine k
552 
553  end subroutine solve_ode_opsplt
554 
555 
556 end module ode_opsplt_mod
integer, parameter ode_os_yoshida
check &#39;ode_opSplt&#39; parameters
type(ode_opsplt) function ode_opsplt_create(method, L_meth, NL_meth, check_overflow)
Constructor for the type ode_opSplt
Type ode_solution: data structure to solve ODE/PDE problems.
BASIC TOOLS
Definition: basic_tools.F90:8
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
integer, parameter ode_os_ruth
character(len=15) function, public name_ode_method_opsplt(method)
Get operator splitting method name.
DERIVED TYPE ode_problem: definition of ODE/PDE problems
subroutine, public solve_ode_opsplt(sol, t0, T, dt, os, pb, kry, out)
solve : operator splitting with constant time step
subroutine, public create_ode_opsplt_sol(sol, pb, os)
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
The type krylov defines the settings of a linear solver.
Definition: krylov_mod.f90:63
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine, public set_solver_ode_nl_1s(slv, method)
set the solver &#39;slv&#39; to a predefined solver being given a method
integer, parameter ode_os_strang
subroutine ode_opsplt_clear(os)
destructor for &#39;ode_opSplt&#39;
The type opSplt defines operator spltting methods.
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
conversion integers or rational to real
Definition: real_type.F90:153
subroutine, public memsize_ode_lin_1s(n_V, method)
required sizes to allocate memory
subroutine kinv(x, bool, b)
Solver for K*x = b.
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
logical function, public check_ode_method_nl_1s(method)
is &#39;method&#39; a one-step non-linear ODE solver ?
subroutine ode_opsplt_print(os)
ONE-STEP SOLVERS FOR LINEAR ODEs
integer function, public order_ode_method_opsplt(method)
Get operator splitting method name.
logical function ode_opsplt_valid(os)
DERIVED TYPE ode_solution: data straucture to solve ODEs
y = x
Definition: R1d_mod.F90:44
subroutine, public set_solver_ode_lin_1s(slv, method)
Setting the solver for diffusion.
CHORAL CONSTANTS
print &#39;ode_opSplt&#39; parameters
allocate memory for real(RP) arrays
Definition: real_type.F90:158
BOTTOM LEVEL MODULE FOR ODEs
Definition: ode_def.f90:12
subroutine, public memsize_ode_nl_1s(n_Y, n_FY, method)
required sizes to allocate memory
integer choral_verb
Verbosity level.
x = x + b*y // OR // z = x + b*y
Definition: R1d_mod.F90:60
subroutine k(y, x)
ONE-STEP SOLVERS FOR NON LINEAR ODEs
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
logical function, public check_ode_method_lin_1s(method)
is &#39;method&#39; a one-step linear ODE solver ?
integer, parameter ode_os_aks3
DERIVED TYPE ode_opSplt: operator splitting methods for ODEs.
Type ode_problem: definition of ODE/PDE problems.
DERIVED TYPE krylov: for the resolution of linear systems
Definition: krylov_mod.f90:25
integer, parameter ode_os_lie
subroutine def_splitting(os)
Define splitting methods.
integer, parameter ode_os_tot_nb
Number of operator splitting methods.