DISPATCH
initial_mod.f90
1 !===============================================================================
2 !> Simple initical condition module for tests
3 !>
4 !> Note that the isothermal and gamma parameters are communicated from the solver
5 !> via the hydro_parameters module.
6 !===============================================================================
7 MODULE initial_mod
8  USE const
9  USE io_mod
10  USE trace_mod
11  USE mesh_mod
12  USE hydro_parameters
13  USE omp_mod
14  USE index_mod
15  implicit none
16  private
17  real(8), parameter:: pi2=8d0*atan(1d0)
18  type, public:: initial_t
19  integer:: seed
20  real:: d0
21  logical:: mhd
22  real(8):: time=0d0
23  character(len=64):: solver
24  procedure(single_solenoidal), pointer:: condition=>null()
25  contains
26  procedure:: init
27  procedure, nopass:: cleanup
28  end type
29  real, save, pointer, dimension(:,:,:,:):: buffer => null()
30  real, save:: csound=1.
31 CONTAINS
32 
33 !===============================================================================
34 !> Setup a uniform initial state, with density=d0 and B=B0
35 !===============================================================================
36 SUBROUTINE init (self, solver, gamma1)
37  class(initial_t):: self
38  character(len=64):: solver
39  character(len=32), save:: type='void'
40  integer:: iostat
41  real, optional:: gamma1
42  logical, save:: first_time=.true.
43  namelist /initial_params/ type
44  !----------------------------------------------------------------------------
45  call trace_begin('initial_t%init')
46  !$omp critical (input_cr)
47  if (first_time) then
48  first_time = .false.
49  rewind(io%input)
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)
53  end if
54  !$omp end critical (input_cr)
55  self%solver = trim(solver)
56  select case (type)
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
63  case ('exa256')
64  self%condition => exa256
65  case ('raw')
66  self%condition => raw
67  case ('void')
68  self%condition => void
69  case default
70  print*,'UNKNOWN INITIAL CONDITION TYPE :'//trim(type)//':'
71  end select
72  call trace_end
73 END SUBROUTINE init
74 
75 !===============================================================================
76 SUBROUTINE cleanup
77  if (associated(buffer)) deallocate(buffer)
78 END SUBROUTINE cleanup
79 
80 !===============================================================================
81 !> Void initializer
82 !===============================================================================
83 SUBROUTINE void (self, m, ff, idx)
84  class(initial_t) :: self
85  type(index_t):: idx
86  class(mesh_t) :: m(3)
87  real, dimension(:,:,:,:) :: ff
88 END SUBROUTINE void
89 
90 !===============================================================================
91 !> Initialize a single 3-D wavenumber
92 !===============================================================================
93 SUBROUTINE single_solenoidal (self, m, ff, idx)
94  class(initial_t) :: self
95  class(mesh_t) :: m(3)
96  real, dimension(:,:,:,:) :: ff
97  type(index_t):: idx
98  !.............................................................................
99  integer :: ix, iy, iz
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
107  !-----------------------------------------------------------------------------
108  call trace_begin('initial_t:single_solenoidal', itimer=itimer)
109  if (first_time) then
110  !$omp critical (input_cr)
111  if (first_time) then
112  first_time = .false.
113  rewind(io%input)
114  read (io%input, initial_solenoidal_params)
115  mhd = mhd .and. any(b0 /= 0.0)
116  if (io%master) write (*, initial_solenoidal_params)
117  end if
118  !$omp end critical (input_cr)
119  end if
120  self%mhd = mhd
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))
132  if (self%mhd) then
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)
136  end if
137  end do
138  end do
139  end do
140  call energy (self, m, idx, ff)
141  end associate
142  call trace_end (itimer)
143 END SUBROUTINE single_solenoidal
144 
145 !===============================================================================
146 !> Initialize a single 3-D wavenumber, with non-random amplitude and zero phase
147 !===============================================================================
148 SUBROUTINE single_compressive (self, m, ff, idx)
149  class(initial_t) :: self
150  type(index_t):: idx
151  class(mesh_t) :: m(3)
152  real, dimension(:,:,:,:) :: ff
153  !.............................................................................
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
162  !-----------------------------------------------------------------------------
163  call trace_begin('initial_t:single_compressive', itimer=itimer)
164  if (first_time) then
165  !$omp critical (input_cr)
166  if (first_time) then
167  first_time = .false.
168  rewind(io%input)
169  read (io%input, initial_compressive_params)
170  mhd = mhd .and. any(b0 /= 0.0)
171  if (io%master) write (*, initial_compressive_params)
172  end if
173  !$omp end critical (input_cr)
174  end if
175  self%mhd = mhd
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)
187  if (self%mhd) then
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)
191  end if
192  end do
193  end do
194  end do
195  call energy (self, m, idx, ff)
196  end associate
197  call trace_end (itimer)
198 END SUBROUTINE single_compressive
199 
200 !===============================================================================
201 !> Initialize a single 3-D wavenumber, with non-random amplitude and phase
202 !===============================================================================
203 SUBROUTINE single_advection (self, m, ff, idx)
204  class(initial_t) :: self
205  type(index_t):: idx
206  class(mesh_t) :: m(3)
207  real, dimension(:,:,:,:) :: ff
208  !.............................................................................
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
216  !-----------------------------------------------------------------------------
217  call trace_begin('initial_t:single_advection', itimer=itimer)
218  if (first_time) then
219  !$omp critical (input_cr)
220  if (first_time) then
221  first_time = .false.
222  rewind(io%input)
223  read (io%input, initial_advection_params)
224  mhd = mhd .and. any(b0 /= 0.0)
225  if (io%master) write (*, initial_advection_params)
226  end if
227  !$omp end critical (input_cr)
228  end if
229  self%mhd = mhd
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)
241  if (self%mhd) then
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)
245  end if
246  end do
247  end do
248  end do
249  call energy (self, m, idx, ff)
250  end associate
251  call trace_end (itimer)
252 END SUBROUTINE single_advection
253 
254 !===============================================================================
255 !> Initialize from a raw data file
256 !===============================================================================
257 SUBROUTINE exa256 (self, m, ff, idx)
258  USE io_unit_mod
259  class(initial_t):: self
260  type(index_t):: idx
261  class(mesh_t):: m(3)
262  real, dimension(:,:,:,:):: ff
263  real, dimension(:), allocatable:: u, d
264  !.............................................................................
265  integer:: ix, iy, iz, jy, jz, rec, iv, li(3), ui(3)
266  integer(8):: nx, ny, nz
267  real(dp):: umax, dt
268  logical, save:: first_time=.true.
269  integer, save:: itimer=0
270  !-----------------------------------------------------------------------------
271  call trace_begin('initial_t::exa256', itimer=itimer)
272  self%mhd = .false.
273  ff(:,:,:,:) = 0.0
274  allocate (u(m(1)%n))
275  allocate (d(m(1)%n))
276  if (first_time) then
277  first_time = .false.
278  open (io_unit%data, file='data/exa256/t0.02.raw', form='unformatted', &
279  access='direct', recl=m(1)%n*4)
280  end if
281  self%time = 0.02
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
288  li = m%lb
289  ui = m%ub
290  umax = 0.0
291  do jz=li(3),ui(3)
292  do jy=li(2),ui(2)
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))))
306  end do
307  end do
308  do iv=idx%px,idx%pz
309  ff(:,:,:,iv) = ff(:,:,:,idx%d)*ff(:,:,:,iv)
310  end do
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
320 
321 !===============================================================================
322 !> Initialize from a raw data file, assumed to contain isothermal data, in the
323 !> order rho, Ux, Uy, Uz
324 !===============================================================================
325 SUBROUTINE raw (self, m, ff, idx)
327  class(initial_t):: self
328  type(index_t):: idx
329  class(mesh_t):: m(3)
330  real, dimension(:,:,:,:):: ff
331  !.............................................................................
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
336  !-----------------------------------------------------------------------------
337  call trace_begin('initial_t::raw', itimer=itimer)
338  !$omp critical (raw_cr)
339  self%mhd = .false.
340  self%time = 0.02
341  !-----------------------------------------------------------------------------
342  ! At first call, allocate a buffer and read the entire file -- this is OK since
343  ! we only use one thread per rank to initialize, computing also how many calls
344  ! will be made, so the buffer can be deallocated when all patches are done
345  !-----------------------------------------------------------------------------
346  if (first_time) then
347  first_time = .false.
348  n = nint(m%b/m%d)
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))
352  do iv=1,mv
353  read (io_unit%data, rec=iv) buffer(:,:,:,iv)
354  end do
355  close (io_unit%data)
356  end if
357  i = 1 + nint(m%p/m%d) - m%n/2 ! index offsets
358  jv(1) = idx%d ! variable index translation
359  jv(2) = idx%px
360  jv(3) = idx%py
361  jv(4) = idx%pz
362  ff = 0.0
363  do iv=1,mv
364  ff(m(1)%li:m(1)%ui, &
365  m(2)%li:m(2)%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)
369  end do
370  do iv=idx%px,idx%pz
371  ff(:,:,:,iv) = ff(:,:,:,idx%d)*ff(:,:,:,iv)
372  end do
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))
378  !$omp end critical (raw_cr)
379  call trace_end (itimer)
380 END SUBROUTINE raw
381 
382 !===============================================================================
383 !> Calculate total energy
384 !===============================================================================
385 SUBROUTINE energy (self, m, idx, ff)
386  class(initial_t) :: self
387  class(mesh_t):: m(3)
388  type(index_t):: idx
389  real, dimension(:,:,:,:):: ff
390  !.............................................................................
391  real(dp):: fgamma
392  !-----------------------------------------------------------------------------
393  call trace%begin ('initial_t%energy')
394  if (isothermal .or. gamma==1.0_dp) then
395  fgamma = 1.0_dp
396  else
397  fgamma = gamma*(gamma-1.0_dp)
398  end if
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
403  if (self%mhd) then
404  ff(:,:,:,idx%e) = ff(:,:,:,idx%e) + &
405  0.5*(ff(:,:,:,idx%bx)**2 + &
406  ff(:,:,:,idx%by)**2 + &
407  ff(:,:,:,idx%bz)**2)
408  end if
409  !print *, 'IIII', idx%e, fgamma, &
410  !minval(ff(m(1)%lb:m(1)%ub,m(2)%lb:m(2)%ub,m(3)%lb:m(3)%ub,idx%d)), &
411  !maxval(ff(m(1)%lb:m(1)%ub,m(2)%lb:m(2)%ub,m(3)%lb:m(3)%ub,idx%d)), &
412  !minval(ff(m(1)%lb:m(1)%ub,m(2)%lb:m(2)%ub,m(3)%lb:m(3)%ub,idx%e)), &
413  !maxval(ff(m(1)%lb:m(1)%ub,m(2)%lb:m(2)%ub,m(3)%lb:m(3)%ub,idx%e))
414  call trace%end
415 END SUBROUTINE energy
416 
417 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