DISPATCH
timestep_mod.f90
1 !===============================================================================
2 !> $Id: 3d926374bcd48d97b7549e8c1ab5ccb6d73881b2 $
3 !> Adams-Bashforth time stepping, with constant time step assumed (for now)
4 !>
5 !> We store results in nt memory slots, starting with number 1, and continuing
6 !> with numbers 2..nt, and then looping back to 1, but this is mapped by the
7 !> iit array so that the last nt-1 time slot are t(iit(1:nt-1)), and t(iit(nt))
8 !> is the next time slot. The (1:nt-1) slots are complete, and may be used
9 !> freely, while the (nt) slot is the one where new values are to be stored.
10 !> The experiment_t%update procedure thus writes new time derivatives
11 !> there, and it is the task of this timestep procedure to update the values
12 !> there, using them time derivatives and the previous values.
13 !===============================================================================
15  USE io_mod
16  USE trace_mod
17  USE kinds_mod
18  USE mesh_mod
19  USE omp_lock_mod
20  implicit none
21  private
22  integer, parameter:: nslice = 4
23  type, public:: timestep_t
24  integer:: verbose=0 ! verbosity
25  integer:: time_order=-1 ! time order of integration
26  real:: wt(nslice) ! time slot weights
27  contains
28  procedure:: init
29  procedure:: update
30  procedure:: update_single
31  end type
32  type(timestep_t), save, public:: timestep
33 CONTAINS
34 
35 !===============================================================================
36 !> Initialize coefficients for Adams-Bashforth time stepping
37 !===============================================================================
38 SUBROUTINE init (self)
39  class(timestep_t):: self
40  integer, save:: verbose=0, time_order=-1
41  namelist /timestep_params/ verbose, time_order
42  !.............................................................................
43  call trace_begin ('timestep_t%init', 1)
44  !$omp critical (timestep_cr)
45  if (time_order==-1) then
46  time_order = 2
47  rewind(io%input)
48  read (io%input,timestep_params)
49  if (io%master) write (*,timestep_params)
50  end if
51  !$omp end critical (timestep_cr)
52  self%verbose = verbose
53  self%time_order = time_order
54  !---------------------------------------------------------------------------
55  select case (time_order)
56  case (0)
57  case (1)
58  self%wt = [ 1.0, 0.0, 0.0, 0.0] ! 1st order
59  case (2)
60  self%wt = [ 1.5, -.5, 0.0, 0.0] ! 2nd order
61  case (3)
62  self%wt = [23./12., -4./3., 5./12., 0.0] ! 3rd order
63  case (4)
64  self%wt = [55./24., -59./24., 37/24., -3./8.] ! 4th order
65  case default
66  if (io%master) print*,'ERROR: illegal timestep time_order =', time_order
67  stop
68  end select
69  call trace_end
70 END SUBROUTINE init
71 
72 !===============================================================================
73 !> Adams-Bashforth linear multi-step update. We consider it to be linear
74 !> multi-step in time index space, which may be just as good -- or better --
75 !> than attempting to map into linear time space. The time slot derivative
76 !> values have already been multiplied with the individual time steps dt.
77 !===============================================================================
78 SUBROUTINE update (self, id, iit, dt, mem, mesh, lock)
79  class(timestep_t):: self
80  class(mesh_t):: mesh(:)
81  type(lock_t):: lock
82  integer:: id, nt
83  integer, dimension(:):: iit
84  real(8), dimension(:):: dt
85  real(kind=KindScalarVar), dimension(:,:,:,:,:,:):: mem ! memory buffer
86  integer:: l(3), u(3)
87  integer, save:: itimer=0
88  !......................................................................
89  call trace_begin ('timestep_t%update', itimer=itimer)
90  if (self%time_order==-1) call io%abort('timestep_mod not initialized!')
91  nt = size(iit) ! number of time slots
92  l = mesh%lb
93  u = mesh%ub
94  where (mesh%lower_boundary)
95  l = mesh%lb
96  end where
97  where (mesh%upper_boundary)
98  u = mesh%ub
99  end where
100  if (size(mem,6)==1.and.self%time_order>0) then
101  if (io%master) print *,id,'WARNING: setting time_order=0 in timestep%update'
102  self%time_order=0
103  end if
104  call lock%set ('timestep')
105  select case (self%time_order)
106  case (0)
107  case (1)
108  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt ),1) = &
109  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),1) + &
110  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),2)*(self%wt(1)*dt(iit(nt-1)))
111  case (2)
112  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt ),1) = &
113  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),1) + &
114  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),2)*(self%wt(1)*dt(iit(nt-1))) + &
115  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-2),2)*(self%wt(2)*dt(iit(nt-2)))
116  case (3)
117  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt ),1) = &
118  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),1) + &
119  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),2)*(self%wt(1)*dt(iit(nt-1))) + &
120  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-2),2)*(self%wt(2)*dt(iit(nt-2))) + &
121  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-3),2)*(self%wt(3)*dt(iit(nt-3)))
122  case (4)
123  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt ),1) = &
124  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),1) + &
125  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-1),2)*(self%wt(1)*dt(iit(nt-1))) + &
126  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-2),2)*(self%wt(2)*dt(iit(nt-2))) + &
127  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-3),2)*(self%wt(3)*dt(iit(nt-3))) + &
128  mem(l(1):u(1),l(2):u(2),l(3):u(3),:,iit(nt-4),2)*(self%wt(4)*dt(iit(nt-4)))
129  end select
130  if (self%verbose>3) then
131  print'(a,i4,2x,a,1p,6e10.2,2x,6i5)', 'timestep, id:', &
132  id, 'minmax:', &
133  minval(mem(l(1):u(1),l(2):u(2),l(3):u(3),1,iit(nt-1),2)), &
134  maxval(mem(l(1):u(1),l(2):u(2),l(3):u(3),1,iit(nt-1),2)), &
135  minval(mem(l(1):u(1),l(2):u(2),l(3):u(3),1,iit(nt-1),1)), &
136  maxval(mem(l(1):u(1),l(2):u(2),l(3):u(3),1,iit(nt-1),1)), &
137  minval(mem(l(1):u(1),l(2):u(2),l(3):u(3),1,iit(nt ),1)), &
138  maxval(mem(l(1):u(1),l(2):u(2),l(3):u(3),1,iit(nt ),1)), shape(mem)
139  end if
140  call lock%unset ('timestep')
141  call trace_end (itimer)
142 END SUBROUTINE update
143 
144 !===============================================================================
145 SUBROUTINE update_single (self, id, iit, dt, mem, dmem, mesh, lock, time_order)
146  class(timestep_t):: self
147  class(mesh_t):: mesh(:)
148  type(lock_t):: lock
149  integer:: id, nt
150  integer, dimension(:):: iit
151  real(8), dimension(:):: dt
152  real(kind=KindScalarVar), dimension(:,:,:,:):: mem, dmem ! memory buffers
153  integer, intent(in):: time_order
154  integer:: l(3), u(3)
155  real:: wt(nslice)
156  integer, save:: itimer=0
157  !......................................................................
158  call trace_begin ('timestep_t%update_single', itimer=itimer)
159  if (time_order <= 0) return
160 
161  nt = size(iit) ! number of time slots
162  l = mesh%lb
163  u = mesh%ub
164  where (mesh%lower_boundary)
165  l = mesh%lb
166  end where
167  where (mesh%upper_boundary)
168  u = mesh%ub
169  end where
170 
171  select case (time_order)
172  case (0)
173  case (1)
174  wt = [ 1.0, 0.0, 0.0, 0.0] ! 1st order
175  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt )) = &
176  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1)) + &
177  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))*(wt(1)*dt(iit(nt-1)))
178  case (2)
179  wt = [ 1.5, -.5, 0.0, 0.0] ! 2nd order
180  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt )) = &
181  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1)) + &
182  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))*(wt(1)*dt(iit(nt-1))) + &
183  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-2))*(wt(2)*dt(iit(nt-2)))
184  case (3)
185  wt = [23./12., -4./3., 5./12., 0.0] ! 3rd order
186  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt )) = &
187  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1)) + &
188  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))*(wt(1)*dt(iit(nt-1))) + &
189  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-2))*(wt(2)*dt(iit(nt-2))) + &
190  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-3))*(wt(3)*dt(iit(nt-3)))
191  case (4)
192  wt = [55./24., -59./24., 37/24., -3./8.] ! 4th order
193  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt )) = &
194  mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1)) + &
195  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))*(wt(1)*dt(iit(nt-1))) + &
196  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-2))*(wt(2)*dt(iit(nt-2))) + &
197  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-3))*(wt(3)*dt(iit(nt-3))) + &
198  dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-4))*(wt(4)*dt(iit(nt-4)))
199  end select
200  if (self%verbose>3) then
201  print'(a,i4,2x,a,1p,6e10.2,2x,6i5)', 'timestep, id:', &
202  id, 'minmax:', &
203  minval(dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))), &
204  maxval(dmem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))), &
205  minval(mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))), &
206  maxval(mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt-1))), &
207  minval(mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt ))), &
208  maxval(mem(l(1):u(1),l(2):u(2),l(3):u(3),iit(nt ))), shape(mem)
209  end if
210  call trace_end (itimer)
211 END SUBROUTINE update_single
212 
213 END MODULE timestep_mod
Template module for mesh.
Definition: mesh_mod.f90:4
Definition: io_mod.f90:4
The lock module uses nested locks, to allow versatile use of locks, where a procedure may want to mak...