Choral
ode_SL_NL_ms_mod.f90
Go to the documentation of this file.
1 
10 
11 
13 
15  use real_type
16  use basic_tools
17  use algebra_set
18  use r1d_mod, only: xpay
19  use ode_def
20  use ode_problem_mod
22  use ode_nl_ms_mod
23  use krylov_mod
24  use ode_sl_ms_mod
25  !$ use OMP_LIB
26 
27  implicit none
28  private
29 
30  ! %----------------------------------------%
31  ! | |
32  ! | PUBLIC DATA |
33  ! | |
34  ! %----------------------------------------%
35 
36  public :: solve_ode_sl_nl_ms
37  public :: create_ode_sl_nl_ms_sol
38 
39 contains
40 
43  subroutine create_ode_sl_nl_ms_sol(sol, pb, SL_meth, NL_meth)
44  type(ode_solution), intent(inout) :: sol
45  type(ode_problem) , intent(in) :: pb
46  integer , intent(in) :: SL_meth, NL_meth
47 
48  integer :: nV, nY, nFY, nFY2
49 
50  call memsize_ode_nl_ms(ny, nfy , nl_meth)
51  call memsize_ode_sl_ms(nv, nfy2, sl_meth)
52  nfy = max(nfy, nfy2)
53  sol = ode_solution(pb, nv=nv, ny=ny, nfy=nfy)
54 
55  end subroutine create_ode_sl_nl_ms_sol
56 
57 
60  subroutine solve_ode_sl_nl_ms(sol, pb, t0, T, dt, &
61  & SL_meth, NL_meth, out, check_overflow, Kinv, kry)
62  type(ode_solution) , intent(inout) :: sol
63  type(ode_problem) , intent(in) :: pb
64  real(RP) , intent(in) :: t0, T, dt
65  integer , intent(in) :: SL_meth, NL_meth
66  procedure(ode_output_proc) :: out
67  logical , intent(in) :: check_overflow
68  procedure(linsystem_solver), optional :: Kinv
69  type(krylov) , optional :: kry
70 
71  procedure(ode_nl_ms_solver), pointer :: slv_NL
72  procedure(ode_lin_solver) , pointer :: slv_SL
73  procedure(linsystem_solver), pointer :: KInv_1
74  type(krylov) :: kry_1
75  real(RP) :: tn, CS
76  real(DP) :: t_NL, t_SL, t_out, t_reac, cpu, t_tot
77  integer :: nS, ii, jj, N, Na, P, n0, n0_F, n0_V
78  logical :: exit_comp, ierr
79 
80  if (choral_verb>1) write(*,*) &
81  &"ode_SL_NL_ms_mod: solve"
82  if (choral_verb>2) write(*,*) &
83  &" SemiLin multistep solver = ",&
84  & name_ode_method(sl_meth)
85  if (choral_verb>2) write(*,*) &
86  &" NonLin multistep solver = ",&
87  & name_ode_method(nl_meth)
88  if (choral_verb>2) write(*,*) &
89  &" K_inv provided =",&
90  & present(kinv)
91  if (.NOT.present(kinv)) then
92  if (choral_verb>2) write(*,*) &
93  &" Krylov settings provided = ",&
94  & present(kry)
95  end if
96 
97  !! Set CS to define K = M + CS*S
98  !!
99  cs = s_prefactor(sl_meth, dt)
100 
101  !! set the linear system solver for the system K*x = rhs
102  !!
103  if (present(kinv)) then
104  kinv_1 => kinv
105  else
106  kinv_1 => kinv_default
107  if (present(kry)) kry_1 = kry
108  end if
109 
110  !! set the elementary solvers
111  !!
112  call set_solver_ode_sl_ms(slv_sl, sl_meth)
113  call set_solver_ode_nl_ms(slv_nl, nl_meth)
114 
115  !! initialise the solution indexes
116  !!
117  call ode_solution_init_indexes(sol)
118 
119  !! timer initialisation
120  t_nl = 0.0_dp
121  t_sl = 0.0_dp
122  t_out = 0.0_dp
123  t_reac = 0.0_dp
124  t_tot = clock()
125 
126  !! ode resolution
127  !!
128  n = pb%N
129  na = pb%Na
130  if (na < n) then
131  p = na
132  else
133  p = n-1
134  end if
135  ns = int( (t-t0) / dt)
136  exit_comp = .false.
137  sol%ierr = 0
138  !!
139  !! Time loop
140  !!
141  do jj=1, ns
142  tn = t0 + re(jj-1)*dt
143 
144  ! output
145  cpu = clock()
146  call out(tn, sol, exit_comp)
147  if (exit_comp) then
148  if (choral_verb>2) write(*,*) &
149  &"ode_SL_NL_ms_mod: solve&
150  & time = ", tn, ': EXIT_COMPp'
151  return
152  end if
153  t_out = t_out + clock() - cpu
154 
155  !! solve the non linear ODE system
156  cpu = clock()
157  call slv_nl(sol, dt, n, p)
158  t_nl = t_nl + clock() - cpu
159 
160  ! solve the semilinear equation
161  cpu = clock()
162  call slv_sl(sol, ierr, dt, pb, kinv_1)
163  t_sl = t_sl + clock() - cpu
164 
165  !! check linear system resolution
166  !!
167  if (ierr) then
168  sol%ierr = 1
169  if (choral_verb>2) write(*,*) &
170  &"ode_SL_NL_ms_mod: solve,&
171  & time = ", tn, ': lin pb inv error'
172  return
173  end if
174 
175  !! updates
176  !!
177  cpu = clock()
178  associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY, &
179  & x=>pb%X, v=>sol%V)
180 
181  call circperm(sol%Y_i)
182  call circperm(sol%F_i)
183  call circperm(sol%V_i)
184 
185  n0 = sol%Y_i(1)
186  n0_f = sol%F_i(1)
187  n0_v = sol%V_i(1)
188 
189  ! 1. recast V into Y
190  !
191  ! 2. compute reaction terms
192  !
193  if (na<n) then
194  !$OMP PARALLEL
195  !$OMP DO
196  do ii=1, sol%dof
197  y(n, n0, ii) = v(ii, n0_v)
198 
199  call pb%AB( ay(:, n0_f, ii), by(:, n0_f, ii), &
200  & x(:,ii), tn+dt, y(:, n0, ii), n, na )
201  end do
202  !$OMP END DO
203  !$OMP END PARALLEL
204  else
205  !$OMP PARALLEL
206  !$OMP DO
207  do ii=1, sol%dof
208  y(n, n0, ii) = v(ii, n0_v)
209 
210  call pb%AB( ay(:, n0_f, ii), by(:, n0_f, ii), &
211  & x(:,ii), tn+dt, y(:, n0, ii), n, na )
212 
213  by(n, n0_f, ii) = by(n, n0_f, ii) &
214  & + ay(n, n0_f, ii)*v(ii, n0_v)
215  end do
216  !$OMP END DO
217  !$OMP END PARALLEL
218  end if
219 
220  !! Check overflow
221  !!
222  if (check_overflow) then
223  ierr = overflow(sol%Y(:,n0,:))
224  if (ierr) then
225  sol%ierr = 10
226  if (choral_verb>2) write(*,*) &
227  &"ode_NL_ms_mod : solve,&
228  & time = ", tn, ': OVERFLOW Y'
229  return
230  end if
231 
232  ierr = overflow(sol%BY(:,n0_f,:))
233  if (ierr) then
234  sol%ierr = 10
235  if (choral_verb>2) write(*,*) &
236  &"ode_NL_ms_mod : solve,&
237  & time = ", tn, ': OVERFLOW BY'
238  return
239  end if
240 
241  ierr = overflow(sol%AY(:,n0_f,:))
242  if (ierr) then
243  sol%ierr = 10
244  if (choral_verb>2) write(*,*) &
245  &"ode_NL_ms_mod : solve,&
246  & time = ", tn, ': OVERFLOW AY'
247  return
248  end if
249  end if
250 
251  end associate
252  t_reac = t_reac + clock() - cpu
253 
254  end do
255 
256  if (choral_verb>1) then
257  write(*,*)"ode_SL_NL_ms_mod: solve, end = timing"
258  write(*,*)" Non Lin. system integration =", &
259  & real(t_nl, sp)
260  write(*,*)" Semilin. eq. integration =", &
261  & real(t_sl, sp)
262  write(*,*)" Reaction term evaluation =", &
263  & real(t_reac, sp)
264  write(*,*)" Output =", &
265  & real(t_out, sp)
266  t_tot = clock() - t_tot
267  write(*,*)" Total CPU =", &
268  & real(t_tot, sp)
269 
270  end if
271 
272  contains
273 
276  subroutine kinv_default(x, bool, b)
277  real(RP), dimension(:), intent(inout) :: x
278  logical , intent(out) :: bool
279  real(RP), dimension(:), intent(in) :: b
280 
281  call solve(x, kry_1, b, k_default)
282  bool = kry_1%ierr
283 
284  end subroutine kinv_default
285 
286  !! matrix-vector product x --> K*x, K = M + CS*S
287  !!
288  subroutine k_default(y, x)
289  real(RP), dimension(:), intent(out) :: y
290  real(RP), dimension(:), intent(in) :: x
291 
292  call pb%M(y , x)
293  call pb%S(sol%aux, x)
294  call xpay(y, cs, sol%aux)
295 
296  end subroutine k_default
297 
298  end subroutine solve_ode_sl_nl_ms
299 
300 
301 end module ode_sl_nl_ms_mod
MULTISTEP SOLVERS FOR A SEMILINEAR ODE coupled with a non-linear ODE system
Type ode_solution: data structure to solve ODE/PDE problems.
BASIC TOOLS
Definition: basic_tools.F90:8
double precision function, public clock()
Returns the internal clock time
Definition: basic_tools.F90:94
subroutine k_default(y, x)
DERIVED TYPE ode_problem: definition of ODE/PDE problems
subroutine, public create_ode_sl_nl_ms_sol(sol, pb, SL_meth, NL_meth)
Create the solution data structure.
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 set_solver_ode_sl_ms(slv, method)
set the resolution solver
integer, parameter, public sp
simple precision for real numbers
Definition: real_type.F90:60
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 solve_ode_sl_nl_ms(sol, pb, t0, T, dt, SL_meth, NL_meth, out, check_overflow, Kinv, kry)
solve : multistep with constant time step
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
DERIVED TYPE ode_solution: data straucture to solve ODEs
subroutine kinv_default(x, bool, b)
Default solver for K*x = b.
subroutine, public memsize_ode_nl_ms(n_Y, n_FY, method)
required sizes to allocate memory
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
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 set_solver_ode_nl_ms(slv, method)
set the resolution solver
MULTISTEP SOLVERS FOR NON LINEAR ODEs
Type ode_problem: definition of ODE/PDE problems.
DERIVED TYPE krylov: for the resolution of linear systems
Definition: krylov_mod.f90:25