DISPATCH
selfgravity_mod.f90
1 !===============================================================================
2 !> Add selfgravity to any solver
3 !===============================================================================
5  USE io_mod
6  USE os_mod
7  USE trace_mod
8  USE patch_mod
9  USE download_mod
10  USE initial_mod
11  USE kinds_mod
12  USE mesh_mod
13  USE units_mod
14  USE poisson_mod
15  USE timer_mod
16  USE omp_timer_mod
17  USE bits_mod
18  USE vector_ops
19  implicit none
20  private
21  type, public:: selfgravity_t
22  type(initial_t) :: initial
23  type(poisson_t) :: poisson
24  logical:: selfgravity=.false.
25  real(8):: d0
26  real(8):: poisson_time=-1.0
27  contains
28  procedure:: pre_init
29  procedure:: init
30  procedure:: pre_update
31  procedure:: post_update
32  procedure:: ahead_of
33  end type
34  real(8), save :: total_mass=0d0, total_volume=0d0
35  real(8), save :: restart_time = 0d0
36  logical, save :: selfgravity=.false., prediction=.true., detailed_timer=.false.
37  integer, save :: verbose=0
38  integer, save :: order=1
39 CONTAINS
40 
41 !===============================================================================
42 !> Further initialisation of the self-gravity module
43 !===============================================================================
44 SUBROUTINE pre_init (self, patch)
45  class(selfgravity_t):: self
46  class(patch_t):: patch
47  integer :: iostat
48  logical, save :: first_time=.true.
49  real, save :: d0=-1.0
50  real(kind=KindScalarVar), pointer :: d(:,:,:)
51  namelist /selfgravity_params/ selfgravity, prediction, order, d0, verbose
52  !-----------------------------------------------------------------------------
53  call trace%begin ('selfgravity_t%pre_init')
54 
55  !$omp critical (input_cr)
56  if (first_time) then
57  first_time = .false.
58  rewind(io%input)
59  read (io%input, selfgravity_params, iostat=iostat)
60  if (io%master) write (io%output, selfgravity_params)
61  end if
62  !$omp end critical (input_cr)
63 
64  self%d0 = d0
65  if (selfgravity) then
66  self%selfgravity = selfgravity
67  if (patch%kind(1:4) /= 'zeus') patch%nv = patch%nv + 1
68  patch%idx%phi = patch%nv
69  if (patch%nw == 1) then
70  patch%nv = patch%nv + 1
71  patch%idx%dphi = patch%nv
72  end if
73  call self%poisson%init
74  else
75  patch%idx%phi = -1
76  end if
77 
78  call trace%end()
79 END SUBROUTINE pre_init
80 
81 !===============================================================================
82 !> Further initialisation of the self-gravity module
83 !===============================================================================
84 SUBROUTINE init (self, patch)
85  class(selfgravity_t):: self
86  class(patch_t):: patch
87  real(kind=KindScalarVar), pointer :: ff(:,:,:,:)
88  logical, save :: first_time=.true., first_poisson=.true.
89  real(kind=KindScalarVar), pointer :: d(:,:,:)
90  !-----------------------------------------------------------------------------
91  call trace%begin ('selfgravity_t%init')
92 
93  if (.not. allocated(patch%force_per_unit_mass)) then
94  if (patch%kind(1:4) /= 'zeus') &
95  allocate (patch%force_per_unit_mass(patch%gn(1),patch%gn(2),patch%gn(3),3))
96  end if
97 
98  !$omp critical (poisson_cr)
99  if (first_poisson) then
100  first_poisson = .false.
101  call download%test
102  end if
103  !$omp end critical (poisson_cr)
104  if (selfgravity) then
105  !-------------------------------------------------------------------------
106  ! Sum up the total mass and total volume, to get the average density d0.
107  ! FIXME: This needs to be modified for MPI
108  !-------------------------------------------------------------------------
109  if (patch%time == 0.0 .and. patch%id==1) then
110  d => patch%mem(:,:,:,patch%idx%d,patch%it,1)
111  !$omp atomic
112  total_mass = total_mass + patch%fsum(d)*product(patch%ds)
113  !$omp atomic
114  total_volume = total_volume + product(patch%size)
115  end if
116  end if
117  call trace%end()
118 END SUBROUTINE init
119 
120 !===============================================================================
121 !> Prepare the external forces for the pde call, in particular selfgravity. At
122 !> this point the nbors are declared up-to-date with respect to everything,
123 !> including their potentials, which we need for guard zones. They may in fact
124 !> be providing only their preliminary potentials in their (nt-1) time slots,
125 !> while their older time slots should by now contain the converged potential.
126 !===============================================================================
127 SUBROUTINE pre_update (self, patch)
128  class(selfgravity_t):: self
129  class(patch_t):: patch
130  !.............................................................................
131  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d, phi
132  real:: dmin, dmax
133  integer:: iter, ix ,iy, iz
134  integer, save:: itimer=0
135  !-----------------------------------------------------------------------------
136  if (.not. selfgravity) return
137  call trace%begin ('selfgravity_t%pre_update', itimer=itimer)
138  d => patch%mem(:,:,:,patch%idx%d ,patch%it,1)
139  phi => patch%mem(:,:,:,patch%idx%phi,patch%it,1)
140  if (self%d0==-1.0) then
141  self%d0 = total_mass/total_volume
142  print *, patch%id, 'setting self%d0 =', self%d0, patch%faver(d)
143  end if
144  dmin = patch%fminval (d)
145  dmax = patch%fmaxval (d)
146  if (dmin==0.0) &
147  print *, patch%id, 'WARNING: selfgravity_t%gravity, dmin, dmax, it =', dmin, dmax, patch%it
148  call self%poisson%update (patch, phi, d, self%d0, detailed_timer)
149  !-----------------------------------------------------------------------------
150  ! Compute the force per unit volume (except for ZEUS solvers, which take the
151  ! gradient of the potential internally.
152  !-----------------------------------------------------------------------------
153  if (patch%kind(1:4) /= 'zeus') &
154  patch%force_per_unit_mass = -grad(patch%mesh, phi)
155  !-----------------------------------------------------------------------------
156  ! compute estimate of dphi/dt
157  !-----------------------------------------------------------------------------
158  call time_derivs (patch)
159  call trace%end (itimer)
160 END SUBROUTINE pre_update
161 
162 !===============================================================================
163 !> Having advanced the MHD state to the new time, we make a preliminary
164 !> call to the Poisson solver, based on extrapolation of the nbor potentials.
165 !> We then provide this preliminary potential to nbors, until we have a
166 !> better offer. The poisson_time controls if an update affects the self%time
167 !> or the poisson_time; the two are updated alternatingly
168 !===============================================================================
169 SUBROUTINE post_update (self, patch)
170  class(selfgravity_t):: self
171  class(patch_t):: patch
172  integer, save:: itimer=0
173  !-----------------------------------------------------------------------------
174  if (.not. selfgravity) return
175  call trace%begin ('selfgravity_t%post_update', itimer=itimer)
176  if (selfgravity) then
177  if (self%poisson_time == patch%time) then
178  if (verbose>0) &
179  print *, patch%id, 'selfgravity_t%update updated patch%time to', patch%time
180  else
181  call patch%dnload (only=patch%idx%d)
182  call patch%dnload (only=patch%idx%phi)
183  call self%pre_update (patch)
184  self%poisson_time = patch%time
185  if (verbose>0) &
186  print *, patch%id, 'selfgravity_t%update updated self%poisson_time to', self%poisson_time
187  end if
188  end if
189  call trace%end(itimer)
190 END SUBROUTINE post_update
191 
192 !===============================================================================
193 !> Compute an approximate, conservative estimate of dphi/dt, which will cause
194 !> the timestep routine to make a good guess at the next phi, at the same time
195 !> as it is updating the other variables. The 1st order prediction is in
196 !> practice better at reducing the number of iterations than the 2nd order.
197 !===============================================================================
198 SUBROUTINE time_derivs (self)
199  class(patch_t):: self
200  integer:: it1, it2, it3, nt, l(3), u(3)
201  real(8):: t1, t2, t3
202  integer, save:: itimer=0
203  real(kind=KindScalarVar), dimension(:,:,:), pointer:: dphidt
204  !-----------------------------------------------------------------------------
205  if (detailed_timer) call trace%begin ('selfgravity_t%time_derivs', itimer=itimer)
206  if (self%nw == 2) then
207  dphidt => self%mem(:,:,:,self%idx%phi,self%it,2)
208  else
209  dphidt => self%mem(:,:,:,self%idx%dphi,self%it,1)
210  end if
211 
212  if (prediction) then
213  nt = self%nt
214  it1 = self%iit(nt-1)
215  it2 = self%iit(nt-2)
216  it3 = self%iit(nt-3)
217  t1 = self%t(it1)
218  t2 = self%t(it2)
219  t3 = self%t(it3)
220  dphidt = 0.0
221  if (it2/=it1 .and. t1/=t2) then
222  l = self%mesh%li
223  u = self%mesh%ui
224  if (order==2 .and. t2/=t3) then
225  dphidt(l(1):u(1),l(2):u(2),l(3):u(3)) = &
226  (self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%phi,it1,1) - &
227  self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%phi,it2,1))/(t1-t2)*1.5 - &
228  (self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%phi,it2,1) - &
229  self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%phi,it3,1))/(t2-t3)*0.5
230  else
231  dphidt(l(1):u(1),l(2):u(2),l(3):u(3)) = &
232  (self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%phi,it1,1) - &
233  self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%phi,it2,1))/(t1-t2)
234  end if
235  if (verbose>2) print *,self%id,'d(phi)/dt min,max,nv =', &
236  self%fminval(dphidt), self%fmaxval(dphidt), size(self%mem,4)
237  end if
238  else
239  dphidt = 0.0
240  end if
241  if (detailed_timer) call trace%end (itimer)
242 END SUBROUTINE time_derivs
243 
244 !===============================================================================
245 !> This functions checks if 'target' is too much ahead of the nbor 'source'.
246 !> If it returns .true. it means the target task is not ready to be updated.
247 !===============================================================================
248 LOGICAL FUNCTION ahead_of (self, target, source)
249  class(selfgravity_t):: self
250  class(patch_t):: target, source
251  !-----------------------------------------------------------------------------
252  ! Level differences larger than 1 should never happen!
253  !-----------------------------------------------------------------------------
254  if (abs(target%level-source%level)>1) then
255  ahead_of = .false.
256  print *,'WARNING: level difference larger than one:', &
257  target%id, target%level, source%id, source%level
258  !-----------------------------------------------------------------------------
259  ! If wee need to solve the Poisson equation, we ignore nbors with higher ot
260  ! the same resolution, and require nbords with lower resolution to be in front
261  ! with their Poisson solutions.
262  !-----------------------------------------------------------------------------
263  else if (self%poisson_time < target%time) then
264  if (target%level <= source%level) then
265  ahead_of = .false.
266  else
267  !ahead_of = target%time > source%poisson_time + source%grace*source%dtime
268  ahead_of = target%time > source%time + source%grace*source%dtime
269  end if
270  if (ahead_of .and. verbose>1) &
271  print 1, target%id, target%level, 'at', target%time, &
272  'cannot update because', source%id, source%level, 'has too old Poisson time', &
273  self%poisson_time + source%grace*source%dtime
274  1 format(i6,i3,2x,a,g15.6,2x,a,i6,i3,2x,a,g15.6)
275  !-----------------------------------------------------------------------------
276  ! If we are going to update the dynamic time, we need the nbor to be in front.
277  !-----------------------------------------------------------------------------
278  else
279  ahead_of = target%time > source%time + source%grace*source%dtime
280  if (ahead_of .and. verbose>1) &
281  print 1, target%id, target%level, 'at', target%time, &
282  'cannot update because', source%id, source%level, 'has too old dynamic time', &
283  source%time + source%grace*source%dtime
284  end if
285 END FUNCTION ahead_of
286 
287 END MODULE selfgravity_mod
Each thread uses a private timer data type, with arrays for start time and total time for each regist...
Definition: timer_mod.f90:11
download_link: takes care of downloads to linktask same: called for patches on the same level differ:...
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Fundamental constants in CGS and SI units.
Definition: units_mod.f90:4
Methods for solving the Poisson equation; partially derived from routines developed by Troels Haugbøll...
Definition: poisson_mod.f90:5
Template module for mesh.
Definition: mesh_mod.f90:4
Simple initical condition module for tests.
Definition: initial_mod.f90:7
Definition: io_mod.f90:4
Add selfgravity to any solver.