DISPATCH
hd_mod.f90
1 !===============================================================================
2 !> RAMSES Godunov solvers
3 !>
4 !> godunov_fine
5 !> unsplit 1 2 3 ui uo ub
6 !> ctoprim +---+---+......+---+---+
7 !> uslope:dq +---+......+---+
8 !> trace3d:qp +---+......+---+
9 !> trace3d:qm +......+---+---+
10 !> cmpflxm +......+---+
11 !> riemann%hllc +......+
12 !>
13 !> Note that it is NECESSARY to have a call to courant_condition here, after
14 !> ctoprim (which calculates an up-to-date self%u_max), and before trace3d,
15 !> which needs the new self%dtime
16 !===============================================================================
17 MODULE hd_mod
18  !.............................................................................
19  USE const
20  USE io_mod
21  USE io_unit_mod
22  USE extras_mod
23  USE link_mod
24  USE riemann_mod
25  USE hydro_parameters
26  USE trace_mod
27  USE mpi_mod
28  USE omp_mod
29  USE omp_lock_mod
30  USE eos_mod
31  USE index_mod
32  USE patch_mod
33  USE validate_mod
34  implicit none
35  PRIVATE
36  type, public, extends(extras_t):: hd_t
37  logical:: mhd=.false.
38  integer:: idim=3
39  character(len=10):: riemann='hllc'
40  real, pointer, dimension(:,:,:,:):: uin ! conserved
41  real, pointer, dimension(:,:,:,:):: unew ! new conserved
42  real, pointer, dimension(:,:,:,:):: q ! primitive
43  real, pointer, dimension(:,:,:,:,:):: qp, qm ! face centered
44  real, pointer, dimension(:,:,:,:,:):: dq ! slopes
45  real, pointer, dimension(:,:,:,:,:):: fl ! 3D fluxes
46  real, pointer, dimension(:,:,:):: c, s ! scratch
47  contains
48  procedure:: init
49  procedure:: dealloc
50  procedure:: update
51  procedure:: gas_pressure
52  end type
53  integer, save:: courant_type=0
54  logical, save:: debug=.false.
55  logical, save:: first_time=.true.
56  logical, save:: detailed_timer=.false.
57  logical, save:: unsigned=.true.
58  character(len=16):: eqn='mhd'
59 CONTAINS
60 
61 !===============================================================================
62 !> Initialize the HD solver
63 !===============================================================================
64 SUBROUTINE init (self)
65  class(hd_t):: self
66  character(len=10):: riemann='hllc'
67  real, save:: csound=1.
68  integer:: i, iv
69  namelist /ramses_params/ gamma, riemann, slope_type, courant_factor, smallr, &
70  smallc, isothermal, csound, detailed_timer, unsigned, courant_type, debug
71  type(index_t):: idx
72  character(len=120):: id = &
73  '$Id: 60e649508a40479a148c99021f42dfa3fe535b77 $ solvers/ramses/hydro/hd_mod.f90'
74 !...............................................................................
75  call trace%begin('hd_mod%init')
76  call trace%print_id (id)
77  self%mhd = .false.
78  self%nw = 1
79  self%ng = 2
80  self%kind = 'ramses_hd_patch'
81  if (self%nv == 0) self%nv = 5
82  call self%idx%init (5, self%mhd)
83  call self%patch_t%init
84  !$omp critical (input_cr)
85  if (first_time) then
86  first_time = .false.
87  rewind(io_unit%input)
88  read (io_unit%input, ramses_params)
89  if (mpi%rank==0) write (*, ramses_params)
90  end if
91  if (gamma == 1d0) &
92  isothermal = .true.
93  if (isothermal) &
94  gamma = 1d0
95  !$omp end critical (input_cr)
96  self%riemann = riemann
97  self%courant = courant_factor
98  self%gamma = gamma
99  self%csound = csound
100  self%staggered = .false.
101  do i=1,3
102  self%mesh(i)%h = 0d0
103  end do
104  self%unsigned(self%idx%d) = unsigned
105  self%unsigned(self%idx%e) = unsigned
106  self%pervolume(self%idx%px) = unsigned
107  self%pervolume(self%idx%py) = unsigned
108  self%pervolume(self%idx%pz) = unsigned
109  call self%gpatch_t%init
110  call self%extras_t%init
111  call trace%end
112 END SUBROUTINE init
113 
114 !===============================================================================
115 !> Deallocate permanently allocated arrays
116 !===============================================================================
117 SUBROUTINE dealloc (self)
118  class(hd_t):: self
119  !-----------------------------------------------------------------------------
120  call trace%begin ('hd_t%dealloc')
121  call self%extras_t%dealloc
122  call trace%end
123 END SUBROUTINE dealloc
124 
125 !===============================================================================
126 !===============================================================================
127 SUBROUTINE update (self)
128  class(hd_t):: self
129  integer:: n(3), m, i, l(3), u(3)
130  integer, save:: itimer=0
131  real, dimension(:,:,:,:), pointer:: p
132  real, dimension(:,:,:), pointer:: d, e, ux=>null(), uy=>null(), uz=>null()
133  real, dimension(:,:,:), pointer:: tmp
134  real:: error
135  !-----------------------------------------------------------------------------
136  call trace%begin ('hd_t%update', itimer=itimer)
137  if (io%smallr /= 0.0) smallr = io%smallr
138  if (io%courant /= 0.0) self%courant = io%courant
139  !-----------------------------------------------------------------------------
140  ! The call to output below catches the state after both guard zone downloads
141  ! and boundary conditions. The call here is a fall back; if another call is
142  ! placed at an earlier point in the task handling the call here will have no
143  ! effect, since out_next has already been updated.
144  !-----------------------------------------------------------------------------
145  call self%output
146  !-----------------------------------------------------------------------------
147  ! Memory mapping and allocations
148  !-----------------------------------------------------------------------------
149  self%uin => self%mem(:,:,:,:,self%it,1)
150  self%unew => self%mem(:,:,:,:,self%new,1)
151  n = self%mesh%gn
152  m = n(1)
153  allocate (self%q (n(1),n(2),n(3),self%nv))
154  allocate (self%qp(n(1),n(2),n(3),self%nv,ndim))
155  allocate (self%qm(n(1),n(2),n(3),self%nv,ndim))
156  allocate (self%dq(n(1),n(2),n(3),self%nv,ndim))
157  allocate (self%fl(n(1),n(2),n(3),self%nv,ndim))
158  allocate (self%c(n(1),n(2),n(3)))
159  allocate (self%s(n(1),n(2),n(3)))
160  !self%q=0.0; self%qp=0.0; self%qm=0.0; self%dq=0.0; self%fl=0.; self%c=0.0; self%s=0.0
161  if (detailed_timer) call trace_end (itimer)
162  if (debug) then
163  if (self%id==229.or.self%id==231) then
164  print *,self%id, 'mk1:px', self%mem(19-2:19+2,19,19,self%idx%px,self%it,1), self%idx%px
165  print *,self%id, 'mk1: d', self%mem(19-2:19+2,19,19,self%idx%d ,self%it,1)
166  print *,self%id, 'mk1: e', self%mem(19-2:19+2,19,19,self%idx%e ,self%it,1), self%idx%e
167  end if
168  if (self%id==230) then
169  print *,self%id, 'mk1:px', self%mem( m-7:m-3 ,19,19,self%idx%px,self%it,1), self%idx%px
170  print *,self%id, 'mk1: d', self%mem( m-7:m-3 ,19,19,self%idx%d ,self%it,1)
171  print *,self%id, 'mk1: e', self%mem( m-7:m-3 ,19,19,self%idx%e ,self%it,1), self%idx%e
172  end if
173  end if
174  d => self%mem(:,:,:, self%idx%d ,self%it,1)
175  e => self%mem(:,:,:, self%idx%e ,self%it,1)
176  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%it,1)
177  !-----------------------------------------------------------------------------
178  ! Godunov solver
179  !-----------------------------------------------------------------------------
180  call godunov (self)
181  if (isothermal) &
182  self%mem(:,:,:,self%idx%e,self%new,1) = &
183  self%mem(:,:,:,self%idx%d,self%new,1)*self%csound**2
184  if (detailed_timer) call trace_begin ('hd_t%update', itimer=itimer)
185  !-----------------------------------------------------------------------------
186  ! Add the external force per unit mass as a source term
187  !-----------------------------------------------------------------------------
188  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
189  if (allocated(self%force_per_unit_mass)) then
190  do i=1,3
191  p(:,:,:,i) = p(:,:,:,i) + self%dtime*self%force_per_unit_mass(:,:,:,i)*d
192  end do
193  end if
194  if (allocated(self%force_per_unit_volume)) then
195  p = p + self%dtime*self%force_per_unit_volume
196  end if
197  call validate%check (self%link, p, 'after force')
198  if (omp_lock%tasks) then
199  call self%lock%unset ('update')
200  end if
201  !-----------------------------------------------------------------------------
202  ! Deallocate
203  !-----------------------------------------------------------------------------
204  call self%counter_update
205  deallocate (self%q,self%qp,self%qm,self%dq,self%fl,self%c,self%s)
206  call trace%end (itimer)
207 END SUBROUTINE update
208 
209 !=======================================================================
210 !> Get velocities from momentum
211 !=======================================================================
212 SUBROUTINE p2u(self, U, it)
213  class(hd_t) :: self
214  real, dimension(:,:,:,:), pointer:: u
215  integer:: it
216  !.....................................................................
217  integer:: i
218  !---------------------------------------------------------------------
219  associate(d => self%mem(:,:,:,self%idx%d, it,1), &
220  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
221  do i=1,size(p,4)
222  u(:,:,:,i) = p(:,:,:,i)/d
223  end do
224  end associate
225 END SUBROUTINE p2u
226 
227 !=======================================================================
228 !> Put velocities to momentum
229 !=======================================================================
230 SUBROUTINE u2p(self, U, it)
231  class(hd_t) :: self
232  real, dimension(:,:,:,:), pointer:: u
233  integer:: it
234  !.....................................................................
235  integer:: i
236  !---------------------------------------------------------------------
237  associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
238  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
239  do i=1,size(p,4)
240  p(:,:,:,i)= u(:,:,:,i)*d
241  end do
242  end associate
243 END SUBROUTINE u2p
244 
245 !=======================================================================
246 !> Get thermal energy from total/specific energy
247 !=======================================================================
248 SUBROUTINE e2e_th(self, E_th, it)
249  class(hd_t) :: self
250  real, dimension(:,:,:), pointer:: e_th
251  integer:: it
252  !.....................................................................
253  real, dimension(:,:,:), pointer:: k
254  integer :: i
255  !---------------------------------------------------------------------
256  associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
257  e => self%mem(:,:,:,self%idx%e ,it,1), &
258  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
259  allocate(k(size(d,1),size(d,2),size(d,3)))
260  k = 0.0
261  do i=1,size(p,4)
262  k = k+p(:,:,:,i)**2
263  end do
264  k = 0.5*k/d
265  e_th = (e-k)/d
266  deallocate(k)
267  end associate
268 END SUBROUTINE e2e_th
269 
270 !=======================================================================
271 !> Get etropy per unit mass from total energy
272 !=======================================================================
273 SUBROUTINE e2s(self, s, it)
274  class(hd_t) :: self
275  real, dimension(:,:,:), pointer:: s
276  integer:: it
277  !.....................................................................
278  real, dimension(:,:,:), pointer:: k
279  integer :: i
280  !---------------------------------------------------------------------
281  if (isothermal) then
282  s = 0.0
283  return
284  end if
285  associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
286  e => self%mem(:,:,:,self%idx%e ,it,1), &
287  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
288  allocate(k(size(d,1),size(d,2),size(d,3)))
289  k = 0.0
290  do i=1,size(p,4)
291  k = k+p(:,:,:,i)**2
292  end do
293  k = 0.5*k/d
294  s = log((e-k)*(gamma-1d0)/d**gamma)/(gamma-1d0)
295  deallocate(k)
296  end associate
297 END SUBROUTINE e2s
298 
299 !=======================================================================
300 !> Put etropy per unit mass to total energy
301 !=======================================================================
302 SUBROUTINE s2e (self, s, it)
303  class(hd_t) :: self
304  real, dimension(:,:,:), pointer:: s
305  integer:: it
306  !.....................................................................
307  real, dimension(:,:,:), pointer:: k
308  integer :: i
309  !---------------------------------------------------------------------
310  associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
311  e => self%mem(:,:,:,self%idx%e ,it,1), &
312  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
313  if (isothermal) then
314  e = d*self%csound**2
315  return
316  end if
317  allocate(k(size(d,1),size(d,2),size(d,3)))
318  k = 0.0
319  do i=1,size(p,4)
320  k = k+p(:,:,:,i)**2
321  end do
322  k = 0.5*k/d
323  e = d**gamma*exp(s*(gamma-1d0)) + k
324  deallocate(k)
325  end associate
326 END SUBROUTINE s2e
327 
328 !=======================================================================
329 !> Put thermal energy to total/specific energy
330 !=======================================================================
331 SUBROUTINE e_th2e(self, E_th, it)
332  class(hd_t) :: self
333  integer:: it
334  real, dimension(:,:,:), pointer:: e_th
335  !.....................................................................
336  real, dimension(:,:,:), pointer:: k
337  integer :: i
338  !---------------------------------------------------------------------
339  associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
340  e => self%mem(:,:,:,self%idx%e ,it,1), &
341  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
342  allocate(k(size(d,1),size(d,2),size(d,3)))
343  k = 0.0
344  do i=1,size(p,4)
345  k = k+p(:,:,:,i)**2
346  end do
347  k = 0.5*k/d
348  e = e_th*d+k
349  deallocate(k)
350  end associate
351 END SUBROUTINE e_th2e
352 
353 !===============================================================================
354 subroutine godunov (self)
355  class(hd_t)::self
356  !-----------------------------------------------------------------------------
357  ! This routine accesses the conserved variables in self. It then calls
358  ! the Godunov solver that computes fluxes. These fluxes are zeroed at
359  ! coarse-fine boundaries, since contribution from finer levels has
360  ! already been taken into account. Conservative variables are updated
361  ! and stored in array unew(:), both at the current level and at the
362  ! coarser level if necessary.
363  !-----------------------------------------------------------------------------
364  integer::i0,j0,k0,i,j,k,idim,ivar,li(3),ui(3)
365  integer,save:: itimer=0
366  associate(unew=>self%unew, uin=>self%uin, fl=>self%fl)
367  !-----------------------------------------------------------------------------
368  ! Compute flux using second-order Godunov method
369  !-----------------------------------------------------------------------------
370  if (.not.detailed_timer) call trace%begin ('hd_t%godunov')
371  li = self%mesh%li
372  ui = self%mesh%ui
373  call check_small (self, uin)
374  call unsplit (self)
375  !-----------------------------------------------------------------------------
376  ! Conservative update, looping over internal cells only
377  !-----------------------------------------------------------------------------
378  if (detailed_timer) call trace%begin ('hd_t%godunov', itimer=itimer)
379  if (omp_lock%tasks) then
380  call self%lock%set ('update')
381  end if
382  unew = uin
383  do ivar=1,self%nv
384  do idim=1,self%ndim
385  i0 = merge(1,0,idim==1)
386  j0 = merge(1,0,idim==2)
387  k0 = merge(1,0,idim==3)
388  do k=li(3),ui(3)
389  do j=li(2),ui(2)
390  do i=li(1),ui(1)
391  ! Update conservative variables
392  unew(i,j,k,ivar) = unew(i,j,k,ivar)+ &
393  & (fl(i ,j ,k ,ivar,idim) &
394  & -fl(i+i0,j+j0,k+k0,ivar,idim))
395  end do
396  end do
397  end do
398  end do
399  end do
400  if (detailed_timer) then
401  call trace%end (itimer)
402  else
403  call trace%end
404  end if
405  end associate
406 end subroutine godunov
407 
408 !===============================================================================
409 !===============================================================================
410 subroutine check_small (self, uu)
411  class(hd_t):: self
412  real, dimension(:,:,:,:):: uu
413  integer:: i,j,k,lo(3),uo(3),n_smalld,n_smallp
414  real(dp):: c,d,u,v,w,p,ekin
415  integer, save:: nprint=10
416  integer, save:: itimer=0
417  integer:: it, jt, iw
418  logical:: panic
419  !-----------------------------------------------------------------------------
420  if (detailed_timer) call trace%begin ('hd_t%check_small', itimer=itimer)
421  n_smalld = 0
422  n_smallp = 0
423  panic = .false.
424  lo = self%mesh%lb
425  uo = self%mesh%ub
426  do k=lo(3),uo(3)
427  do j=lo(2),uo(2)
428  do i=lo(1),uo(1)
429  if (uu(i,j,k,self%idx%d) < smallr) then
430  panic = .true.
431  if (nprint>0) then
432  nprint = nprint-1
433  write (io_unit%log,'(i7,3i5,1p,2e12.3)') self%id,i,j,k,uu(i,j,k,self%idx%d),smallr
434  end if
435  n_smalld = n_smalld+1
436  uu(i,j,k,self%idx%d) = smallr
437  end if
438  d = uu(i,j,k,self%idx%d)
439  c = 1./d
440  u = uu(i,j,k,self%idx%px)*c
441  v = uu(i,j,k,self%idx%py)*c
442  w = uu(i,j,k,self%idx%pz)*c
443  ekin = 0.5*d*(u**2+v**2+w**2)
444  if (isothermal) then
445  uu(i,j,k,self%idx%e) = max(uu(i,j,k,self%idx%e), smallr*smallc**2)
446  else
447  p = (uu(i,j,k,self%idx%e)-ekin)*(gamma-one)
448  if (p < smallr*smallc**2) then
449  p = smallr*smallc**2
450  uu(i,j,k,self%idx%e) = p/(gamma-one)+ekin
451  n_smallp = n_smallp+1
452  end if
453  end if
454  end do
455  end do
456  end do
457  if (n_smalld+n_smallp > 0) then
458  write (io_unit%log,*) 'small: id,nr,np', self%id,n_smalld,n_smallp
459  end if
460  if (panic .and. self%n_dump > 0) then
461  !$omp critical (panic_cr)
462  self%n_dump = self%n_dump-1
463  write (io_unit%dump) self%id, self%gn, self%nv, self%nt, self%nw, self%istep
464  write (io_unit%dump) self%time, self%dtime
465  do iw=1,self%nw
466  do it=1,self%nt
467  !$omp atomic read
468  jt = self%iit(it)
469  write (io_unit%dump) self%mem(:,:,:,:,jt,iw)
470  end do
471  end do
472  write (io_unit%mpi,*) 'density < SMALLR, dumped', self%id, self%time
473  !$omp end critical (panic_cr)
474  end if
475  if (detailed_timer) call trace%end (itimer)
476 end subroutine check_small
477 
478 !===============================================================================
479 !> UNSPLIT Unsplit second order Godunov integrator for
480 !> polytropic gas dynamics using MUSCL-HANCOCK
481 !===============================================================================
482 subroutine unsplit(self)
483  class(hd_t):: self
484  real(dp):: dtds(3)
485  !-----------------------------------------------------------------------------
486  call trace%begin ('hd_t%unsplit')
487  call ctoprim(self) ! Primitive variables, sound speeds
488  call self%courant_condition (detailed_timer)
489  call uslope(self) ! TVD slopes
490  call trace3d(self) ! predictor step
491  call cmpflxm(self) ! Solve for 1D flux in X direction
492  call trace%end
493 end subroutine unsplit
494 
495 !===============================================================================
496 !> Conservative (uin) to primitive (q) variables, stored in the order
497 !> 1:d, 2:u, 3:v, 4:w, 5:p
498 !===============================================================================
499 subroutine ctoprim (self)
500  class(hd_t):: self
501  real(dp):: smallp, dtxhalf
502  integer,save:: itimer=0
503  integer:: l(3), u(3)
504  real, dimension(:,:,:), pointer:: d, e, pg
505  associate(uin=>self%uin, q=>self%q, c=>self%c, s=>self%s)
506  !-----------------------------------------------------------------------------
507  if (detailed_timer) then
508  call trace%begin ('hd_mod::ctoprim', itimer=itimer)
509  else
510  call trace%begin ('hd_mod::ctoprim')
511  end if
512  smallp = smallr*smallc**2/gamma
513  q(:,:,:,1) = max(uin(:,:,:,self%idx%d),smallr)
514  s = 1.d0/q(:,:,:,1)
515  q(:,:,:,2) = uin(:,:,:,self%idx%px)*s
516  q(:,:,:,3) = uin(:,:,:,self%idx%py)*s
517  q(:,:,:,4) = uin(:,:,:,self%idx%pz)*s
518  s = q(:,:,:,2)**2+q(:,:,:,3)**2+q(:,:,:,4)**2
519  !-----------------------------------------------------------------------------
520  ! Equation-of-state: s is the internal energy, and pg is the gas pressure
521  !-----------------------------------------------------------------------------
522  s = max((uin(:,:,:,self%idx%e)-half*q(:,:,:,1)*s(:,:,:)),smallp)
523  d => q(:,:,:,1)
524  e => s
525  pg => q(:,:,:,5)
526  if (isothermal) then
527  pg = d*self%csound**2
528  else
529  pg = e*(gamma-1d0)
530  end if
531  !call eos%lookup (shape(e), ee=e, d=d, pg=pg, gamma=gamma)
532  !
533  if (self%do_pebbles) then
534  self%vgas(:,:,:,1:3,self%it) = q(:,:,:,2:4)
535  self%pressure(:,:,:, self%it) = pg
536  end if
537  c(:,:,:)=sqrt(gamma*q(:,:,:,5)/q(:,:,:,1))
538  if (courant_type==1) then
539  c = c + max(abs(q(:,:,:,2)), abs(q(:,:,:,3)), abs(q(:,:,:,4)))
540  else
541  c = c + (abs(q(:,:,:,2)) + abs(q(:,:,:,3)) + abs(q(:,:,:,4)))/3.0
542  end if
543  l = self%mesh%lo
544  u = self%mesh%uo
545  self%u_max = maxval(c(l(1):u(1),l(2):u(2),l(3):u(3)))
546  if (detailed_timer) then
547  call trace%end (itimer)
548  else
549  call trace%end
550  end if
551  end associate
552 end subroutine ctoprim
553 
554 !===============================================================================
555 !> Slope limiters
556 !===============================================================================
557 subroutine uslope(self)
558  class(hd_t):: self
559  integer::ngrid
560  integer::i, j, k, n, l(3), u(3), islope_type
561  real(dp)::dsgn, dlim, dcen, dlft, drgt, slop, c1, c2
562  real(dp)::dfll,dflm,dflr,dfml,dfmm,dfmr,dfrl,dfrm,dfrr
563  real(dp)::dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl
564  real(dp)::dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm
565  real(dp)::dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr
566  real(dp)::vmin,vmax,dfx,dfy,dfz,dff,idff,xslope_type
567  integer,save:: itimer=0
568  associate(q=>self%q, dq=>self%dq)
569  !-----------------------------------------------------------------------------
570  if (detailed_timer) then
571  call trace%begin ('hd_mod::uslope', itimer=itimer)
572  else
573  call trace%begin ('hd_mod::uslope')
574  end if
575  !
576  islope_type = int(slope_type)
577  l = self%mesh%lo
578  u = self%mesh%uo
579  if(islope_type==-1)then
580  do n = 1, self%nv
581  do k = l(3),u(3)
582  do j = l(2),u(2)
583  do i = l(1),u(1)
584  dq(i,j,k,n,1) = half*(q(i+1,j,k,n) - q(i-1,j,k,n))
585  dq(i,j,k,n,2) = half*(q(i,j+1,k,n) - q(i,j-1,k,n))
586  dq(i,j,k,n,3) = half*(q(i,j,k+1,n) - q(i,j,k-1,n))
587  end do
588  end do
589  end do
590  end do
591  else if(islope_type==0)then
592  dq=0.0_dp
593  if (detailed_timer) then
594  call trace%end (itimer)
595  else
596  call trace%end
597  end if
598  return
599  else if (islope_type==1) then ! minmod
600  do n = 1, self%nv
601  do k = l(3),u(3)
602  do j = l(2),u(2)
603  do i = l(1),u(1)
604  ! slopes in first coordinate direction
605  dlft = q(i ,j,k,n) - q(i-1,j,k,n)
606  drgt = q(i+1,j,k,n) - q(i ,j,k,n)
607  dq(i,j,k,n,1) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
608  ! slopes in second coordinate direction
609  dlft = q(i,j ,k,n) - q(i,j-1,k,n)
610  drgt = q(i,j+1,k,n) - q(i,j ,k,n)
611  dq(i,j,k,n,2) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
612  ! slopes in third coordinate direction
613  dlft = q(i,j,k ,n) - q(i,j,k-1,n)
614  drgt = q(i,j,k+1,n) - q(i,j,k ,n)
615  dq(i,j,k,n,3) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
616  end do
617  end do
618  end do
619  end do
620  else if (islope_type==2) then ! moncen
621  c1 = slope_type
622  c2 = 0.5_dp/c1
623  do n = 1, self%nv
624  do k = l(3),u(3)
625  do j = l(2),u(2)
626  do i = l(1),u(1)
627  ! slopes in first coordinate direction
628  dlft = c1*(q(i ,j,k,n) - q(i-1,j,k,n))
629  drgt = c1*(q(i+1,j,k,n) - q(i ,j,k,n))
630  dcen = c2*(dlft+drgt)
631  dsgn = sign(one, dcen)
632  slop = min(abs(dlft),abs(drgt))
633  dlim = merge(zero,slop,(dlft*drgt)<=zero)
634  dq(i,j,k,n,1) = dsgn*min(dlim,abs(dcen))
635 ! end do
636 ! do i = l(1),u(1)
637  ! slopes in second coordinate direction
638  dlft = c1*(q(i,j ,k,n) - q(i,j-1,k,n))
639  drgt = c1*(q(i,j+1,k,n) - q(i,j ,k,n))
640  dcen = c2*(dlft+drgt)
641  dsgn = sign(one,dcen)
642  slop = min(abs(dlft),abs(drgt))
643  dlim = merge(zero,slop,(dlft*drgt)<=zero)
644  dq(i,j,k,n,2) = dsgn*min(dlim,abs(dcen))
645 ! end do
646 ! do i = l(1),u(1)
647  ! slopes in third coordinate direction
648  dlft = c1*(q(i,j,k ,n) - q(i,j,k-1,n))
649  drgt = c1*(q(i,j,k+1,n) - q(i,j,k ,n))
650  dcen = c2*(dlft+drgt)
651  dsgn = sign(one,dcen)
652  slop = min(abs(dlft),abs(drgt))
653  dlim = merge(zero,slop,(dlft*drgt)<=zero)
654  dq(i,j,k,n,3) = dsgn*min(dlim,abs(dcen))
655  end do
656  end do
657  end do
658  end do
659  else if (islope_type==3) then ! positivity preserving 3d unsplit slope
660  xslope_type = real(slope_type - 2_dp,kind=dp)
661  do n = 1, self%nv
662  do k = l(3),u(3)
663  do j = l(2),u(2)
664  do i = l(1),u(1)
665  dflll = q(i-1,j-1,k-1,n)-q(i,j,k,n)
666  dflml = q(i-1,j ,k-1,n)-q(i,j,k,n)
667  dflrl = q(i-1,j+1,k-1,n)-q(i,j,k,n)
668  dfmll = q(i ,j-1,k-1,n)-q(i,j,k,n)
669  dfmml = q(i ,j ,k-1,n)-q(i,j,k,n)
670  dfmrl = q(i ,j+1,k-1,n)-q(i,j,k,n)
671  dfrll = q(i+1,j-1,k-1,n)-q(i,j,k,n)
672  dfrml = q(i+1,j ,k-1,n)-q(i,j,k,n)
673  dfrrl = q(i+1,j+1,k-1,n)-q(i,j,k,n)
674  !
675  dfllm = q(i-1,j-1,k ,n)-q(i,j,k,n)
676  dflmm = q(i-1,j ,k ,n)-q(i,j,k,n)
677  dflrm = q(i-1,j+1,k ,n)-q(i,j,k,n)
678  dfmlm = q(i ,j-1,k ,n)-q(i,j,k,n)
679  dfmmm = q(i ,j ,k ,n)-q(i,j,k,n)
680  dfmrm = q(i ,j+1,k ,n)-q(i,j,k,n)
681  dfrlm = q(i+1,j-1,k ,n)-q(i,j,k,n)
682  dfrmm = q(i+1,j ,k ,n)-q(i,j,k,n)
683  dfrrm = q(i+1,j+1,k ,n)-q(i,j,k,n)
684  !
685  dfllr = q(i-1,j-1,k+1,n)-q(i,j,k,n)
686  dflmr = q(i-1,j ,k+1,n)-q(i,j,k,n)
687  dflrr = q(i-1,j+1,k+1,n)-q(i,j,k,n)
688  dfmlr = q(i ,j-1,k+1,n)-q(i,j,k,n)
689  dfmmr = q(i ,j ,k+1,n)-q(i,j,k,n)
690  dfmrr = q(i ,j+1,k+1,n)-q(i,j,k,n)
691  dfrlr = q(i+1,j-1,k+1,n)-q(i,j,k,n)
692  dfrmr = q(i+1,j ,k+1,n)-q(i,j,k,n)
693  dfrrr = q(i+1,j+1,k+1,n)-q(i,j,k,n)
694  !
695  vmin = min(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
696  & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
697  & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
698  vmax = max(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
699  & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
700  & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
701  !
702  dfx = half*(q(i+1,j,k,n)-q(i-1,j,k,n))
703  dfy = half*(q(i,j+1,k,n)-q(i,j-1,k,n))
704  dfz = half*(q(i,j,k+1,n)-q(i,j,k-1,n))
705  idff = xslope_type / sqrt(dfx**2+dfy**2+dfz**2+tiny(1.0_dp)) ! xslope_type = 1 or 2 dep. on slope_type
706  !
707  ! vmin and vmax are the lower and upper limit of one-cell differences
708  ! so a constant slope with step dq will give vmin=-dq and vmax=+dq,
709  ! while dfx (dff for any direction) will also be = dq, unless we keep
710  ! the 'half' in front of the sqrt. With that kept, the dlim factor
711  ! will remain unity until the slopes differ with a factor of 2. Then,
712  ! if any of the 26 slopes is steeper than twice the gradient slope,
713  ! the gradient will be diminished.
714  !
715  ! idff > 0, so 1/dff = idff < 0.001 * huge is equivalent to dff > 0
716  slop = merge(min(one,min(abs(vmin),abs(vmax))*idff), one, idff<0.001_dp * huge(1.0_dp))
717  dlim = slop
718  dq(i,j,k,n,1) = dlim*dfx
719  dq(i,j,k,n,2) = dlim*dfy
720  dq(i,j,k,n,3) = dlim*dfz
721  end do
722  end do
723  end do
724  end do
725  else if (islope_type==4) then ! positivity preserving 3d unsplit slope
726  xslope_type = real(slope_type - 3_dp,kind=dp)
727  do k = l(3),u(3)
728  do j = l(2),u(2)
729  do i = l(1),u(1)
730  block
731  real(dp):: adfx(self%nv), adfy(self%nv), adfz(self%nv)
732  dlim = 1.0
733  do n = 1, self%nv
734  dflll = q(i-1,j-1,k-1,n)-q(i,j,k,n)
735  dflml = q(i-1,j ,k-1,n)-q(i,j,k,n)
736  dflrl = q(i-1,j+1,k-1,n)-q(i,j,k,n)
737  dfmll = q(i ,j-1,k-1,n)-q(i,j,k,n)
738  dfmml = q(i ,j ,k-1,n)-q(i,j,k,n)
739  dfmrl = q(i ,j+1,k-1,n)-q(i,j,k,n)
740  dfrll = q(i+1,j-1,k-1,n)-q(i,j,k,n)
741  dfrml = q(i+1,j ,k-1,n)-q(i,j,k,n)
742  dfrrl = q(i+1,j+1,k-1,n)-q(i,j,k,n)
743  !
744  dfllm = q(i-1,j-1,k ,n)-q(i,j,k,n)
745  dflmm = q(i-1,j ,k ,n)-q(i,j,k,n)
746  dflrm = q(i-1,j+1,k ,n)-q(i,j,k,n)
747  dfmlm = q(i ,j-1,k ,n)-q(i,j,k,n)
748  dfmmm = q(i ,j ,k ,n)-q(i,j,k,n)
749  dfmrm = q(i ,j+1,k ,n)-q(i,j,k,n)
750  dfrlm = q(i+1,j-1,k ,n)-q(i,j,k,n)
751  dfrmm = q(i+1,j ,k ,n)-q(i,j,k,n)
752  dfrrm = q(i+1,j+1,k ,n)-q(i,j,k,n)
753  !
754  dfllr = q(i-1,j-1,k+1,n)-q(i,j,k,n)
755  dflmr = q(i-1,j ,k+1,n)-q(i,j,k,n)
756  dflrr = q(i-1,j+1,k+1,n)-q(i,j,k,n)
757  dfmlr = q(i ,j-1,k+1,n)-q(i,j,k,n)
758  dfmmr = q(i ,j ,k+1,n)-q(i,j,k,n)
759  dfmrr = q(i ,j+1,k+1,n)-q(i,j,k,n)
760  dfrlr = q(i+1,j-1,k+1,n)-q(i,j,k,n)
761  dfrmr = q(i+1,j ,k+1,n)-q(i,j,k,n)
762  dfrrr = q(i+1,j+1,k+1,n)-q(i,j,k,n)
763  !
764  vmin = min(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
765  & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
766  & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
767  vmax = max(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
768  & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
769  & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
770  !
771  dfx = half*(q(i+1,j,k,n)-q(i-1,j,k,n))
772  dfy = half*(q(i,j+1,k,n)-q(i,j-1,k,n))
773  dfz = half*(q(i,j,k+1,n)-q(i,j,k-1,n))
774  idff = xslope_type / sqrt(dfx**2+dfy**2+dfz**2+tiny(1.0_dp)) ! xslope_type = 1 or 2 dep. on slope_type
775  !-------------------------------------------------------------------
776  ! vmin and vmax are the lower and upper limit of one-cell differences
777  ! so a constant slope with step dq will give vmin=-dq and vmax=+dq,
778  ! while dfx (dff for any direction) will also be = dq, unless we keep
779  ! the 'half' in front of the sqrt. With that kept, the dlim factor
780  ! will remain unity until the slopes differ with a factor of 2. Then,
781  ! if any of the 26 slopes is steeper than twice the gradient slope,
782  ! the gradient will be diminished.
783  !
784  ! idff > 0, so 1/dff = idff < 0.001 * huge is equivalent to dff > 0
785  !-------------------------------------------------------------------
786  slop = merge(min(one,min(abs(vmin),abs(vmax))*idff), one, idff<0.001_dp * huge(1.0_dp))
787  !-------------------------------------------------------------------
788  ! Combine the dlim for different variables into one (experimental)
789  !-------------------------------------------------------------------
790  !dlim = slop
791  !if (n>=2.and.n<=5) dlim = min(dlim,slop)
792  dlim = min(dlim,slop)
793  ! save fluxes for each variable
794  adfx(n) = dfx
795  adfy(n) = dfy
796  adfz(n) = dfz
797  end do
798  do n = 1, self%nv
799  ! use saved fluxes, with common dlim
800  dq(i,j,k,n,1) = dlim*adfx(n)
801  dq(i,j,k,n,2) = dlim*adfy(n)
802  dq(i,j,k,n,3) = dlim*adfz(n)
803  end do
804  end block
805  end do
806  end do
807  end do
808  else
809  write(*,*)'Unknown slope type'
810  stop
811  endif
812  if (detailed_timer) then
813  call trace%end (itimer)
814  else
815  call trace%end
816  end if
817  end associate
818 end subroutine uslope
819 
820 !===============================================================================
821 !> Predictor step, using primitive variables
822 !===============================================================================
823 subroutine trace3d(self)
824  class(hd_t):: self
825  integer ::i, j, k, lb(3), ub(3)
826  integer ::ir=1, iu=2, iv=3, iw=4, ip=5
827  real(dp)::dtdx, dtdy, dtdz
828  real(dp)::r, u, v, w, p, a
829  real(dp)::drx, dux, dvx, dwx, dpx, dax
830  real(dp)::dry, duy, dvy, dwy, dpy, day
831  real(dp)::drz, duz, dvz, dwz, dpz, daz
832  real(dp)::sr0, su0, sv0, sw0, sp0, sa0
833  real(dp)::ff(self%gn(1),3)
834  integer,save:: itimer=0
835  associate(q=>self%q, dq=>self%dq, qp=>self%qp, qm=>self%qm)
836  !-----------------------------------------------------------------------------
837  if (detailed_timer) then
838  call trace%begin ('hd_mod::trace3d', itimer=itimer)
839  else
840  call trace%begin ('hd_mod::trace3d')
841  end if
842  dtdx = self%dtime/self%mesh(1)%d
843  dtdy = self%dtime/self%mesh(2)%d
844  dtdz = self%dtime/self%mesh(3)%d
845  !-----------------------------------------------------------------------------
846  lb = self%mesh%lo
847  ub = self%mesh%uo
848  do k = lb(3),ub(3)
849  do j = lb(2),ub(2)
850  if (allocated(self%force_per_unit_mass)) then
851  do i=lb(1),ub(1)
852  ff(i,1) = self%force_per_unit_mass(i,j,k,1)*self%mesh(1)%d
853  ff(i,2) = self%force_per_unit_mass(i,j,k,2)*self%mesh(2)%d
854  ff(i,3) = self%force_per_unit_mass(i,j,k,3)*self%mesh(3)%d
855  end do
856  else
857  ff(:,:) = 0.0_dp
858  end if
859  do i = lb(1),ub(1)
860  ! Cell centered values
861  r = q(i,j,k,ir)
862  u = q(i,j,k,iu)
863  v = q(i,j,k,iv)
864  w = q(i,j,k,iw)
865  p = q(i,j,k,ip)
866 
867  ! TVD slopes in all 3 directions
868  drx = dq(i,j,k,ir,1)
869  dpx = dq(i,j,k,ip,1)
870  dux = dq(i,j,k,iu,1)
871  dvx = dq(i,j,k,iv,1)
872  dwx = dq(i,j,k,iw,1)
873 
874  dry = dq(i,j,k,ir,2)
875  dpy = dq(i,j,k,ip,2)
876  duy = dq(i,j,k,iu,2)
877  dvy = dq(i,j,k,iv,2)
878  dwy = dq(i,j,k,iw,2)
879 
880  drz = dq(i,j,k,ir,3)
881  dpz = dq(i,j,k,ip,3)
882  duz = dq(i,j,k,iu,3)
883  dvz = dq(i,j,k,iv,3)
884  dwz = dq(i,j,k,iw,3)
885 
886  ! Source terms (including transverse derivatives)
887  sr0 = -u*drx-v*dry-w*drz - (dux+dvy+dwz)*r
888  sp0 = -u*dpx-v*dpy-w*dpz - (dux+dvy+dwz)*gamma*p
889  su0 = -u*dux-v*duy-w*duz - dpx/r + ff(i,1)
890  sv0 = -u*dvx-v*dvy-w*dvz - dpy/r + ff(i,2)
891  sw0 = -u*dwx-v*dwy-w*dwz - dpz/r + ff(i,3)
892 
893  ! Right state at left interface
894  qp(i,j,k,ir,1) = r - half*drx + sr0*dtdx*half
895  qp(i,j,k,ip,1) = p - half*dpx + sp0*dtdx*half
896  qp(i,j,k,iu,1) = u - half*dux + su0*dtdx*half
897  qp(i,j,k,iv,1) = v - half*dvx + sv0*dtdx*half
898  qp(i,j,k,iw,1) = w - half*dwx + sw0*dtdx*half
899  qp(i,j,k,ir,1) = max(smallr, qp(i,j,k,ir,1))
900 
901  ! Left state at right interface
902  qm(i+1,j,k,ir,1) = r + half*drx + sr0*dtdx*half
903  qm(i+1,j,k,ip,1) = p + half*dpx + sp0*dtdx*half
904  qm(i+1,j,k,iu,1) = u + half*dux + su0*dtdx*half
905  qm(i+1,j,k,iv,1) = v + half*dvx + sv0*dtdx*half
906  qm(i+1,j,k,iw,1) = w + half*dwx + sw0*dtdx*half
907  qm(i+1,j,k,ir,1) = max(smallr, qm(i+1,j,k,ir,1))
908 
909  ! Top state at bottom interface
910  qp(i,j,k,ir,2) = r - half*dry + sr0*dtdy*half
911  qp(i,j,k,ip,2) = p - half*dpy + sp0*dtdy*half
912  qp(i,j,k,iu,2) = u - half*duy + su0*dtdy*half
913  qp(i,j,k,iv,2) = v - half*dvy + sv0*dtdy*half
914  qp(i,j,k,iw,2) = w - half*dwy + sw0*dtdy*half
915  qp(i,j,k,ir,2) = max(smallr, qp(i,j,k,ir,2))
916 
917  ! Bottom state at top interface
918  qm(i,j+1,k,ir,2) = r + half*dry + sr0*dtdy*half
919  qm(i,j+1,k,ip,2) = p + half*dpy + sp0*dtdy*half
920  qm(i,j+1,k,iu,2) = u + half*duy + su0*dtdy*half
921  qm(i,j+1,k,iv,2) = v + half*dvy + sv0*dtdy*half
922  qm(i,j+1,k,iw,2) = w + half*dwy + sw0*dtdy*half
923  qm(i,j+1,k,ir,2) = max(smallr, qm(i,j+1,k,ir,2))
924 
925  ! Back state at front interface
926  qp(i,j,k,ir,3) = r - half*drz + sr0*dtdz*half
927  qp(i,j,k,ip,3) = p - half*dpz + sp0*dtdz*half
928  qp(i,j,k,iu,3) = u - half*duz + su0*dtdz*half
929  qp(i,j,k,iv,3) = v - half*dvz + sv0*dtdz*half
930  qp(i,j,k,iw,3) = w - half*dwz + sw0*dtdz*half
931  qp(i,j,k,ir,3) = max(smallr, qp(i,j,k,ir,3))
932 
933  ! Front state at back interface
934  qm(i,j,k+1,ir,3) = r + half*drz + sr0*dtdz*half
935  qm(i,j,k+1,ip,3) = p + half*dpz + sp0*dtdz*half
936  qm(i,j,k+1,iu,3) = u + half*duz + su0*dtdz*half
937  qm(i,j,k+1,iv,3) = v + half*dvz + sv0*dtdz*half
938  qm(i,j,k+1,iw,3) = w + half*dwz + sw0*dtdz*half
939  qm(i,j,k+1,ir,3) = max(smallr, qm(i,j,k+1,ir,3))
940  end do
941  end do
942  end do
943  if (detailed_timer) then
944  call trace%end (itimer)
945  else
946  call trace%end
947  end if
948  end associate
949 end subroutine trace3d
950 
951 !===============================================================================
952 !> Compute fluxes
953 !===============================================================================
954 subroutine cmpflxm(self)
955  class(hd_t):: self
956  ! local variables
957  integer:: i, j, k, lb(3), ub(3)
958  real(dp), dimension(self%mesh(1)%gn,self%nv):: qleft, qright, fgdnv
959  real(dp):: dtds(3)
960  logical:: error
961  integer,save:: itimer=0
962  associate(qp=>self%qp, qm=>self%qm, fl=>self%fl)
963  !-----------------------------------------------------------------------------
964  if (detailed_timer) then
965  call trace%begin ('hd_mod::cmpflxm', itimer=itimer)
966  else
967  call trace%begin ('hd_mod::cmpflxm')
968  end if
969  error = .false.
970  dtds = self%dtime/self%mesh%d
971  lb = self%mesh%li
972  ub = self%mesh%uo
973  qleft=0.0; qright=0.0
974  do k = lb(3),ub(3)
975  do j = lb(2),ub(2)
976  do i = lb(1),ub(1)
977  qleft(i,1) = qm(i,j,k,1,1); qright(i,1) = qp(i,j,k,1,1)
978  qleft(i,2) = qm(i,j,k,2,1); qright(i,2) = qp(i,j,k,2,1)
979  qleft(i,3) = qm(i,j,k,5,1); qright(i,3) = qp(i,j,k,5,1)
980  qleft(i,4) = qm(i,j,k,3,1); qright(i,4) = qp(i,j,k,3,1)
981  qleft(i,5) = qm(i,j,k,4,1); qright(i,5) = qp(i,j,k,4,1)
982  end do
983  if (detailed_timer) call trace_end (itimer)
984  call solver (self,self%mesh(1)%gn,qleft,qright,fgdnv,error,detailed_timer)
985  if (detailed_timer) call trace_begin ('hd_mod::cmpflxm', itimer=itimer)
986  if (error) call error_handler('x', j, k, lb(1), ub(1))
987  do i = lb(1),ub(1)
988  fl(i,j,k,self%idx%d ,1) = fgdnv(i,1)*dtds(1)
989  fl(i,j,k,self%idx%px,1) = fgdnv(i,2)*dtds(1)
990  fl(i,j,k,self%idx%e ,1) = fgdnv(i,3)*dtds(1)
991  fl(i,j,k,self%idx%py,1) = fgdnv(i,4)*dtds(1)
992  fl(i,j,k,self%idx%pz,1) = fgdnv(i,5)*dtds(1)
993  end do
994  end do
995  end do
996  if (any(self%mesh%gn/=self%mesh(1)%gn)) then
997  call cmpflxm_y
998  call cmpflxm_z
999  else
1000  do k = lb(3),ub(3)
1001  do i = lb(1),ub(1)
1002  do j = lb(2),ub(2)
1003  qleft(j,1) = qm(i,j,k,1,2); qright(j,1) = qp(i,j,k,1,2)
1004  qleft(j,2) = qm(i,j,k,3,2); qright(j,2) = qp(i,j,k,3,2)
1005  qleft(j,3) = qm(i,j,k,5,2); qright(j,3) = qp(i,j,k,5,2)
1006  qleft(j,4) = qm(i,j,k,2,2); qright(j,4) = qp(i,j,k,2,2)
1007  qleft(j,5) = qm(i,j,k,4,2); qright(j,5) = qp(i,j,k,4,2)
1008  end do
1009  if (detailed_timer) call trace_end (itimer)
1010  call solver (self,self%mesh(2)%gn,qleft,qright,fgdnv,error,detailed_timer)
1011  if (detailed_timer) call trace_begin ('hd_mod::cmpflxm', itimer=itimer)
1012  if (error) call error_handler('y', i, k, lb(2), ub(2))
1013  do j = lb(2),ub(2)
1014  fl(i,j,k,self%idx%d ,2) = fgdnv(j,1)*dtds(2)
1015  fl(i,j,k,self%idx%py,2) = fgdnv(j,2)*dtds(2)
1016  fl(i,j,k,self%idx%e ,2) = fgdnv(j,3)*dtds(2)
1017  fl(i,j,k,self%idx%px,2) = fgdnv(j,4)*dtds(2)
1018  fl(i,j,k,self%idx%pz,2) = fgdnv(j,5)*dtds(2)
1019  end do
1020  end do
1021  end do
1022  do j = lb(2),ub(2)
1023  do i = lb(1),ub(1)
1024  do k = lb(3),ub(3)
1025  qleft(k,1) = qm(i,j,k,1,3); qright(k,1) = qp(i,j,k,1,3)
1026  qleft(k,2) = qm(i,j,k,4,3); qright(k,2) = qp(i,j,k,4,3)
1027  qleft(k,3) = qm(i,j,k,5,3); qright(k,3) = qp(i,j,k,5,3)
1028  qleft(k,4) = qm(i,j,k,2,3); qright(k,4) = qp(i,j,k,2,3)
1029  qleft(k,5) = qm(i,j,k,3,3); qright(k,5) = qp(i,j,k,3,3)
1030  end do
1031  if (detailed_timer) call trace_end (itimer)
1032  call solver (self,self%mesh(3)%gn,qleft,qright,fgdnv,error,detailed_timer)
1033  if (detailed_timer) call trace_begin ('hd_mod::cmpflxm', itimer=itimer)
1034  if (error) call error_handler('z', i, j, lb(3), ub(3))
1035  do k = lb(3),ub(3)
1036  fl(i,j,k,self%idx%d ,3) = fgdnv(k,1)*dtds(3)
1037  fl(i,j,k,self%idx%pz,3) = fgdnv(k,2)*dtds(3)
1038  fl(i,j,k,self%idx%e ,3) = fgdnv(k,3)*dtds(3)
1039  fl(i,j,k,self%idx%px,3) = fgdnv(k,4)*dtds(3)
1040  fl(i,j,k,self%idx%py,3) = fgdnv(k,5)*dtds(3)
1041  end do
1042  end do
1043  end do
1044  end if
1045  end associate
1046  if (detailed_timer) then
1047  call trace%end (itimer)
1048  else
1049  call trace%end
1050  end if
1051 contains
1052  subroutine cmpflxm_y
1053  ! local variables
1054  real(dp), dimension(self%mesh(2)%gn,self%nv):: qleft, qright, fgdnv
1055  associate(qp=>self%qp, qm=>self%qm, fl=>self%fl)
1056  !-----------------------------------------------------------------------------
1057  !-----------------------------------------------------------------------------
1058  dtds = self%dtime/self%mesh%d ! needed here
1059  lb = self%mesh%li
1060  ub = self%mesh%uo
1061  qleft=0.0; qright=0.0
1062  do k = lb(3),ub(3)
1063  do i = lb(1),ub(1)
1064  qleft(:,1) = qm(i,:,k,1,2); qright(:,1) = qp(i,:,k,1,2)
1065  qleft(:,2) = qm(i,:,k,3,2); qright(:,2) = qp(i,:,k,3,2)
1066  qleft(:,3) = qm(i,:,k,5,2); qright(:,3) = qp(i,:,k,5,2)
1067  qleft(:,4) = qm(i,:,k,2,2); qright(:,4) = qp(i,:,k,2,2)
1068  qleft(:,5) = qm(i,:,k,4,2); qright(:,5) = qp(i,:,k,4,2)
1069  if (detailed_timer) call trace_end (itimer)
1070  call solver (self,self%mesh(2)%gn,qleft,qright,fgdnv,error,detailed_timer)
1071  if (error) call error_handler('y', i, k, lb(2), ub(2))
1072  if (detailed_timer) call trace_begin ('hd_mod::cmpflxm', itimer=itimer)
1073  fl(i,:,k,self%idx%d ,2) = fgdnv(:,1)*dtds(2)
1074  fl(i,:,k,self%idx%py,2) = fgdnv(:,2)*dtds(2)
1075  fl(i,:,k,self%idx%e ,2) = fgdnv(:,3)*dtds(2)
1076  fl(i,:,k,self%idx%px,2) = fgdnv(:,4)*dtds(2)
1077  fl(i,:,k,self%idx%pz,2) = fgdnv(:,5)*dtds(2)
1078  end do
1079  end do
1080  end associate
1081  end subroutine cmpflxm_y
1082  !-----------------------------------------------------------------------------
1083  !-----------------------------------------------------------------------------
1084  subroutine cmpflxm_z
1085  ! local variables
1086  real(dp), dimension(self%mesh(3)%gn,self%nv):: qleft, qright, fgdnv
1087  associate(qp=>self%qp, qm=>self%qm, fl=>self%fl)
1088  !-----------------------------------------------------------------------------
1089  dtds = self%dtime/self%mesh%d ! needed here
1090  lb = self%mesh%li
1091  ub = self%mesh%uo
1092  qleft=0.0; qright=0.0
1093  do j = lb(2),ub(2)
1094  do i = lb(1),ub(1)
1095  qleft(:,1) = qm(i,j,:,1,3); qright(:,1) = qp(i,j,:,1,3)
1096  qleft(:,2) = qm(i,j,:,4,3); qright(:,2) = qp(i,j,:,4,3)
1097  qleft(:,3) = qm(i,j,:,5,3); qright(:,3) = qp(i,j,:,5,3)
1098  qleft(:,4) = qm(i,j,:,2,3); qright(:,4) = qp(i,j,:,2,3)
1099  qleft(:,5) = qm(i,j,:,3,3); qright(:,5) = qp(i,j,:,3,3)
1100  if (detailed_timer) call trace_end (itimer)
1101  call solver (self,self%mesh(3)%gn,qleft,qright,fgdnv,error,detailed_timer)
1102  if (error) call error_handler('z', i, j, lb(3), ub(3))
1103  if (detailed_timer) call trace_begin ('hd_mod::cmpflxm', itimer=itimer)
1104  fl(i,j,:,self%idx%d ,3) = fgdnv(:,1)*dtds(3)
1105  fl(i,j,:,self%idx%pz,3) = fgdnv(:,2)*dtds(3)
1106  fl(i,j,:,self%idx%e ,3) = fgdnv(:,3)*dtds(3)
1107  fl(i,j,:,self%idx%px,3) = fgdnv(:,4)*dtds(3)
1108  fl(i,j,:,self%idx%py,3) = fgdnv(:,5)*dtds(3)
1109  end do
1110  end do
1111  end associate
1112  end subroutine cmpflxm_z
1113  !-----------------------------------------------------------------------------
1114  !-----------------------------------------------------------------------------
1115  subroutine error_handler (label, i1, i2, lb, ub)
1116  class(*), pointer:: link
1117  character(len=*):: label
1118  integer:: i1, i2, lb, ub
1119  class(link_t), pointer:: nbor
1120  !.............................................................................
1121  write (io_unit%log,*) 'solver error: id, i[12] =', self%id, i1, i2
1122  !-----------------------------------------------------------------------------
1123  ! Print the left and right states that give errors
1124  !-----------------------------------------------------------------------------
1125  do i=lb,ub
1126  write (io_unit%log,'(a,1p,2(5e12.3,2x))') &
1127  trim(label)//'solver error:', qleft(i,1:5), qright(i,1:5)
1128  end do
1129  !-----------------------------------------------------------------------------
1130  ! Print a list of neighbors and their time differences
1131  !-----------------------------------------------------------------------------
1132  link => self%link
1133  select type (link)
1134  class is (link_t)
1135  nbor => link%nbor
1136  class default
1137  nullify(nbor)
1138  end select
1139  do while (associated(nbor))
1140  write (io_unit%log,'(a,i9,f10.3,3x,3f6.2)') &
1141  trim(label)//'solver error: nbor, time diff', &
1142  nbor%task%id, (nbor%task%time-self%time)/self%dtime, &
1143  (nbor%task%position-self%position)/self%size
1144  nbor => nbor%next
1145  end do
1146  call mpi%abort ('HLLC')
1147  end subroutine
1148 end subroutine cmpflxm
1149 
1150 !===============================================================================
1151 !> Selection of 1-D solvers
1152 !===============================================================================
1153 subroutine solver (self,ngrid,qleft,qright,fgdnv,error,detailed_timer)
1154  class(hd_t):: self
1155  integer:: ngrid
1156  real(dp),dimension(:,:):: qleft,qright,fgdnv
1157  logical:: error, detailed_timer
1158  !.............................................................................
1159  if(self%riemann.eq.'acoustic')then
1160  call rieman%acoustic(qleft,qright,fgdnv,ngrid,self%idx%e)
1161  else if (self%riemann.eq.'exact')then
1162  call rieman%approx (qleft,qright,fgdnv,ngrid,self%idx%e)
1163  else if (self%riemann.eq.'llf')then
1164  call rieman%llf (qleft,qright,fgdnv,ngrid,self%idx%e)
1165  else if (self%riemann.eq.'hllc')then
1166  call rieman%hllc (qleft,qright,fgdnv,ngrid,self%idx%e,isothermal,error,detailed_timer)
1167  else if (self%riemann.eq.'hllc_a')then
1168  call rieman%hllc_a (qleft,qright,fgdnv,ngrid,self%idx%e)
1169  else if (self%riemann.eq.'hllc_v')then
1170  call rieman%hllc_v (qleft,qright,fgdnv,ngrid,self%idx%e)
1171  else if (self%riemann.eq.'hll')then
1172  call rieman%hll (qleft,qright,fgdnv,ngrid,self%idx%e)
1173  else
1174  write(*,*)'unknown Riemann solver'
1175  stop
1176  end if
1177 end subroutine solver
1178 
1179 !===============================================================================
1180 !===============================================================================
1181 FUNCTION gas_pressure (self) RESULT (pg)
1182  class(hd_t):: self
1183  real, dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
1184  !-----------------------------------------------------------------------------
1185  associate(d => self%mem(:,:,:,self%idx%d ,self%it,1), &
1186  px => self%mem(:,:,:,self%idx%px,self%it,1), &
1187  py => self%mem(:,:,:,self%idx%py,self%it,1), &
1188  pz => self%mem(:,:,:,self%idx%pz,self%it,1), &
1189  e => self%mem(:,:,:,self%idx%e ,self%it,1))
1190  if (isothermal) then
1191  pg = d*self%csound**2
1192  else
1193  pg = (gamma-1d0)*(e-0.5*((px**2+py**2+pz**2)/d))
1194  end if
1195  end associate
1196 END FUNCTION gas_pressure
1197 
1198 !===============================================================================
1199 !===============================================================================
1200 SUBROUTINE gravitational_force (self)
1201  class(hd_t):: self
1202  real, dimension(:,:,:), pointer:: phi, d, et, gradphi1, gradphi2, gradphi3
1203  real, dimension(:,:,:,:), pointer:: p
1204  integer:: i,j,k
1205  real, dimension(3):: wtp1, wtm1
1206  real:: etemp
1207  !.............................................................................
1208  associate(l => self%li, u => self%ui)
1209  wtp1 = 0.50 * 4.0 / 3.0 / self%ds
1210  wtm1 = 0.25 * 1.0 / 3.0 / self%ds
1211  phi => self%mem(:,:,:,self%idx%phi,self%it,1)
1212  d => self%mem(:,:,:,self%idx%d,self%it,1)
1213  et => self%mem(:,:,:,self%idx%e,self%new,1)
1214  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
1215  allocate (gradphi1(self%gn(1),self%gn(2),self%gn(3)), &
1216  gradphi2(self%gn(1),self%gn(2),self%gn(3)), &
1217  gradphi3(self%gn(1),self%gn(2),self%gn(3)))
1218 
1219  do k=l(3),u(3)
1220  do j=l(2),u(2)
1221  do i=l(1),u(1)
1222  ! calculate energy minus the kinetic energy
1223  etemp = et(i,j,k) - 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1224  - 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1225  - 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1226  ! calculate gradient of phi
1227  gradphi1(i,j,k) = wtp1(1) * (phi(i+1,j,k) - phi(i-1,j,k)) &
1228  - wtm1(1) * (phi(i+2,j,k) - phi(i-2,j,k))
1229  gradphi2(i,j,k) = wtp1(2) * (phi(i,j+1,k) - phi(i,j-1,k)) &
1230  - wtm1(2) * (phi(i,j+2,k) - phi(i,j-2,k))
1231  gradphi3(i,j,k) = wtp1(3) * (phi(i,j,k+1) - phi(i,j,k-1)) &
1232  - wtm1(3) * (phi(i,j,k+2) - phi(i,j,k-2))
1233  ! update momentum
1234  p(i,j,k,1) = p(i,j,k,1) - self%dtime * d(i,j,k) * gradphi1(i,j,k)
1235  p(i,j,k,2) = p(i,j,k,2) - self%dtime * d(i,j,k) * gradphi2(i,j,k)
1236  p(i,j,k,3) = p(i,j,k,3) - self%dtime * d(i,j,k) * gradphi3(i,j,k)
1237  ! update total energy with new momentum
1238  et(i,j,k) = etemp + 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1239  + 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1240  + 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1241  end do
1242  end do
1243  end do
1244 
1245  deallocate (gradphi1, gradphi2, gradphi3)
1246  end associate
1247 END SUBROUTINE gravitational_force
1248 
1249 END MODULE
Template version of a module intended for adding extra features, at a level between the basic task an...
Definition: extras_mod.f90:54
RAMSES Godunov solvers.
Definition: hd_mod.f90:17
Generic validation module. The general idea is to be able to compare two runs at critical points in t...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Equation of state module for any sort of tables, provided by a reader.
Definition: eos_mod.f90:4
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
The lock module uses nested locks, to allow versatile use of locks, where a procedure may want to mak...