DISPATCH
initial_mod.f90
1 !===============================================================================
2 !> This is a start on a module that could become a general turbulence forcing
3 !> module, where the standard forcing we have been using in RAMSES is implemented
4 !===============================================================================
5 MODULE initial_mod
6  USE io_mod
7  USE trace_mod
8  USE vector_mod
9  USE mesh_mod
10  USE random_mod
11  USE stagger_mod
12  USE index_mod
13  implicit none
14  private
15  real(8), parameter:: pi2=8d0*atan(1d0)
16  type, public:: initial_t
17  integer:: id=0
18  integer:: seed
19  logical:: mhd
20  real(8):: gamma=0d0
21  real(8):: a(3), b(3), phi(3), d0, p0, e0
22  real(8):: a0(3), k(3), b0(3), u0(3)
23  real(8):: time
24  character(len=64):: solver
25  procedure(single_solenoidal), pointer:: condition=>null()
26  contains
27  procedure:: init
28  procedure, nopass:: cleanup
29  end type
30  logical, save:: first_time=.true.
31  real, save, pointer, dimension(:,:,:,:):: buffer => null()
32 CONTAINS
33 
34 !===============================================================================
35 !> Setup a uniform initial state, with density=d0 and B=B0
36 !===============================================================================
37 SUBROUTINE init (self, solver, gamma)
38  class(initial_t):: self
39  character(len=64):: solver
40  real, optional:: gamma
41  !............................................................................
42  integer:: iostat
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
48  !----------------------------------------------------------------------------
49  call trace_begin('initial_t%init')
50  if (present(gamma)) self%gamma = gamma
51  !$omp critical (input_cr)
52  if (first_time) then
53  first_time = .false.
54  mhd = self%mhd
55  rewind(io%input)
56  read (io%input, initial_params, iostat=iostat)
57  mhd = mhd .or. any(b0 /= 0.0)
58  if (io%master) write (*, initial_params)
59  end if
60  self%solver = trim(solver)
61  !$omp end critical (input_cr)
62  self%seed = seed
63  self%k = k
64  self%a0 = a0
65  self%b0 = b0
66  self%u0 = u0
67  self%d0 = d0
68  self%p0 = p0
69  self%gamma = max(self%gamma,1d0+1d-5)
70  self%e0 = p0/(self%gamma-1d0)
71  self%mhd = mhd
72  self%time = 0d0
73  select case (type)
74  case ('void')
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
82  case ('exa256')
83  self%condition => exa256
84  case ('raw')
85  self%condition => raw
86  case default
87  print*,'UNKNOWN INITIAL CONDITION TYPE'
88  end select
89  call trace_end
90 END SUBROUTINE init
91 
92 !===============================================================================
93 !> Do-nothing IC; when initializing in experiment_t%init()
94 !===============================================================================
95 SUBROUTINE void (self, m, ff, idx)
96  class(initial_t) :: self
97  type(index_t):: idx
98  class(mesh_t), pointer, dimension(:):: m
99  real, dimension(:,:,:,:) :: ff
100 END SUBROUTINE void
101 
102 !===============================================================================
103 !> Initialize a single 3-D wavenumber, with random amplitude and phase
104 !===============================================================================
105 SUBROUTINE single_solenoidal (self, m, ff, idx)
106  class(initial_t) :: self
107  type(index_t):: idx
108  class(mesh_t), pointer, dimension(:):: m
109  real, dimension(:,:,:,:) :: ff
110  !.............................................................................
111  type(random_t) :: random
112  integer :: ix, iy, iz
113  real(8) :: kx, ky, kz
114  !-----------------------------------------------------------------------------
115  call trace_begin('initial_t:single_solenoidal')
116  self%phi = 0.0
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))
129  if (self%mhd) then
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))
133  end if
134  end do
135  end do
136  end do
137  end associate
138 
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)
142  end associate
143  end if
144  call trace_end
145 END SUBROUTINE single_solenoidal
146 
147 !===============================================================================
148 !> Initialize a single 3-D wavenumber, with random amplitude and phase
149 !===============================================================================
150 SUBROUTINE single_compressive (self, m, ff, idx)
151  class(initial_t) :: self
152  type(index_t):: idx
153  class(mesh_t), pointer, dimension(:):: m
154  real, dimension(:,:,:,:) :: ff
155  !.............................................................................
156  type(random_t) :: random
157  integer :: ix, iy, iz
158  real(8) :: kx, ky, kz
159  !-----------------------------------------------------------------------------
160  call trace_begin('initial_t:single_compressive')
161  self%phi = 0.0
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)
174  if (self%mhd) then
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)
178  end if
179  end do
180  end do
181  end do
182  end associate
183 
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)
187  end associate
188  end if
189  call trace_end
190 END SUBROUTINE single_compressive
191 
192 !===============================================================================
193 !> Initialize a single 3-D wavenumber, with random amplitude and phase
194 !===============================================================================
195 SUBROUTINE single_advection (self, m, ff, idx)
196  class(initial_t) :: self
197  type(index_t):: idx
198  class(mesh_t), pointer, dimension(:):: m
199  real, dimension(:,:,:,:) :: ff
200  !.............................................................................
201  type(random_t) :: random
202  integer :: ix, iy, iz
203  real(8) :: kx, ky, kz
204  !-----------------------------------------------------------------------------
205  call trace_begin('initial_t:single_compressive')
206  self%phi = 0.0
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)
219  if (self%mhd) then
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)
223  end if
224  end do
225  end do
226  end do
227  end associate
228 
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)
232  end associate
233  end if
234  call trace_end
235 END SUBROUTINE single_advection
236 
237 !===============================================================================
238 !> Initialize from a raw data file
239 !===============================================================================
240 SUBROUTINE exa256 (self, m, ff, idx)
241  USE io_unit_mod
242  class(initial_t) :: self
243  type(index_t):: idx
244  class(mesh_t), pointer, dimension(:):: m
245  real, dimension(:,:,:,:) :: ff
246  real, dimension(:), allocatable:: u, d
247  !.............................................................................
248  integer :: ix, iy, iz, jy, jz, rec
249  integer(8) :: nx, ny, nz
250  !-----------------------------------------------------------------------------
251  call trace_begin('initial_t::exa256')
252  self%mhd = .false.
253  ff(:,:,:,:) = 0.0
254  allocate (u(m(1)%n))
255  allocate (d(m(1)%n))
256  open (io_unit%data, file='data/exa256/t0.02.stag', form='unformatted', &
257  access='direct', recl=m(1)%n*4)
258  self%time = 0.02
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
270  !print*,ix,iy,iz,rec
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
279  end do
280  end do
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))
285 
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)
289  end associate
290  end if
291  call trace_end
292 END SUBROUTINE exa256
293 
294 !===============================================================================
295 !> Initialize from a raw data file
296 !===============================================================================
297 SUBROUTINE raw (self, m, ff, idx)
299  class(initial_t):: self
300  type(index_t):: idx
301  class(mesh_t), pointer, dimension(:):: m
302  real, dimension(:,:,:,:):: ff
303  !.............................................................................
304  integer:: n(3), i(3), iv, iw, iostat
305  integer, save:: mv=4
306  logical, save:: first_time=.true.
307  character(len=64), save:: filename = 'data/exa256/t0.02.stag'
308  namelist /raw_params/ filename
309  !-----------------------------------------------------------------------------
310  call trace_begin('initial_t::raw')
311  !$omp critical (raw_cr)
312  self%mhd = .false.
313  self%time = 0.02
314  !-----------------------------------------------------------------------------
315  ! At first call, allocate a buffer and read the entire file -- this is OK since
316  ! we only use one thread per rank to initialize, computing also how many calls
317  ! will be made, so the buffer can be deallocated when all patches are done
318  !-----------------------------------------------------------------------------
319  if (first_time) then
320  first_time = .false.
321  rewind(io%input); read (io%input,raw_params,iostat=iostat)
322  if (io%master) write (*,raw_params)
323  n = nint(m%b/m%d)
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))
327  do iv=1,mv
328  read (io_unit%data, rec=iv) buffer(:,:,:,iv)
329  end do
330  close (io_unit%data)
331  end if
332  i = 1 + nint(m%p/m%d) - m%n/2 ! index offsets
333  ff = 0.0
334  do iw=1,mv
335  iv = merge(1,iw+1,iw==1)
336  ff(m(1)%li:m(1)%ui, &
337  m(2)%li:m(2)%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)
341  end do
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))
347  !$omp end critical (raw_cr)
348 
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)
352  end associate
353  end if
354  call trace_end
355 END SUBROUTINE raw
356 
357 !===============================================================================
358 SUBROUTINE cleanup
359  if (associated(buffer)) deallocate(buffer)
360 END SUBROUTINE cleanup
361 
362 !=======================================================================
363 !> Get entropy per unit volume from energy per unit volume
364 !> Note that energy and entropy are assumed to occupy the same memory
365 !=======================================================================
366 SUBROUTINE e2s(s, d, gamma)
367  real(8):: gamma
368  !.....................................................................
369  real, dimension(:,:,:):: s, d
370  real:: g1
371  !---------------------------------------------------------------------
372  g1 = gamma-1.0
373  s = d*log(s*g1/d**gamma)/g1
374 END SUBROUTINE e2s
375 
376 END MODULE
Template module for mesh.
Definition: mesh_mod.f90:4
Simple initical condition module for tests.
Definition: initial_mod.f90:7
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Definition: index_mod.f90:7
Definition: io_mod.f90:4
6th order stagger operators, with self-test procedure
Definition: stagger_2nd.f90:4