Choral
ode_Lin_1s_mod.f90
Go to the documentation of this file.
1 
11 
13 
16  use real_type
17  use basic_tools
18  use r1d_mod, only: xpay, scale, axpby
19  use ode_def
20  use ode_problem_mod
22  use krylov_mod
23 
24  implicit none
25  private
26 
27  public :: solve_ode_lin_1s
28  public :: create_ode_lin_1s_sol
29  public :: set_solver_ode_lin_1s
30  public :: memsize_ode_lin_1s
31  public :: check_ode_method_lin_1s
32 
33 contains
34 
36  function check_ode_method_lin_1s(method) result(b)
37  logical :: b
38  integer, intent(in) :: method
39 
40  b = .false.
41 
42  select case(method)
43  case(ode_be, ode_cn, ode_sdirk4)
44  b = .true.
45 
46  end select
47 
48  end function check_ode_method_lin_1s
49 
50 
57  subroutine memsize_ode_lin_1s(n_V, method)
58  integer, intent(out) :: n_V
59  integer, intent(in) :: method
60 
61  select case(method)
62 
63  case(ode_be, ode_cn)
64  n_v = 1
65 
66  case(ode_sdirk4)
67  n_v = 4
68 
69  case default
70  call quit("ode_Lin_1s_mod: memSize_ode_Lin_1s;&
71  & uncorrect method")
72 
73  end select
74 
75  end subroutine memsize_ode_lin_1s
76 
79  subroutine set_solver_ode_lin_1s(slv, method)
80  procedure(ode_lin_solver), pointer :: slv
81  integer , intent(in) :: method
82 
83  select case(method)
84  case(ode_be)
85  slv => lin_1s_be
86  case(ode_cn)
87  slv => lin_1s_cn
88  case(ode_sdirk4)
89  slv => lin_1s_sdirk4
90 
91  case default
92  call quit("ode_Lin_1s_mod: set_solver_ode_Lin_1s;&
93  & uncorrect method")
94  end select
95 
96  end subroutine set_solver_ode_lin_1s
97 
98 
101  subroutine create_ode_lin_1s_sol(sol, pb, method)
102  type(ode_solution), intent(inout) :: sol
103  type(ode_problem) , intent(in) :: pb
104  integer , intent(in) :: method
105 
106  integer :: n_V
107  logical :: bool
108 
109  bool = check_ode_method_lin_1s(method)
110  if (.NOT.bool) call quit(&
111  & "ode_Lin_1s_mod: create_ode_Lin_1s_sol: uncorrect method")
112  call memsize_ode_lin_1s(n_v, method)
113 
114  sol = ode_solution(pb, nv=n_v)
115 
116  end subroutine create_ode_lin_1s_sol
117 
118 
121  subroutine solve_ode_lin_1s(sol, pb, t0, T, dt, method, &
122  & out, Kinv, kry)
123  type(ode_solution) , intent(inout) :: sol
124  type(ode_problem) , intent(in) :: pb
125  real(RP) , intent(in) :: t0, T, dt
126  integer , intent(in) :: method
127  procedure(ode_output_proc) :: out
128  procedure(linsystem_solver), optional :: Kinv
129  type(krylov) , optional :: kry
130 
131  procedure(ode_lin_solver) , pointer :: slv=>null()
132  procedure(linsystem_solver), pointer :: KInv_1
133  type(krylov) :: kry_1
134  logical :: ierr, exit_comp
135  integer :: ns, jj
136  real(RP) :: tn, CS
137 
138 
139  if (choral_verb>1) write(*,*) &
140  & "ode_Lin_1s_mod : solve_ode_Lin_1s"
141  if (choral_verb>2) write(*,*) &
142  & " Linear onestep solver = ",&
143  & name_ode_method(method)
144  if (choral_verb>2) write(*,*) &
145  & " K_inv provided =",&
146  & present(kinv)
147  if (.NOT.present(kinv)) then
148  if (choral_verb>2) write(*,*) &
149  &" Krylov settings provided = ",&
150  & present(kry)
151  end if
152 
153  !! Set CS to define K = M + CS*S
154  !!
155  cs = s_prefactor(method, dt)
156 
157  !! set the linear system solver for the system K*x = rhs
158  !!
159  if (present(kinv)) then
160  kinv_1 => kinv
161  else
162  kinv_1 => kinv_default
163  if (present(kry)) kry_1 = kry
164  end if
165 
166  !! set the ode solver depending on the method
167  !!
168  call set_solver_ode_lin_1s(slv, method)
169 
170  !! ode resolution
171  !!
172  ns = int( (t-t0) / dt)
173  sol%ierr = 0
174  exit_comp = .false.
175  do jj=1, ns
176  tn = t0 + re(jj-1)*dt
177 
178  ! output
179  call out(tn, sol, exit_comp)
180  if (exit_comp) then
181  if (choral_verb>2) write(*,*) &
182  &"ode_NL_ms_mod : solve&
183  & time = ", tn, ': EXIT_COMPp'
184  return
185  end if
186 
187  call slv(sol, ierr, dt, pb, kinv_1)
188 
189  if (ierr) then
190  sol%ierr = 1
191  return
192  end if
193  end do
194 
195  contains
196 
199  subroutine kinv_default(x, bool, b)
200  real(RP), dimension(:), intent(inout) :: x
201  logical , intent(out) :: bool
202  real(RP), dimension(:), intent(in) :: b
203 
204  call solve(x, kry_1, b, k_default)
205  bool = kry_1%ierr
206 
207  end subroutine kinv_default
208 
209  !! matrix-vector product x --> K*x, K = M + CS*S
210  !!
211  subroutine k_default(y, x)
212  real(RP), dimension(:), intent(out) :: y
213  real(RP), dimension(:), intent(in) :: x
214 
215  call pb%M(y , x)
216  call pb%S(sol%aux, x)
217  call xpay(y, cs, sol%aux)
218 
219  end subroutine k_default
220 
221  end subroutine solve_ode_lin_1s
222 
223 
226  subroutine lin_1s_be(sol, ierr, dt, pb, KInv)
227  type(ode_solution) , intent(inout) :: sol
228  logical , intent(out) :: ierr
229  real(RP) , intent(in) :: dt
230  type(ode_problem) , intent(in) :: pb
231  procedure(linsystem_solver) :: KInv
232 
233  call pb%M(sol%rhs, sol%V(:,1))
234  call kinv(sol%V(:,1), ierr, sol%rhs)
235 
236  end subroutine lin_1s_be
237 
238 
241  subroutine lin_1s_cn(sol, ierr, dt, pb, KInv)
242  type(ode_solution) , intent(inout) :: sol
243  logical , intent(out) :: ierr
244  real(RP) , intent(in) :: dt
245  type(ode_problem) , intent(in) :: pb
246  procedure(linsystem_solver) :: KInv
247 
248  call pb%M(sol%rhs, sol%V(:,1))
249  call pb%S(sol%aux, sol%V(:,1))
250  call xpay(sol%rhs, -dt/2._rp, sol%aux)
251 
252  call kinv(sol%V(:,1), ierr, sol%rhs)
253 
254  end subroutine lin_1s_cn
255 
256 
259  subroutine lin_1s_sdirk4(sol, ierr, dt, pb, KInv)
260  type(ode_solution) , intent(inout) :: sol
261  logical , intent(out) :: ierr
262  real(RP) , intent(in) :: dt
263  type(ode_problem) , intent(in) :: pb
264  procedure(linsystem_solver) :: KInv
265 
266  real(RP) :: gamma, a
267 
268  gamma = 1._rp/sqrt(3._rp) * cos(pi/18._rp)
269  gamma = gamma + 0.5_rp
270 
271  ! sol%V(:,2) := k1
272  call scale(sol%aux, -1._rp, sol%V(:,1))
273  call pb%S(sol%rhs, sol%aux)
274  call kinv(sol%V(:,2), ierr, sol%rhs)
275  if (ierr) return
276 
277  ! sol%V(:,3) := k2
278  a = dt * (0.5_rp - gamma)
279  call axpby(sol%aux, -1._rp, sol%V(:,1), -a, sol%V(:,2))
280  call pb%S(sol%rhs, sol%aux)
281  call kinv(sol%V(:,3), ierr, sol%rhs)
282  if (ierr) return
283 
284  ! sol%V(:,4) := k3
285  a = dt * 2._rp * gamma
286  call axpby(sol%aux, -1._rp, sol%V(:,1), -a, sol%V(:,2))
287  a = dt * (1._rp - 4._rp*gamma)
288  call xpay(sol%aux, -a, sol%V(:,3))
289  call pb%S(sol%rhs, sol%aux)
290  call kinv(sol%V(:,4), ierr, sol%rhs)
291  if (ierr) return
292 
293  a = 6._rp * ( 2._rp*gamma - 1._rp )**2
294  a = dt / a
295  call xpay(sol%V(:,1), a, sol%V(:,2) )
296  call xpay(sol%V(:,1), a, sol%V(:,4) )
297 
298  a = dt - 2._rp*a
299  call xpay(sol%V(:,1), a, sol%V(:,3) )
300 
301  end subroutine lin_1s_sdirk4
302 
303 
304 end module ode_lin_1s_mod
integer, parameter ode_be
Backward Euler.
Type ode_solution: data structure to solve ODE/PDE problems.
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine k_default(y, x)
DERIVED TYPE ode_problem: definition of ODE/PDE problems
subroutine lin_1s_sdirk4(sol, ierr, dt, pb, KInv)
SDIRK4 HAIRER II p. 100.
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
x = a*x // OR // y = a*x
Definition: R1d_mod.F90:49
The type krylov defines the settings of a linear solver.
Definition: krylov_mod.f90:63
integer, parameter ode_cn
Crank Nicholson.
subroutine, public quit(message)
Stop program execution, display an error messahe.
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
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
subroutine, public solve_ode_lin_1s(sol, pb, t0, T, dt, method, out, Kinv, kry)
solve with constant time-step
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
ONE-STEP SOLVERS FOR LINEAR ODEs
DERIVED TYPE ode_solution: data straucture to solve ODEs
subroutine kinv_default(x, bool, b)
Default solver for K*x = b.
subroutine, public set_solver_ode_lin_1s(slv, method)
Setting the solver for diffusion.
CHORAL CONSTANTS
real(rp), parameter, public pi
REAL CONSTANT Pi.
Definition: real_type.F90:96
subroutine lin_1s_cn(sol, ierr, dt, pb, KInv)
Crank - Nicolson.
subroutine lin_1s_be(sol, ierr, dt, pb, KInv)
BacWard Euler.
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
z = a*x + b*y
Definition: R1d_mod.F90:66
subroutine, public create_ode_lin_1s_sol(sol, pb, method)
allocate memory for the ode_solution structure 'sol'
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
logical function, public check_ode_method_lin_1s(method)
is 'method' a one-step linear ODE solver ?
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_sdirk4
SDIRK4.