15 real(8),
parameter:: pi2=8d0*atan(1d0)
21 real(8):: a(3), b(3), phi(3), d0, p0, e0
22 real(8):: a0(3), k(3), b0(3), u0(3)
24 character(len=64):: solver
25 procedure(single_solenoidal),
pointer:: condition=>null()
28 procedure,
nopass:: cleanup
30 logical,
save:: first_time=.true.
31 real,
save,
pointer,
dimension(:,:,:,:):: buffer => null()
37 SUBROUTINE init (self, solver, gamma)
39 character(len=64):: solver
40 real,
optional:: gamma
43 integer,
save:: seed=-33
44 logical,
save:: mhd=.true.
45 real,
save:: a0(3)=0.1, k(3)=1.0, d0=1.0, p0=1.0, b0(3)=0.0, u0(3)=0.0
46 character(len=32),
save:: type=
'void' 47 namelist /initial_params/
type, seed, k, a0, u0, b0, d0, p0, mhd
49 call trace_begin(
'initial_t%init')
50 if (
present(gamma)) self%gamma = gamma
56 read (io%input, initial_params, iostat=iostat)
57 mhd = mhd .or. any(b0 /= 0.0)
58 if (io%master)
write (*, initial_params)
60 self%solver = trim(solver)
69 self%gamma = max(self%gamma,1d0+1d-5)
70 self%e0 = p0/(self%gamma-1d0)
75 self%condition => void
76 case (
'single_solenoidal')
77 self%condition => single_solenoidal
78 case (
'single_compressive')
79 self%condition => single_compressive
80 case (
'single_advection')
81 self%condition => single_advection
83 self%condition => exa256
87 print*,
'UNKNOWN INITIAL CONDITION TYPE' 95 SUBROUTINE void (self, m, ff, idx)
98 class(
mesh_t),
pointer,
dimension(:):: m
99 real,
dimension(:,:,:,:) :: ff
105 SUBROUTINE single_solenoidal (self, m, ff, idx)
108 class(
mesh_t),
pointer,
dimension(:):: m
109 real,
dimension(:,:,:,:) :: ff
112 integer :: ix, iy, iz
113 real(8) :: kx, ky, kz
115 call trace_begin(
'initial_t:single_solenoidal')
117 associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
118 do iz=m(3)%lb,m(3)%ub
119 kz = pi2*self%k(3)*(r3(iz)+m(3)%p-0.5d0*m(3)%s)/m(3)%b + self%phi(3)
120 do iy=m(2)%lb,m(2)%ub
121 ky = pi2*self%k(2)*(r2(iy)+m(2)%p-0.5d0*m(2)%s)/m(2)%b + self%phi(2)
122 do ix=m(1)%lb,m(1)%ub
123 kx = pi2*self%k(1)*(r1(ix)+m(1)%p-0.5d0*m(1)%s)/m(1)%b + self%phi(1)
124 ff(ix,iy,iz,idx%d ) = self%d0*(1.0 + self%a0(1)*sin(kx) + self%a0(2)*sin(ky) + self%a0(3)*sin(kz))
125 ff(ix,iy,iz,idx%e ) = self%e0*(1.0 + self%a0(1)*sin(kx) + self%a0(2)*sin(ky) + self%a0(3)*sin(kz))
126 ff(ix,iy,iz,idx%px) = self%u0(1)*(sin(ky)+sin(kz))
127 ff(ix,iy,iz,idx%py) = self%u0(2)*(sin(kz)+sin(kx))
128 ff(ix,iy,iz,idx%pz) = self%u0(3)*(sin(kx)+sin(ky))
130 ff(ix,iy,iz,idx%bx) = self%b0(1)*(sin(ky)+sin(kz))
131 ff(ix,iy,iz,idx%by) = self%b0(2)*(sin(kz)+sin(kx))
132 ff(ix,iy,iz,idx%bz) = self%b0(3)*(sin(kx)+sin(ky))
139 if (self%solver(1:9) ==
'stagger2_')
then 140 associate(e=>ff(:,:,:,idx%e), d=>ff(:,:,:,idx%d))
141 call e2s(e, d, self%gamma)
145 END SUBROUTINE single_solenoidal
150 SUBROUTINE single_compressive (self, m, ff, idx)
153 class(mesh_t),
pointer,
dimension(:):: m
154 real,
dimension(:,:,:,:) :: ff
156 type(random_t) :: random
157 integer :: ix, iy, iz
158 real(8) :: kx, ky, kz
160 call trace_begin(
'initial_t:single_compressive')
162 associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
163 do iz=m(3)%lb,m(3)%ub
164 kz = pi2*self%k(3)*(r3(iz)+m(3)%p-0.5d0*m(3)%s)/m(3)%b + self%phi(3)
165 do iy=m(2)%lb,m(2)%ub
166 ky = pi2*self%k(2)*(r2(iy)+m(2)%p-0.5d0*m(2)%s)/m(2)%b + self%phi(2)
167 do ix=m(1)%lb,m(1)%ub
168 kx = pi2*self%k(1)*(r1(ix)+m(1)%p-0.5d0*m(1)%s)/m(1)%b + self%phi(1)
169 ff(ix,iy,iz,idx%d ) = self%d0
170 ff(ix,iy,iz,idx%e ) = self%e0
171 ff(ix,iy,iz,idx%px) = self%u0(1)*sin(kx)
172 ff(ix,iy,iz,idx%py) = self%u0(2)*sin(ky)
173 ff(ix,iy,iz,idx%pz) = self%u0(3)*sin(kz)
175 ff(ix,iy,iz,idx%bx) = self%b0(1)
176 ff(ix,iy,iz,idx%by) = self%b0(2)
177 ff(ix,iy,iz,idx%bz) = self%b0(3)
184 if (self%solver(1:9) ==
'stagger2_')
then 185 associate(e=>ff(:,:,:,idx%e), d=>ff(:,:,:,idx%d))
186 call e2s(e, d, self%gamma)
190 END SUBROUTINE single_compressive
195 SUBROUTINE single_advection (self, m, ff, idx)
198 class(mesh_t),
pointer,
dimension(:):: m
199 real,
dimension(:,:,:,:) :: ff
201 type(random_t) :: random
202 integer :: ix, iy, iz
203 real(8) :: kx, ky, kz
205 call trace_begin(
'initial_t:single_compressive')
207 associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
208 do iz=m(3)%lb,m(3)%ub
209 kz = pi2*self%k(3)*(r3(iz)+m(3)%p-0.5d0*m(3)%s)/m(3)%b + self%phi(3)
210 do iy=m(2)%lb,m(2)%ub
211 ky = pi2*self%k(2)*(r2(iy)+m(2)%p-0.5d0*m(2)%s)/m(2)%b + self%phi(2)
212 do ix=m(1)%lb,m(1)%ub
213 kx = pi2*self%k(1)*(r1(ix)+m(1)%p-0.5d0*m(1)%s)/m(1)%b + self%phi(1)
214 ff(ix,iy,iz,idx%d ) = self%d0*(1.+self%a0(1)*sin(kx)+self%a0(2)*sin(ky)+self%a0(3)*sin(kz))
215 ff(ix,iy,iz,idx%e ) = self%e0
216 ff(ix,iy,iz,idx%px) = self%d0*self%u0(1)
217 ff(ix,iy,iz,idx%py) = self%d0*self%u0(2)
218 ff(ix,iy,iz,idx%pz) = self%d0*self%u0(3)
220 ff(ix,iy,iz,idx%bx) = self%b0(3)*sin(kz)
221 ff(ix,iy,iz,idx%by) = self%b0(1)*sin(kx)
222 ff(ix,iy,iz,idx%bz) = self%b0(2)*sin(ky)
229 if (self%solver(1:9) ==
'stagger2_')
then 230 associate(e=>ff(:,:,:,idx%e), d=>ff(:,:,:,idx%d))
231 call e2s(e, d, self%gamma)
235 END SUBROUTINE single_advection
240 SUBROUTINE exa256 (self, m, ff, idx)
244 class(mesh_t),
pointer,
dimension(:):: m
245 real,
dimension(:,:,:,:) :: ff
246 real,
dimension(:),
allocatable:: u, d
248 integer :: ix, iy, iz, jy, jz, rec
249 integer(8) :: nx, ny, nz
251 call trace_begin(
'initial_t::exa256')
256 open (io_unit%data, file=
'data/exa256/t0.02.stag', form=
'unformatted', &
257 access=
'direct', recl=m(1)%n*4)
259 ix=m(1)%p/m(1)%d+0.5-m(1)%n/2
260 iy=m(2)%p/m(2)%d+0.5-m(2)%n/2
261 iz=m(3)%p/m(3)%d+0.5-m(3)%n/2
262 nx = m(1)%b/m(1)%d + 0.5
263 ny = m(2)%b/m(2)%d + 0.5
264 nz = m(3)%b/m(3)%d + 0.5
265 do jz=m(3)%li,m(3)%ui
266 do jy=m(2)%li,m(2)%ui
267 rec = 1 + (ix ) /m(1)%n
268 rec = rec + (iy+jy-m(2)%li)*nx /m(1)%n
269 rec = rec + (iz+jz-m(3)%li)*nx*ny/m(1)%n
271 read (io_unit%data,rec=rec) d; ff(m(1)%li:m(1)%ui,jy,jz,idx%d) = d
272 ff(m(1)%li:m(1)%ui,jy,jz,idx%e) = d
273 rec = rec + nx*ny*nz/m(1)%n
274 read (io_unit%data,rec=rec) u; ff(m(1)%li:m(1)%ui,jy,jz,idx%px) = u
275 rec = rec + nx*ny*nz/m(1)%n
276 read (io_unit%data,rec=rec) u; ff(m(1)%li:m(1)%ui,jy,jz,idx%py) = u
277 rec = rec + nx*ny*nz/m(1)%n
278 read (io_unit%data,rec=rec) u; ff(m(1)%li:m(1)%ui,jy,jz,idx%pz) = u
281 if (io%verbose > 0) &
282 print *,
'ix,iy,iz,dmin,dmax =', ix, iy, iz, &
283 minval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1)), &
284 maxval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1))
286 if (self%solver(1:9) ==
'stagger2_')
then 287 associate(e=>ff(:,:,:,idx%e), d=>ff(:,:,:,idx%d))
288 call e2s(e, d, self%gamma)
292 END SUBROUTINE exa256
297 SUBROUTINE raw (self, m, ff, idx)
301 class(mesh_t),
pointer,
dimension(:):: m
302 real,
dimension(:,:,:,:):: ff
304 integer:: n(3), i(3), iv, iw, iostat
306 logical,
save:: first_time=.true.
307 character(len=64),
save:: filename =
'data/exa256/t0.02.stag' 308 namelist /raw_params/ filename
310 call trace_begin(
'initial_t::raw')
321 rewind(io%input);
read (io%input,raw_params,iostat=iostat)
322 if (io%master)
write (*,raw_params)
324 allocate (buffer(n(1),n(2),n(3),mv))
325 open (io_unit%data, file=filename, form=
'unformatted', &
326 access=
'direct', recl=4*product(n))
328 read (io_unit%data, rec=iv) buffer(:,:,:,iv)
332 i = 1 + nint(m%p/m%d) - m%n/2
335 iv = merge(1,iw+1,iw==1)
336 ff(m(1)%li:m(1)%ui, &
338 m(3)%li:m(3)%ui,iv) = buffer(i(1):i(1)+m(1)%n-1, &
339 i(2):i(2)+m(2)%n-1, &
340 i(3):i(3)+m(3)%n-1,iw)
342 ff(:,:,:,idx%e) = ff(:,:,:,idx%d)
343 if (io%verbose > 2) &
344 write (io_unit%log,
'(a,3i4,1p,2g12.3)')
'ix,iy,iz,dmin,dmax =', i, &
345 minval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1)), &
346 maxval(ff(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui,1))
349 if (self%solver(1:9) ==
'stagger2_')
then 350 associate(e=>ff(:,:,:,idx%e), d=>ff(:,:,:,idx%d))
351 call e2s(e, d, self%gamma)
359 if (
associated(buffer))
deallocate(buffer)
360 END SUBROUTINE cleanup
366 SUBROUTINE e2s(s, d, gamma)
369 real,
dimension(:,:,:):: s, d
373 s = d*log(s*g1/d**gamma)/g1
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...
6th order stagger operators, with self-test procedure