24 logical:: selfgravity=.false.
26 real(8):: poisson_time=-1.0
30 procedure:: pre_update
31 procedure:: post_update
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
44 SUBROUTINE pre_init (self, patch)
48 logical,
save :: first_time=.true.
50 real(kind=KindScalarVar),
pointer :: d(:,:,:)
51 namelist /selfgravity_params/ selfgravity, prediction, order, d0, verbose
53 call trace%begin (
'selfgravity_t%pre_init')
59 read (io%input, selfgravity_params, iostat=iostat)
60 if (io%master)
write (io%output, selfgravity_params)
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
73 call self%poisson%init
79 END SUBROUTINE pre_init
84 SUBROUTINE init (self, patch)
87 real(kind=KindScalarVar),
pointer :: ff(:,:,:,:)
88 logical,
save :: first_time=.true., first_poisson=.true.
89 real(kind=KindScalarVar),
pointer :: d(:,:,:)
91 call trace%begin (
'selfgravity_t%init')
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))
99 if (first_poisson)
then 100 first_poisson = .false.
104 if (selfgravity)
then 109 if (patch%time == 0.0 .and. patch%id==1)
then 110 d => patch%mem(:,:,:,patch%idx%d,patch%it,1)
112 total_mass = total_mass + patch%fsum(d)*product(patch%ds)
114 total_volume = total_volume + product(patch%size)
127 SUBROUTINE pre_update (self, patch)
131 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: d, phi
133 integer:: iter, ix ,iy, iz
134 integer,
save:: itimer=0
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)
144 dmin = patch%fminval (d)
145 dmax = patch%fmaxval (d)
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)
153 if (patch%kind(1:4) /=
'zeus') &
154 patch%force_per_unit_mass = -
grad(patch%mesh, phi)
158 call time_derivs (patch)
159 call trace%end (itimer)
160 END SUBROUTINE pre_update
169 SUBROUTINE post_update (self, patch)
172 integer,
save:: itimer=0
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 179 print *, patch%id,
'selfgravity_t%update updated patch%time to', patch%time
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
186 print *, patch%id,
'selfgravity_t%update updated self%poisson_time to', self%poisson_time
189 call trace%end(itimer)
190 END SUBROUTINE post_update
198 SUBROUTINE time_derivs (self)
200 integer:: it1, it2, it3, nt, l(3), u(3)
202 integer,
save:: itimer=0
203 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: dphidt
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)
209 dphidt => self%mem(:,:,:,self%idx%dphi,self%it,1)
221 if (it2/=it1 .and. t1/=t2)
then 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
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)
235 if (verbose>2) print *,self%id,
'd(phi)/dt min,max,nv =', &
236 self%fminval(dphidt), self%fmaxval(dphidt),
size(self%mem,4)
241 if (detailed_timer)
call trace%end (itimer)
242 END SUBROUTINE time_derivs
248 LOGICAL FUNCTION ahead_of (self, target, source)
250 class(
patch_t)::
target, source
254 if (abs(
target%level-source%level)>1)
then 256 print *,
'WARNING: level difference larger than one:', &
257 target%id,
target%level, source%id, source%level
263 else if (self%poisson_time <
target%time)
then 264 if (
target%level <= source%level)
then 268 ahead_of =
target%time > source%time + source%grace*source%dtime
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)
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
285 END FUNCTION ahead_of
Each thread uses a private timer data type, with arrays for start time and total time for each regist...
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...
Fundamental constants in CGS and SI units.
Methods for solving the Poisson equation; partially derived from routines developed by Troels Haugbøll...
Template module for mesh.
Simple initical condition module for tests.
Add selfgravity to any solver.