Choral
ode_NL_DC_mod.f90
Go to the documentation of this file.
1 
11 
12 
14 
17  use real_type
18  use basic_tools
19  use algebra_set
20  use r1d_mod, only: xpay
21  use ode_def
22  use ode_problem_mod
24  use krylov_mod
25  !$ use OMP_LIB
26 
27  implicit none
28  private
29 
30  ! %----------------------------------------%
31  ! | |
32  ! | PUBLIC DATA |
33  ! | |
34  ! %----------------------------------------%
35 
36  public :: solve_ode_nl_dc
37  public :: create_ode_nl_dc_sol
38  public :: check_ode_method_nl_dc
39 
40 contains
41 
42 
44  function check_ode_method_nl_dc(method) result(b)
45  logical :: b
46  integer, intent(in) :: method
47 
48  b = .false.
49 
50  select case(method)
51  case(ode_dc_2, ode_dc_3)
52  b = .true.
53 
54  end select
55 
56  end function check_ode_method_nl_dc
57 
60  subroutine memsize_ode_nl_dc(n_Y, n_FY, method)
61  integer, intent(out) :: n_FY, n_Y
62  integer, intent(in) :: method
63 
64  select case(method)
65 
66  case(ode_dc_2)
67  n_y = 3; n_fy = 1
68 
69  case(ode_dc_3)
70  n_y = 6; n_fy = 3
71 
72  case default
73  call quit("ode_NL_DC_mod: memSize_ode_NL_DC: uncorrect method")
74 
75  end select
76 
77  end subroutine memsize_ode_nl_dc
78 
79 
82  subroutine create_ode_nl_dc_sol(sol, pb, method)
83  type(ode_solution), intent(inout) :: sol
84  type(ode_problem) , intent(in) :: pb
85  integer , intent(in) :: method
86 
87  integer :: n_Y, n_FY
88 
89  if (pb%Na/=0) call quit("ode_NL_DC_mod: create_ode_NL_DC:&
90  & Na muste be ezero")
91 
92  call memsize_ode_nl_dc(n_y, n_fy, method)
93  sol = ode_solution(pb, ny=n_y, nfy=n_fy)
94 
95  end subroutine create_ode_nl_dc_sol
96 
97 
98  subroutine init_nl_dc(sol, pb, meth, t0, dt)
99  type(ode_solution), intent(inout) :: sol
100  type(ode_problem) , intent(in) :: pb
101  integer , intent(in) :: meth
102  real(RP) , intent(in) :: t0, dt
103 
104  select case(meth)
105  case(ode_dc_2)
106  call init_dc_2(sol, pb, t0)
107 
108  case(ode_dc_3)
109  call init_dc_3(sol, pb, t0, dt)
110 
111  case default
112  call quit("ode_NL_DC_mod: init_NL_DC: wrong method")
113 
114  end select
115 
116  end subroutine init_nl_dc
117 
118 
119  subroutine init_dc_2(sol, pb, t0)
120  type(ode_solution), intent(inout) :: sol
121  type(ode_problem) , intent(in) :: pb
122  real(RP) , intent(in) :: t0
123 
124  real(RP), dimension(0) :: ay
125 
126  integer :: ii, N
127 
128  n = pb%N
129 
130  associate( y=>sol%Y, by=>sol%BY, x=>pb%X)
131 
132  !$OMP PARALLEL PRIVATE(ay)
133  !$OMP DO
134  do ii=1, sol%dof
135  y(:,2,ii) = y(:,1,ii)
136  y(:,3,ii) = 0.0_rp
137  call pb%AB( ay, by(:,1,ii), &
138  & x(:,ii), t0, y(:,1,ii), n, 0 )
139  end do
140  !$OMP END DO
141  !$OMP END PARALLEL
142 
143  end associate
144  end subroutine init_dc_2
145 
146 
147  subroutine init_dc_3(sol, pb, t0, dt)
148  type(ode_solution), intent(inout) :: sol
149  type(ode_problem) , intent(in) :: pb
150  real(RP) , intent(in) :: t0, dt
151 
152  integer :: ii, N
153  real(RP), dimension(0) :: ay
154 
155  n = pb%N
156 
157  associate( y=>sol%Y, fy=>sol%BY, x=>pb%X)
158 
159  !$OMP PARALLEL PRIVATE(ay)
160  !$OMP DO
161  do ii=1, sol%dof
162 
163  y(:,2,ii) = y(:,1,ii) ! Y_0^0
164  y(:,3,ii) = 0.0_rp ! Y_1^0
165  y(:,4,ii) = 0.0_rp ! Y_2^0
166 
167  ! F_0^0
168  call pb%AB( ay, fy(:,1,ii), &
169  & x(:,ii), t0, y(:,2,ii), n, 0 )
170 
171  ! Y_0^1
172  y(:,5,ii) = y(:,2,ii) + dt*fy(:,1,ii)
173 
174  ! F_0^1
175  call pb%AB( ay, fy(:,2,ii), &
176  & x(:,ii), t0+dt, y(:,5,ii), n, 0 )
177 
178  ! F_1^1
179  fy(:,3,ii) = fy(:,2,ii)
180 
181  ! Y_1^1
182  y(:,6,ii) = 0.5_rp*(fy(:,2,ii)-fy(:,1,ii))
183 
184  end do
185  !$OMP END DO
186  !$OMP END PARALLEL
187 
188  end associate
189  end subroutine init_dc_3
190 
193  subroutine solve_ode_nl_dc(sol, pb, t0, T, dt, &
194  & meth, out, check_overflow)
195  type(ode_solution), intent(inout) :: sol
196  type(ode_problem) , intent(in) :: pb
197  real(RP) , intent(in) :: t0, T, dt
198  integer , intent(in) :: meth
199  procedure(ode_output_proc) :: out
200  logical , intent(in) :: check_overflow
201 
202  real(RP) :: tn
203  integer :: nS, jj
204  logical :: exit_comp, b
205 
206  if (choral_verb>1) write(*,*) &
207  &"ode_NL_DC_mod: solve"
208  if (choral_verb>2) write(*,*) &
209  &" Method = ",&
210  & name_ode_method(meth)
211 
212  b = check_ode_method_nl_dc(meth)
213  if (.NOT.b) call quit("ode_NL_DC_mod: solve_ode_NL_DC:&
214  & wrong method")
215 
216  if (pb%Na/=0) call quit("ode_NL_DC_mod: solve_ode_NL_DC:&
217  & Na muste be ezero")
218 
219  !! Initialise
220  !!
221  call init_nl_dc(sol, pb, meth, t0, dt)
222 
223  !! ode resolution
224  !!
225  ns = int( (t-t0) / dt)
226  exit_comp = .false.
227  sol%ierr = 0
228  !!
229  !! Time loop
230  !!
231  do jj=1, ns
232  tn = t0 + re(jj-1)*dt
233 
234  ! output
235  call out(tn, sol, exit_comp)
236  if (exit_comp) then
237  if (choral_verb>2) write(*,*) &
238  &"ode_NL_DC_mod: solve&
239  & time = ", tn, ': EXIT_COMPp'
240  return
241  end if
242 
243  select case(meth)
244  case(ode_dc_2)
245  call nl_dc_2(sol, dt, tn, pb)
246 
247  case(ode_dc_3)
248  call nl_dc_3(sol, dt, tn, pb)
249 
250  end select
251 
252  end do
253 
254  end subroutine solve_ode_nl_dc
255 
256 
257 
260  subroutine nl_dc_2(sol, dt, t, pb)
261  type(ode_solution), intent(inout) :: sol
262  real(RP) , intent(in) :: dt, t
263  type(ode_problem) , intent(in) :: pb
264 
265  integer :: ii, N
266  real(RP), dimension(pb%N) :: y0, F0, y1, F1
267  real(RP), dimension(0) :: ay
268 
269  n = pb%N
270 
271  associate( y=>sol%Y, fy=>sol%BY, x=>pb%X)
272 
273  !$OMP PARALLEL PRIVATE(y0, F0, y1, F1, ay)
274  !$OMP DO
275  do ii=1, sol%dof
276 
277  !! y0 = y_0^{n+1} := y_0^{n} + dt * F(t_n, y_0^{n})
278  y0 = y(:,2,ii) + dt*fy(:,1,ii)
279 
280  !! F0 = F_0^{n+1} := F(y0^{n+1}, t_{n+1})
281  call pb%AB(ay, f0, x(:,ii), t+dt, y0, n, 0 )
282 
283  !! F1 = F_1^{n+1} := F(t_{n+1}, y_0^{n+1} + dt*y_1^{n})
284  y1 = y0 + dt*y(:,3,ii)
285  call pb%AB(ay, f1, x(:,ii), t+dt, y1, n, 0 )
286 
287  !! y1 y_1^{n+1}
288  y1 = y(:,3,ii) + f1 - 0.5_rp*( f0 + fy(:,1,ii) )
289 
290  !! end, updates
291  y(:,1,ii) = y0 + dt*y1
292  y(:,2,ii) = y0
293  y(:,3,ii) = y1
294  fy(:,1,ii) = f0
295 
296  end do
297  !$OMP END DO
298  !$OMP END PARALLEL
299 
300  end associate
301 
302  end subroutine nl_dc_2
303 
304 
305 
308  subroutine nl_dc_3(sol, dt, t, pb)
309  type(ode_solution), intent(inout) :: sol
310  real(RP) , intent(in) :: dt, t
311  type(ode_problem) , intent(in) :: pb
312 
313  real(RP), dimension(pb%N) :: y0, y1, F0, F1
314  real(RP), dimension(pb%N) :: F2, y2
315  real(RP), dimension(0) :: ay
316  real(RP), parameter :: C1_6 = 1.0_rp/6.0_rp
317  integer :: ii, N
318  real(RP) :: t1, t2, dt_inv
319 
320  n = pb%N
321  t1 = t+dt
322  t2 = t+dt+dt
323  dt_inv = 1.0_rp / dt
324 
325  associate( y=>sol%Y, fy=>sol%BY, x=>pb%X)
326 
327  !$OMP PARALLEL PRIVATE(y0, F0, y1, F1, ay, F2, y2)
328  !$OMP DO
329  do ii=1, sol%dof
330 
331  ! !!!!!!!!!!!!!!!!!!!!!!!
332  ! START DC-2 ar stage n+2
333 
334  ! y0 := Y_0^{n+2}
335  y0 = y(:,5,ii) + dt*fy(:,2,ii)
336 
337  ! F0 := F_0^{n+2}
338  call pb%AB(ay, f0, x(:,ii), t2, y0, n, 0 )
339 
340  ! F1 = F_1^{n+2}
341  y1 = y0 + dt*y(:,6,ii)
342  call pb%AB(ay, f1, x(:,ii), t2, y1, n, 0 )
343 
344  ! y1 = Y_1^{n+2}
345  y1 = y(:,6,ii) + f1 - (f0 + fy(:,2,ii))*0.5_rp
346 
347  ! END DC-2
348  ! !!!!!!!!!!!!!!!!!!!!!!!
349 
350  ! F2 = F_2^{n+1}
351  y2 = y(:,5,ii) + y(:,6,ii)*dt + y(:,4,ii)*dt**2
352  call pb%AB(ay, f2, x(:,ii), t1, y2, n, 0 )
353 
354  ! Computation of y2 = y_2^{n+1}
355  !
356  ! y2 = (1/6) * d^3Y_0^{n+3} * dt**2
357  y2 = f0 - 2.0_rp*fy(:,2,ii) + fy(:,1,ii)
358  y2 = y2 * c1_6
359  !
360  ! y2 = [ (1/6) * d^3Y_0^{n+3} - (1/2) * d^2Y_1^{n+2} ]
361  ! * dt**2
362  y2 = y2 - ( y1 - 2.0_rp*y(:,6,ii) + y(:,3,ii) ) * 0.5_rp
363  !
364  ! y2 = [ (1/6) * d^3Y_0^{n+3} - (1/2) * d^2Y_1^{n+2}
365  ! + F_2^{n+1} - F_1^{n+1} ] * dt**2
366  y2 = y2 + f2 - fy(:,3,ii)
367  !
368  ! y2 = y_2^{n+1}
369  y2 = y(:,4,ii) + y2 * dt_inv
370 
371  !! end, updates
372  y(:,1,ii) = y(:,5,ii) + dt*y(:,6,ii) + dt**2*y2
373  y(:,2,ii) = y(:,5,ii)
374  y(:,3,ii) = y(:,6,ii)
375  y(:,4,ii) = y2
376  y(:,5,ii) = y0
377  y(:,6,ii) = y1
378 
379  fy(:,1,ii) = fy(:,2,ii)
380  fy(:,2,ii) = f0
381  fy(:,3,ii) = f1
382 
383  end do
384  !$OMP END DO
385  !$OMP END PARALLEL
386 
387  end associate
388 
389  end subroutine nl_dc_3
390 
391 end module ode_nl_dc_mod
integer, parameter ode_dc_3
Deferred corrections 3.
Type ode_solution: data structure to solve ODE/PDE problems.
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine init_dc_3(sol, pb, t0, dt)
DERIVED TYPE ode_problem: definition of ODE/PDE problems
subroutine memsize_ode_nl_dc(n_Y, n_FY, method)
required sizes to allocate memory
logical function, public check_ode_method_nl_dc(method)
is 'method' a DC method for NL problems ?
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine nl_dc_2(sol, dt, t, pb)
DC order 2.
subroutine init_nl_dc(sol, pb, meth, t0, dt)
integer, parameter ode_dc_2
Deferred corrections 2.
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
subroutine, public solve_ode_nl_dc(sol, pb, t0, T, dt, meth, out, check_overflow)
solve : DC with constant time step
subroutine nl_dc_3(sol, dt, t, pb)
DC order 3.
DERIVED TYPE ode_solution: data straucture to solve ODEs
CHORAL CONSTANTS
ALGEBRAIC OPERATIONS ON SETS
Definition: algebra_set.F90:14
BOTTOM LEVEL MODULE FOR ODEs
Definition: ode_def.f90:12
integer choral_verb
Verbosity level.
x = x + b*y // OR // z = x + b*y
Definition: R1d_mod.F90:60
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
Type ode_problem: definition of ODE/PDE problems.
DEFERRED CORRECTION SOLVERS FOR NON LINEAR ODEs
DERIVED TYPE krylov: for the resolution of linear systems
Definition: krylov_mod.f90:25
subroutine init_dc_2(sol, pb, t0)
subroutine, public create_ode_nl_dc_sol(sol, pb, method)
Create the solution data structure.