17 real(8),
parameter:: pi2=8d0*atan(1d0)
23 character(len=64):: solver
24 procedure(single_solenoidal),
pointer:: condition=>null()
27 procedure,
nopass:: cleanup
29 real,
save,
pointer,
dimension(:,:,:,:):: buffer => null()
30 real,
save:: csound=1.
36 SUBROUTINE init (self, solver, gamma1)
38 character(len=64):: solver
39 character(len=32),
save:: type=
'void' 41 real,
optional:: gamma1
42 logical,
save:: first_time=.true.
43 namelist /initial_params/
type 45 call trace_begin(
'initial_t%init')
50 read (io%input, initial_params, iostat=iostat)
51 if (iostat/=0)
call io%namelist_warning (
'ramses/hydro/initial_mod.f90: initial_params')
52 if (io%master)
write (*, initial_params)
55 self%solver = trim(solver)
57 case (
'single_solenoidal')
58 self%condition => single_solenoidal
59 case (
'single_compressive')
60 self%condition => single_compressive
61 case (
'single_advection')
62 self%condition => single_advection
64 self%condition => exa256
68 self%condition => void
70 print*,
'UNKNOWN INITIAL CONDITION TYPE :'//trim(type)//
':' 77 if (
associated(buffer))
deallocate(buffer)
78 END SUBROUTINE cleanup
83 SUBROUTINE void (self, m, ff, idx)
87 real,
dimension(:,:,:,:) :: ff
93 SUBROUTINE single_solenoidal (self, m, ff, idx)
96 real,
dimension(:,:,:,:) :: ff
100 real(8) :: kx, ky, kz
101 logical,
save :: mhd=.false.
102 real,
save :: a0(3)=0.1, k(3)=1.0, d0=1.0, s0=0.0, b0(3)=0.0, u0(3)=0.0, phi(3)=0.0
103 real(8),
save :: t_turn=0.1
104 logical,
save :: first_time=.true.
105 namelist /initial_solenoidal_params/ k, a0, u0, b0, d0, s0, mhd, gamma, t_turn
106 integer,
save :: itimer=0
108 call trace_begin(
'initial_t:single_solenoidal', itimer=itimer)
114 read (io%input, initial_solenoidal_params)
115 mhd = mhd .and. any(b0 /= 0.0)
116 if (io%master)
write (*, initial_solenoidal_params)
121 associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
122 do iz=m(3)%lb,m(3)%ub
123 kz = pi2*k(3)*(r3(iz)+m(3)%p-0.5d0*m(3)%s)/m(3)%b + phi(3)
124 do iy=m(2)%lb,m(2)%ub
125 ky = pi2*k(2)*(r2(iy)+m(2)%p-0.5d0*m(2)%s)/m(2)%b + phi(2)
126 do ix=m(1)%lb,m(1)%ub
127 kx = pi2*k(1)*(r1(ix)+m(1)%p-0.5d0*m(1)%s)/m(1)%b + phi(1)
128 ff(ix,iy,iz,idx%d ) = d0
129 ff(ix,iy,iz,idx%px) = d0*a0(1)*(sin(ky)+sin(kz))
130 ff(ix,iy,iz,idx%py) = d0*a0(2)*(sin(kz)+sin(kx))
131 ff(ix,iy,iz,idx%pz) = d0*a0(3)*(sin(kx)+sin(ky))
133 ff(ix,iy,iz,idx%bx) = b0(1)
134 ff(ix,iy,iz,idx%by) = b0(2)
135 ff(ix,iy,iz,idx%bz) = b0(3)
140 call energy (self, m, idx, ff)
142 call trace_end (itimer)
143 END SUBROUTINE single_solenoidal
148 SUBROUTINE single_compressive (self, m, ff, idx)
152 real,
dimension(:,:,:,:) :: ff
154 integer :: ix, iy, iz
155 real(8) :: kx, ky, kz
156 logical,
save :: mhd=.false.
157 real,
save :: a0(3)=0.1, k(3)=1.0, d0=1.0, s0=0.0, b0(3)=0.0, u0(3)=0.0, phi(3)=0.0
158 real(8),
save :: t_turn=0.1
159 logical,
save :: first_time=.true.
160 namelist /initial_compressive_params/ k, a0, u0, b0, d0, s0, mhd, gamma, t_turn
161 integer,
save :: itimer=0
163 call trace_begin(
'initial_t:single_compressive', itimer=itimer)
169 read (io%input, initial_compressive_params)
170 mhd = mhd .and. any(b0 /= 0.0)
171 if (io%master)
write (*, initial_compressive_params)
176 associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
177 do iz=m(3)%lb,m(3)%ub
178 kz = pi2*k(3)*(r3(iz)+m(3)%p-0.5d0*m(3)%s)/m(3)%b + phi(3)
179 do iy=m(2)%lb,m(2)%ub
180 ky = pi2*k(2)*(r2(iy)+m(2)%p-0.5d0*m(2)%s)/m(2)%b + phi(2)
181 do ix=m(1)%lb,m(1)%ub
182 kx = pi2*k(1)*(r1(ix)+m(1)%p-0.5d0*m(1)%s)/m(1)%b + phi(1)
183 ff(ix,iy,iz,idx%d ) = d0
184 ff(ix,iy,iz,idx%px) = d0*u0(1)*sin(kx)
185 ff(ix,iy,iz,idx%py) = d0*u0(2)*sin(ky)
186 ff(ix,iy,iz,idx%pz) = d0*u0(3)*sin(kz)
188 ff(ix,iy,iz,idx%bx) = b0(1)
189 ff(ix,iy,iz,idx%by) = b0(2)
190 ff(ix,iy,iz,idx%bz) = b0(3)
195 call energy (self, m, idx, ff)
197 call trace_end (itimer)
198 END SUBROUTINE single_compressive
203 SUBROUTINE single_advection (self, m, ff, idx)
206 class(mesh_t) :: m(3)
207 real,
dimension(:,:,:,:) :: ff
209 integer :: ix, iy, iz
210 real(8) :: kx, ky, kz
211 logical,
save :: mhd=.false.
212 logical,
save :: first_time=.true.
213 real,
save :: a0(3)=0.1, k(3)=1.0, d0=1.0, s0=0.0, b0(3)=0.0, u0(3)=0.0, phi(3)=0.0
214 namelist /initial_advection_params/ k, a0, u0, b0, d0, s0, mhd, gamma
215 integer,
save :: itimer=0
217 call trace_begin(
'initial_t:single_advection', itimer=itimer)
223 read (io%input, initial_advection_params)
224 mhd = mhd .and. any(b0 /= 0.0)
225 if (io%master)
write (*, initial_advection_params)
230 associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
231 do iz=m(3)%lb,m(3)%ub
232 kz = pi2*k(3)*(r3(iz)+m(3)%p-0.5d0*m(3)%s)/m(3)%b + phi(3)
233 do iy=m(2)%lb,m(2)%ub
234 ky = pi2*k(2)*(r2(iy)+m(2)%p-0.5d0*m(2)%s)/m(2)%b + phi(2)
235 do ix=m(1)%lb,m(1)%ub
236 kx = pi2*k(1)*(r1(ix)+m(1)%p-0.5d0*m(1)%s)/m(1)%b + phi(1)
237 ff(ix,iy,iz,idx%d ) = d0*(1.+a0(1)*sin(kx)+a0(2)*sin(ky)+a0(3)*sin(kz))
238 ff(ix,iy,iz,idx%px) = d0*u0(1)
239 ff(ix,iy,iz,idx%py) = d0*u0(2)
240 ff(ix,iy,iz,idx%pz) = d0*u0(3)
242 ff(ix,iy,iz,idx%bx) = b0(3)*sin(kz)
243 ff(ix,iy,iz,idx%by) = b0(1)*sin(kx)
244 ff(ix,iy,iz,idx%bz) = b0(2)*sin(ky)
249 call energy (self, m, idx, ff)
251 call trace_end (itimer)
252 END SUBROUTINE single_advection
257 SUBROUTINE exa256 (self, m, ff, idx)
262 real,
dimension(:,:,:,:):: ff
263 real,
dimension(:),
allocatable:: u, d
265 integer:: ix, iy, iz, jy, jz, rec, iv, li(3), ui(3)
266 integer(8):: nx, ny, nz
268 logical,
save:: first_time=.true.
269 integer,
save:: itimer=0
271 call trace_begin(
'initial_t::exa256', itimer=itimer)
278 open (io_unit%data, file=
'data/exa256/t0.02.raw', form=
'unformatted', &
279 access=
'direct', recl=m(1)%n*4)
282 ix=m(1)%p/m(1)%d+0.5-m(1)%n/2
283 iy=m(2)%p/m(2)%d+0.5-m(2)%n/2
284 iz=m(3)%p/m(3)%d+0.5-m(3)%n/2
285 nx = m(1)%b/m(1)%d + 0.5
286 ny = m(2)%b/m(2)%d + 0.5
287 nz = m(3)%b/m(3)%d + 0.5
293 rec = 1 + (ix ) /m(1)%n
294 rec = rec + (iy+jy-m(2)%li)*nx /m(1)%n
295 rec = rec + (iz+jz-m(3)%li)*nx*ny/m(1)%n
296 read (io_unit%data,rec=rec) d; ff(m(1)%li:m(1)%ui,jy,jz,idx%d ) = d
297 rec = rec + nx*ny*nz/m(1)%n
298 read (io_unit%data,rec=rec) u; ff(m(1)%li:m(1)%ui,jy,jz,idx%px) = u
299 rec = rec + nx*ny*nz/m(1)%n
300 read (io_unit%data,rec=rec) u; ff(m(1)%li:m(1)%ui,jy,jz,idx%py) = u
301 rec = rec + nx*ny*nz/m(1)%n
302 read (io_unit%data,rec=rec) u; ff(m(1)%li:m(1)%ui,jy,jz,idx%pz) = u
303 umax=max(umax,3.+maxval(abs(ff(li(1):ui(1),jy,jz,idx%px)) &
304 +abs(ff(li(1):ui(1),jy,jz,idx%py)) &
305 +abs(ff(li(1):ui(1),jy,jz,idx%pz))))
309 ff(:,:,:,iv) = ff(:,:,:,idx%d)*ff(:,:,:,iv)
311 call energy (self, m, idx, ff)
312 dt = courant_factor*m(1)%d/umax
313 if (io%verbose > 2)
write (io_unit%log,*)
'initial dt:',dt,courant_factor,gamma,umax
314 if (io%verbose > 2) &
315 write (io_unit%log,
'(a,3i4,1p,2g12.3)')
'ix,iy,iz,dmin,dmax =', ix, iy, iz, &
316 minval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1)), &
317 maxval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1))
318 call trace_end (itimer)
319 END SUBROUTINE exa256
325 SUBROUTINE raw (self, m, ff, idx)
330 real,
dimension(:,:,:,:):: ff
332 integer,
parameter:: mv=4
333 integer:: n(3), i(3), iv, jv(mv)
334 logical,
save:: first_time=.true.
335 integer,
save :: itimer=0
337 call trace_begin(
'initial_t::raw', itimer=itimer)
349 allocate (buffer(n(1),n(2),n(3),mv))
350 open (io_unit%data, file=
'data/exa256/t0.02.raw', form=
'unformatted', &
351 access=
'direct', recl=4*product(n))
353 read (io_unit%data, rec=iv) buffer(:,:,:,iv)
357 i = 1 + nint(m%p/m%d) - m%n/2
364 ff(m(1)%li:m(1)%ui, &
366 m(3)%li:m(3)%ui,jv(iv)) = buffer(i(1):i(1)+m(1)%n-1, &
367 i(2):i(2)+m(2)%n-1, &
368 i(3):i(3)+m(3)%n-1,iv)
371 ff(:,:,:,iv) = ff(:,:,:,idx%d)*ff(:,:,:,iv)
373 call energy (self, m, idx, ff)
374 if (io%verbose > 2) &
375 write (io_unit%log,
'(a,3i4,1p,2g12.3)')
'ix,iy,iz,dmin,dmax =', i, &
376 minval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1)), &
377 maxval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1))
379 call trace_end (itimer)
385 SUBROUTINE energy (self, m, idx, ff)
389 real,
dimension(:,:,:,:):: ff
393 call trace%begin (
'initial_t%energy')
394 if (isothermal .or. gamma==1.0_dp)
then 397 fgamma = gamma*(gamma-1.0_dp)
399 ff(:,:,:,idx%e) = 0.5*(ff(:,:,:,idx%px)**2+ &
400 ff(:,:,:,idx%py)**2+ &
401 ff(:,:,:,idx%pz)**2)/ff(:,:,:,idx%d) &
402 + ff(:,:,:,idx%d)/fgamma
404 ff(:,:,:,idx%e) = ff(:,:,:,idx%e) + &
405 0.5*(ff(:,:,:,idx%bx)**2 + &
406 ff(:,:,:,idx%by)**2 + &
415 END SUBROUTINE energy
Template module for mesh.
Simple initical condition module for tests.
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...