DISPATCH
mhd_mod.f90
1 !===============================================================================
2 !> RAMSES Godunov solvers, use of guard zones; specifically in HLLD
3 !>
4 !> godunov
5 !> unsplit 1 2 3 4 ui uo ub
6 !> ctoprim:uin +---+---+---+......+---+---+---+ ! in, with guards
7 !> ctoprim:bx +---+---+---+......+---+---+---+ ! down-staggered
8 !> ctoprim:q +---+---+---+......+---+---+ ! up-staggered in b
9 !> uslope:dq +---+---+......+---+ ! both ends loose
10 !> trace3d:qm +---+---+......+---+ ! down-staggered
11 !> trace3d:qp +---+......+---+---+ ! up-shifted
12 !> cmpflxm +---+......+---+ ! need qm & qp
13 !> cmp_mag_flx +---+......+---+ ! is down-staggered
14 !> unew +---+......+ ! up-staggered
15 !===============================================================================
16 MODULE mhd_mod
17  !.............................................................................
18  USE const
19  USE io_mod
20  USE io_unit_mod
21  USE extras_mod
22  USE link_mod
23  USE riemann_mod
24  use hydro_parameters
25  USE trace_mod
26  USE mpi_mod
27  USE omp_mod
28  USE omp_lock_mod
29  USE index_mod
30  USE non_ideal_mod
31  USE omp_timer_mod
32  implicit none
33  PRIVATE
34  type, public, extends(extras_t):: mhd_t
35  logical:: mhd=.true.
36  integer:: idim=3
37  character(len=10):: riemann='hlld'
38  real, pointer, dimension(:,:,:,:):: q ! primitive
39  real, pointer, dimension(:,:,:,:):: bf, em ! face centered B, emf
40  real, pointer, dimension(:,:,:,:,:):: qp, qm ! face centered
41  real, pointer, dimension(:,:,:,:,:):: dq, dq2, dbf ! slopes
42  real, pointer, dimension(:,:,:,:,:):: qlb,qrb,qlt,qrt
43  real, pointer, dimension(:,:,:,:):: uin ! conserved
44  real, pointer, dimension(:,:,:,:):: unew ! new conserved
45  real, pointer, dimension(:,:,:):: ex,ey,ez
46  real, pointer, dimension(:,:,:):: c
47  contains
48  procedure:: init
49  procedure:: update
50  procedure:: gas_pressure
51  end type
52  logical, save:: first_time=.true.
53  logical, save:: detailed_timer=.false.
54  logical, save:: unsigned=.true.
55  integer, save:: verbose=0
56  integer, save:: divb_cleaner=2
57  integer, save:: b_slope_type=0
58  integer, save:: u_slope_type=0
59  character(len=16):: eqn='mhd'
60 CONTAINS
61 
62 !===============================================================================
63 SUBROUTINE init (self)
64  class(mhd_t):: self
65  character(len=10):: riemann='hlld'
66  real, save:: max_dlogd=20.
67  integer, save:: n_dump=0
68  real:: csound=1.
69  integer:: i, iv
70  namelist /ramses_params/ gamma, riemann, slope_type, b_slope_type, &
71  u_slope_type, courant_factor, smallr, smallc, isothermal, csound, &
72  verbose, detailed_timer, unsigned, courant_type, max_dlogd, n_dump, &
73  divb_cleaner
74  character(len=120):: id = &
75  '$Id: 4b4ff75d6255b6cfcb2ea17e1d438c0ccadfb8ba $ solvers/ramses/mhd/mhd_mod.f90'
76 !...............................................................................
77  call trace%begin('mhd_mod%init')
78  call trace%print_id (id)
79  self%nw = 1
80  self%ng = 3
81  self%nv = 8
82  self%kind = 'ramses_mhd_patch'
83  call self%idx%init (5, self%mhd)
84  call self%patch_t%init
85  !-----------------------------------------------------------------------------
86  ! Read solver-relevant parameters and store in the task instance
87  !-----------------------------------------------------------------------------
88  !$omp critical (pde)
89  if (first_time) then
90  first_time = .false.
91  rewind(io_unit%input)
92  read (io_unit%input, ramses_params)
93  isothermal = isothermal .or. gamma==1d0
94  !---------------------------------------------------------------------------
95  ! The HLLD solver throws a divide fault if gamma==1, hence this workaround:
96  !---------------------------------------------------------------------------
97  if (isothermal) gamma=1.0d0
98  if (mpi%master) write (*, ramses_params)
99  end if
100  !$omp end critical (pde)
101  rieman%max_dlogd = max_dlogd
102  self%n_dump = n_dump
103  self%riemann = riemann
104  self%courant = courant_factor
105  self%gamma = gamma
106  self%csound = csound
107  self%staggered = .false.
108  do i=1,3
109  self%mesh(i)%h = 0d0
110  self%mesh(i)%h(self%idx%bx+i-1) = -0.5
111  end do
112  call non_ideal%init
113  self%unsigned(self%idx%d) = unsigned
114  self%unsigned(self%idx%e) = unsigned
115  self%pervolume(self%idx%px) = unsigned
116  self%pervolume(self%idx%py) = unsigned
117  self%pervolume(self%idx%pz) = unsigned
118  call self%gpatch_t%init
119  call self%extras_t%init
120  call trace%end
121 END SUBROUTINE init
122 
123 !===============================================================================
124 !===============================================================================
125 SUBROUTINE update (self)
126  class(mhd_t):: self
127  real, pointer, dimension(:,:,:,:):: p
128  real, allocatable, dimension(:,:,:,:):: em
129  real, allocatable, dimension(:,:,:,:,:):: fl
130  real, allocatable, dimension(:,:,:,:)::fluxambdiff,emfambdiff
131  real, pointer, dimension(:,:,:):: d, e, ux=>null(), uy=>null(), uz=>null()
132  real(dp):: dtdx, dtdy, dtdz
133  integer:: n(3), nv3, i
134  integer, parameter:: ndim=3
135  integer, save:: itimer=0
136  !-----------------------------------------------------------------------------
137  call trace%begin ('mhd_t%update', itimer=itimer)
138  n = self%mesh%ui-self%mesh%li+1
139  if (io%smallr > 0) smallr = io%smallr
140  if (io%courant > 0) self%courant = io%courant
141  !----------------------------------------------------------------------------
142  ! Checks whether to do output
143  !----------------------------------------------------------------------------
144  call self%output
145  !-----------------------------------------------------------------------------
146  ! Memory mapping and allocations. All alloc and dealloc should be done here,
147  ! even if it would be possible to do some of them locally, in the procedures
148  ! being called. Set all values to zero initially, in case some loop extends
149  ! beyond defined values.
150  !-----------------------------------------------------------------------------
151  self%uin => self%mem(:,:,:,:,self%it ,1)
152  self%unew => self%mem(:,:,:,:,self%new,1)
153 
154  n = self%mesh%gn
155  allocate (self%c (n(1),n(2),n(3))) ! self%c = 0_dp
156  allocate (self%Ex (n(1),n(2),n(3))) ! self%Ex = 0_dp
157  allocate (self%Ey (n(1),n(2),n(3))) ! self%Ey = 0_dp
158  allocate (self%Ez (n(1),n(2),n(3))) ! self%Ez = 0_dp
159  allocate (self%q (n(1),n(2),n(3),self%nv)) ! self%q = 0_dp
160  allocate (self%qp (n(1),n(2),n(3),self%nv,ndim)) ! self%qp = 0_dp
161  allocate (self%qm (n(1),n(2),n(3),self%nv,ndim)) ! self%qm = 0_dp
162  allocate (self%dq (n(1),n(2),n(3),self%nv,ndim)) ! self%dq = 0_dp
163  allocate (self%dq2(n(1),n(2),n(3),self%nv,ndim)) ! self%dq2 = 0_dp
164  allocate (self%bf (n(1),n(2),n(3),3)) ! self%bf = 0_dp
165  allocate (self%dbf(n(1),n(2),n(3),3,2)) ! self%dbf = 0_dp
166  allocate (self%qLB(n(1),n(2),n(3),self%nv,ndim)) ! self%qLB = 0_dp
167  allocate (self%qLT(n(1),n(2),n(3),self%nv,ndim)) ! self%qLT = 0_dp
168  allocate (self%qRB(n(1),n(2),n(3),self%nv,ndim)) ! self%qRB = 0_dp
169  allocate (self%qRT(n(1),n(2),n(3),self%nv,ndim)) ! self%qRT = 0_dp
170  d => self%mem(:,:,:, self%idx%d ,self%it,1)
171  e => self%mem(:,:,:, self%idx%e ,self%it,1)
172  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%it,1)
173  !-----------------------------------------------------------------------------
174  ! Compute update using second-order Godunov method
175  !-----------------------------------------------------------------------------
176  call check_small (self) ! Check smallr, smallp
177  call ctoprim(self) ! Primitive variables
178  !self%dq = 0d0 ; self%dq2 = 0d0 ; self%dbf = 0d0
179  call uslope(self,self%bf,self%q,self%dq,self%dbf,slope_type) ! TVD slopes
180  if (verbose>0) print *, 'uslope', self%id, minval(self%dbf), maxval(self%dbf)
181  !-----------------------------------------------------------------------------
182  ! For slope types in the interval 3-4, use the conservative 3.0 for corners
183  ! (cf. RAMSES/CPH)
184  !-----------------------------------------------------------------------------
185  if (slope_type <= 3.0) then
186  self%dq2=self%dq
187  else
188  call uslope(self,self%bf,self%q,self%dq2,self%dbf,3.0) ! TVD slopes
189  endif
190 
191  if (verbose>0) print *, 'uslope', self%id, minval(self%dq), maxval(self%dq)
192  if (verbose>0) print *, 'uslope', self%id, minval(self%dq2), maxval(self%dq2)
193 
194  if (non_ideal%is_used) then
195  allocate (fluxambdiff(self%gn(1),self%gn(2),self%gn(3),3))
196  allocate (emfambdiff(self%gn(1),self%gn(2),self%gn(3),3))
197  call non_ideal%non_ideal_comp(self%mesh,self%idx,self%bf,self%uin,self%q,&
198  fluxambdiff,emfambdiff,self%u_max,self%ds,self%dtime,self%gamma)
199  end if
200  call self%courant_condition ! compute dt
201 
202  call trace3d(self) ! predictor step
203 
204  allocate (fl(self%gn(1),self%gn(2),self%gn(3),self%nv,3))
205  allocate (em(self%gn(1),self%gn(2),self%gn(3),3))
206  fl = 0.0d0
207  em = 0.0d0
208  call cmpflxm(self,self%gn(1),2,3,4,6,7,8,fl(:,:,:,:,1))
209  call cmpflxm(self,self%gn(2),3,2,4,7,6,8,fl(:,:,:,:,2))
210  call cmpflxm(self,self%gn(3),4,2,3,8,6,7,fl(:,:,:,:,3))
211 
212  call cmp_mag_flx(self, self%qRT, self%qRB, self%qLT, self%qLB,2,3,4,6,7,8,em(:,:,:,3),self%gn(3)-1,self%nv)
213  call cmp_mag_flx(self, self%qRT, self%qLT, self%qRB, self%qLB,4,2,3,8,6,7,em(:,:,:,2),self%gn(2)-1,self%nv)
214  call cmp_mag_flx(self, self%qRT, self%qRB, self%qLT, self%qLB,3,4,2,7,8,6,em(:,:,:,1),self%gn(1)-1,self%nv)
215 
216  dtdx = self%dtime/self%ds(1)
217  dtdy = self%dtime/self%ds(2)
218  dtdz = self%dtime/self%ds(3)
219  em(:,:,:,1) = em(:,:,:,1)*dtdy
220  em(:,:,:,2) = em(:,:,:,2)*dtdz
221  em(:,:,:,3) = em(:,:,:,3)*dtdx
222 
223  if (non_ideal%is_used) then
224  call non_ideal%flux_upd(fl,fluxambdiff,self%ds(1),self%dtime,1)
225  call non_ideal%flux_upd(fl,fluxambdiff,self%ds(2),self%dtime,2)
226  call non_ideal%flux_upd(fl,fluxambdiff,self%ds(3),self%dtime,3)
227  call non_ideal%non_ideal_emf_up(em(:,:,:,3),emfambdiff,self%ds(3),self%dtime,3)
228  call non_ideal%non_ideal_emf_up(em(:,:,:,2),emfambdiff,self%ds(2),self%dtime,2)
229  call non_ideal%non_ideal_emf_up(em(:,:,:,1),emfambdiff,self%ds(1),self%dtime,1)
230  end if
231  ! Compute fluxes with appropriate time step
232  fl(:,:,:,:,1) = fl(:,:,:,:,1) * dtdx
233  fl(:,:,:,:,2) = fl(:,:,:,:,2) * dtdy
234  fl(:,:,:,:,3) = fl(:,:,:,:,3) * dtdz
235 
236  !-----------------------------------------------------------------------------
237  ! Update conservative variables
238  !-----------------------------------------------------------------------------
239  associate(i0=>self%mesh(1)%li, j0=>self%mesh(2)%li, k0=>self%mesh(3)%li, &
240  i1=>self%mesh(1)%ui, j1=>self%mesh(2)%ui, k1=>self%mesh(3)%ui, &
241  uin=>self%uin,unew=>self%unew)
242 
243  call self%lock%set ('update')
244  unew = uin
245  nv3 = 5
246  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3 ) = &
247  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3 ) + &
248  fl(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3,1) - &
249  fl(i0+1:i1+1,j0 :j1 ,k0 :k1 ,1:nv3,1) + &
250  fl(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3,2) - &
251  fl(i0 :i1 ,j0+1:j1+1,k0 :k1 ,1:nv3,2) + &
252  fl(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3,3) - &
253  fl(i0 :i1 ,j0 :j1 ,k0+1:k1+1,1:nv3,3)
254 
255  !!!! TO DO !!!!
256  !!! Instead of applying the update in time explicitly in mhd_mod
257  !!! only compute the fluxes and provide them to call timestep%update
258  !!! analogous to what is done in stagger2/mhd_mod.f90
259 
260  ! Compute emf fluxes with appropriate time steps
261 
262 
263  !-----------------------------------------------------------------------------
264  ! Update Bx
265  !-----------------------------------------------------------------------------
266  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,6) = &
267  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,6) + &
268  ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,2) - &
269  em(i0 :i1 ,j0 :j1 ,k0+1:k1+1,2) ) - &
270  ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,3) - &
271  em(i0 :i1 ,j0+1:j1+1,k0 :k1 ,3) )
272 
273  !-----------------------------------------------------------------------------
274  ! Update By
275  !-----------------------------------------------------------------------------
276  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,7) = &
277  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,7) + &
278  ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,3) - &
279  em(i0+1:i1+1,j0 :j1 ,k0 :k1 ,3) ) - &
280  ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,1) - &
281  em(i0 :i1 ,j0 :j1 ,k0+1:k1+1,1) )
282 
283  !-----------------------------------------------------------------------------
284  ! Update Bz
285  !-----------------------------------------------------------------------------
286  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,8) = &
287  unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,8) + &
288  ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,1) - &
289  em(i0 :i1 ,j0+1:j1+1,k0 :k1 ,1) ) - &
290  ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,2) - &
291  em(i0+1:i1+1,j0 :j1 ,k0 :k1 ,2) )
292 
293  if (verbose>0) print *, 'unew', self%id, minval(unew), maxval(unew)
294  end associate
295  !-----------------------------------------------------------------------------
296  ! Reset the energy variable in isothermal cases
297  !-----------------------------------------------------------------------------
298  if (isothermal) &
299  self%mem(:,:,:,self%idx%e,self%new,1) = &
300  self%mem(:,:,:,self%idx%d,self%new,1)
301  !-----------------------------------------------------------------------------
302  ! Add external force as source term
303  !-----------------------------------------------------------------------------
304  if (allocated(self%force_per_unit_mass)) then
305  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
306  do i=1,3
307  p(:,:,:,i) = p(:,:,:,i) + self%dtime*self%force_per_unit_mass(:,:,:,i)*d
308  end do
309  end if
310  !-----------------------------------------------------------------------------
311  ! Deallocate
312  !-----------------------------------------------------------------------------
313  call self%counter_update
314  if (non_ideal%is_used) &
315  deallocate (fluxambdiff,emfambdiff)
316  deallocate (fl,em,self%q,self%qp,self%qm,self%dq,self%dq2,self%bf,self%dbf, &
317  self%qLT,self%qLB,self%qRT,self%qRB,self%c, &
318  self%Ex,self%Ey,self%Ez)
319  !----------------------------------------------------------------------------
320  ! Use one of several methods to clean div(B) from guard zones
321  !----------------------------------------------------------------------------
322  select case (divb_cleaner)
323  case (1)
324  call divb_clean1 (self)
325  case (2)
326  call divb_clean2 (self)
327  end select
328  !.............................................................................
329  call self%lock%unset ('update')
330  call trace%end (itimer)
331  if (rieman%ndiff > 0) then
332  write(io_unit%mpi,*) 'id,ndiff', self%id, rieman%ndiff, rieman%nsolu
333  !$omp atomic write
334  rieman%ndiff = 0
335  !$omp atomic write
336  rieman%nsolu = 0
337  end if
338 END SUBROUTINE update
339 
340 !===============================================================================
341 !> Perform a single iteration of a div(B) clean. The limit of stability is
342 !> eps=1./6. This is intended to prevent build-up of layers of large div(B)
343 !> between patches. An exact solution would be to solve an internal Laplace
344 !> equation, with boundary conditiosn = the boundary value of div(B), and then
345 !> subtract the gradient of the resulting potential.
346 !>
347 !> The problem can also be greatly reduced by using a higher order time
348 !> interpoation scheme, thus making the electric field computed from the
349 !> guard zone values of B and U a better match to the neighbors E-fields.
350 !===============================================================================
351 SUBROUTINE divb_clean1 (self)
352  class(mhd_t):: self
353  !.............................................................................
354  real, dimension(:,:,:), pointer:: bx, by, bz
355  real, dimension(:,:,:), allocatable:: divb
356  real, parameter:: eps=0.11
357  integer:: itimer=0
358  !-----------------------------------------------------------------------------
359  call trace%begin('mhd_t%divb_clean1', itimer=itimer)
360  bx => self%mem(:,:,:,self%idx%bx,self%it,1)
361  by => self%mem(:,:,:,self%idx%by,self%it,1)
362  bz => self%mem(:,:,:,self%idx%bz,self%it,1)
363  allocate (divb(self%mesh(1)%gn,self%mesh(2)%gn,self%mesh(3)%gn))
364  divb = ddup(bx,1) + ddup(by,2) + ddup(bz,3)
365  bx = bx + eps*dddn(divb,1)
366  by = by + eps*dddn(divb,2)
367  bz = bz + eps*dddn(divb,3)
368  deallocate (divb)
369  call trace%end (itimer)
370  return
371 contains
372  !-----------------------------------------------------------------------------
373  function ddup(f,i) result (g)
374  real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
375  integer:: i
376  !.............................................................................
377  g = cshift(f,1,i) - f
378  end function
379  !-----------------------------------------------------------------------------
380  function dddn(f,i) result (g)
381  real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
382  integer:: i
383  !.............................................................................
384  g = f - cshift(f,-1,i)
385  end function
386 END SUBROUTINE divb_clean1
387 
388 !===============================================================================
389 !> Enforce div(B)=0 in the guard zones, by succesively imposing the condition,
390 !> layer-by-layer, on the component perpendicular to each of the patch faces.
391 !===============================================================================
392 SUBROUTINE divb_clean2 (self)
393  class(mhd_t):: self
394  !.............................................................................
395  real, dimension(:,:,:), pointer:: bx, by, bz
396  real, dimension(:,:,:), allocatable:: divb
397  real, parameter:: eps=0.11
398  integer:: l(3), u(3), ix, iy, iz
399  integer:: itimer=0
400  !-----------------------------------------------------------------------------
401  call trace%begin('mhd_t%divb_clean2', itimer=itimer)
402  bx => self%mem(:,:,:,self%idx%bx,self%it,1)
403  by => self%mem(:,:,:,self%idx%by,self%it,1)
404  bz => self%mem(:,:,:,self%idx%bz,self%it,1)
405  l = self%mesh%lb
406  u = self%mesh%ub-1
407  !--- x-direction -------------------------------------------------------------
408  do ix=self%mesh(1)%lo,self%mesh(1)%lb,-1
409  bx(ix ,l(2):u(2),l(3):u(3)) = &
410  bx(ix+1,l(2):u(2),l(3):u(3)) + ( &
411  by(ix,l(2)+1:u(2)+1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
412  bz(ix,l(2) :u(2) ,l(3)+1:u(3)+1) - bz(ix,l(2):u(2),l(3):u(3)))
413  end do
414  do ix=self%mesh(1)%ui,self%mesh(1)%ub-1
415  bx(ix+1,l(2):u(2),l(3):u(3)) = &
416  bx(ix ,l(2):u(2),l(3):u(3)) - ( &
417  by(ix,l(2)+1:u(2)+1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
418  bz(ix,l(2) :u(2) ,l(3)+1:u(3)+1) - bz(ix,l(2):u(2),l(3):u(3)))
419  end do
420  !--- y-direction -------------------------------------------------------------
421  do iy=self%mesh(2)%lo,self%mesh(2)%lb,-1
422  by(l(1):u(1),iy ,l(3):u(3)) = &
423  by(l(1):u(1),iy+1,l(3):u(3)) + ( &
424  bx(l(1)+1:u(1)+1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
425  bz(l(1) :u(1) ,iy,l(3)+1:u(3)+1) - bz(l(1):u(1),iy,l(3):u(3)))
426  end do
427  do iy=self%mesh(2)%ui,self%mesh(2)%ub-1
428  by(l(1):u(1),iy+1,l(3):u(3)) = &
429  by(l(1):u(1),iy ,l(3):u(3)) - ( &
430  bx(l(1)+1:u(1)+1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
431  bz(l(1) :u(1) ,iy,l(3)+1:u(3)+1) - bz(l(1):u(1),iy,l(3):u(3)))
432  end do
433  !--- z-direction -------------------------------------------------------------
434  do iz=self%mesh(3)%lo,self%mesh(3)%lb,-1
435  bz(l(1):u(1),l(2):u(2),iz ) = &
436  bz(l(1):u(1),l(2):u(2),iz+1) + ( &
437  bx(l(1)+1:u(1)+1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
438  by(l(1) :u(1) ,l(2)+1:u(2)+1,iz) - by(l(1):u(1),l(2):u(2),iz))
439  end do
440  do iz=self%mesh(3)%ui,self%mesh(3)%ub-1
441  bz(l(1):u(1),l(2):u(2),iz+1) = &
442  bz(l(1):u(1),l(2):u(2),iz ) - ( &
443  bx(l(1)+1:u(1)+1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
444  by(l(1) :u(1) ,l(2)+1:u(2)+1,iz) - by(l(1):u(1),l(2):u(2),iz))
445  end do
446  call trace%end (itimer)
447 END SUBROUTINE divb_clean2
448 
449 !===============================================================================
450 !===============================================================================
451 subroutine check_small (self)
452  class(mhd_t):: self
453  integer:: i,j,k,lo(3),uo(3),n_smalld,n_smallp
454  real(dp):: c,d,u,v,w,p,ekin,b1,b2,b3,emag
455  integer, save:: nprint=10
456  integer, save:: itimer=0
457  integer:: it, jt, iw
458  logical:: panic
459  associate(uu=>self%uin)
460  !-----------------------------------------------------------------------------
461  call trace%begin('mhd_t%check_small', itimer=itimer, detailed_timer=detailed_timer)
462  n_smalld = 0
463  n_smallp = 0
464  panic = .false.
465  lo = self%mesh%lb
466  uo = self%mesh%ub-1
467  do k=lo(3),uo(3)
468  do j=lo(2),uo(2)
469  do i=lo(1),uo(1)
470  if (uu(i,j,k,1) < smallr) then
471  panic = .true.
472  if (nprint>0) then
473  nprint = nprint-1
474  write (io_unit%log,'(i7,3i5,1p,2e12.3)') self%id,i,j,k,uu(i,j,k,1),smallr
475  end if
476  n_smalld = n_smalld+1
477  uu(i,j,k,1) = smallr
478  end if
479  d = uu(i,j,k,1)
480  c = 1./d
481  u = uu(i,j,k,2)*c
482  v = uu(i,j,k,3)*c
483  w = uu(i,j,k,4)*c
484  ekin = 0.5*d*(u**2+v**2+w**2)
485  b1 = 0.5*(uu(i,j,k,6)+uu(i+1,j,k,6))
486  b2 = 0.5*(uu(i,j,k,7)+uu(i,j+1,k,7))
487  b3 = 0.5*(uu(i,j,k,8)+uu(i,j,k+1,8))
488  emag = 0.5*(b1**2+b2**2+b3**2)
489  if (isothermal) then
490  p = d*self%csound**2
491  uu(i,j,k,5) = p+ekin+emag
492  else
493  p = (uu(i,j,k,5)-ekin)*(gamma-one)
494  if (p < smallr*smallc**2) then
495  p = smallr*smallc**2
496  uu(i,j,k,5) = p/(gamma-one)+ekin+emag
497  n_smallp = n_smallp+1
498  end if
499  end if
500  end do
501  end do
502  end do
503  if (n_smalld+n_smallp > 0) then
504  write (io_unit%log,*) 'small: id,nr,np', &
505  wallclock(),self%id,self%time,n_smalld,n_smallp
506  end if
507  end associate
508  if (panic .and. self%n_dump > 0) then
509  !$omp critical (panic_cr)
510  self%n_dump = self%n_dump-1
511  write (io_unit%dump) self%id, self%gn, self%nv, self%nt, self%nw, self%istep
512  write (io_unit%dump) self%time, self%dtime
513  do iw=1,self%nw
514  do it=1,self%nt
515  jt = self%iit(it)
516  write (io_unit%dump) self%mem(:,:,:,:,jt,iw)
517  end do
518  end do
519  write (io_unit%mpi,*) 'density < SMALLR, dumped', self%id, self%time
520  !$omp end critical (panic_cr)
521  end if
522  call trace%end (itimer, detailed_timer=detailed_timer)
523 end subroutine check_small
524 
525 !===============================================================================
526 !> Conservative to primitive variables. Using scalar scratch variables reduces
527 !> teh need for a large (non-default) OMP_STACKSIZE
528 !===============================================================================
529 subroutine ctoprim(self)
530  class(mhd_t):: self
531  integer::n(4)
532  real(dp):: smallp, c, entho
533  real(dp):: eken, emag, etot, eint, b1, b2, b3
534  integer, save:: itimer=0
535  integer:: i, j, k, loc(3), l(3), u(3)
536  associate(uin=>self%uin,q=>self%q,bf=>self%bf,dt=>self%dtime,c=>self%c)
537  !-----------------------------------------------------------------------------
538  call trace%begin('mhd_t%ctoprim', itimer=itimer, detailed_timer=detailed_timer)
539  if (isothermal) then
540  entho = one
541  else
542  entho = one/max(gamma-one,1e-6)
543  end if
544  smallp = smallr*smallc**2/self%gamma
545  n = shape(uin)-1
546  !
547  !self%bf => uin(:,:,:,6:8) ! Left face centered B
548  uin(:,:,:,1) = max(uin(:,:,:,1),smallr) ! Density
549  q(:,:,:,1) = uin(:,:,:,1) ! Density
550  q(:,:,:,2) = uin(:,:,:,2)/q(:,:,:,1) ! vx
551  q(:,:,:,3) = uin(:,:,:,3)/q(:,:,:,1) ! vy
552  q(:,:,:,4) = uin(:,:,:,4)/q(:,:,:,1) ! vz
553  q(:,:,:,5) = 0.0 ! p_t; initialised for safety
554  q(:,:,:,6) = uin(:,:,:,6) ! Bx
555  q(:,:,:,7) = uin(:,:,:,7) ! By
556  q(:,:,:,8) = uin(:,:,:,8) ! Bz
557  bf(:,:,:,1) = uin(:,:,:,6)
558  bf(:,:,:,2) = uin(:,:,:,7)
559  bf(:,:,:,3) = uin(:,:,:,8)
560  if (verbose>0) print *, 'ctoprim:rh', self%id, minval(q(:,:,:,1)), maxval(q(:,:,:,1))
561  if (verbose>0) print *, 'ctoprim:bf', self%id, minval(bf), maxval(bf)
562 
563  do k=1,self%gn(3)-1
564  do j=1,self%gn(2)-1
565  do i=1,self%gn(1)-1
566  eken = half*(q(i,j,k,2)**2 + q(i,j,k,3)**2 + q(i,j,k,4)**2)*q(i,j,k,1)
567  b1 = half*(q(i,j,k,6)+q(i+1,j,k,6))
568  b2 = half*(q(i,j,k,7)+q(i,j+1,k,7))
569  b3 = half*(q(i,j,k,8)+q(i,j,k+1,8))
570  emag = half*(b1**2 + b2**2 + b3**2)
571  q(i,j,k,5) = (uin(i,j,k,5) - emag - eken)*(gamma-1d0)
572  q(i,j,k,5) = merge(q(i,j,k,1)*self%csound**2, max(q(i,j,k,5),smallp), isothermal)
573  c(i,j,k) = sqrt((q(i,j,k,5)*gamma + 2.0*emag)/q(i,j,k,1))
574  end do
575  end do
576  end do
577  !-----------------------------------------------------------------------------
578  ! Alternative Courant condition formulations
579  !-----------------------------------------------------------------------------
580  if (courant_type==1) then
581  c(:,:,:) = c(:,:,:) + max(abs(q(:,:,:,2)), abs(q(:,:,:,3)), abs(q(:,:,:,4)))
582  else
583  c(:,:,:) = c(:,:,:) + (abs(q(:,:,:,2)) + abs(q(:,:,:,3)) + abs(q(:,:,:,4)))/3.0
584  end if
585  ! Gravity predictor step
586  !q(:,:,:,2:4) = q(:,:,:,idim+1) + Grho*gravin(:,:,:,1:3)*dt*half
587  l = self%mesh%lo
588  u = self%mesh%uo
589  self%u_max = maxval(c(l(1):u(1),l(2):u(2),l(3):u(3)))
590  if (self%get_umax_location) then
591  loc = maxloc(c(l(1):u(1),l(2):u(2),l(3):u(3))) + l-1
592  write (io_unit%log,'(a,i6,f12.6,2(1x,3i4),3f9.4,2(1x,a,1p,4g10.2))') 'u_max:', &
593  self%id, self%time, self%ipos, loc, &
594  (self%mesh(1)%p + self%mesh(1)%r(loc(1))-self%mesh(2)%llc_cart)/self%mesh(1)%b, &
595  (self%mesh(2)%p + self%mesh(2)%r(loc(2))-self%mesh(2)%llc_cart)/self%mesh(2)%b, &
596  (self%mesh(3)%p + self%mesh(3)%r(loc(3))-self%mesh(2)%llc_cart)/self%mesh(3)%b, &
597  'P,B:', q(loc(1),loc(2),loc(3),5:8), &
598  'D,U:', q(loc(1),loc(2),loc(3),1:4)
599  flush (io_unit%log)
600  end if
601  end associate
602  call trace%end (itimer, detailed_timer=detailed_timer)
603 end subroutine ctoprim
604 
605 !===============================================================================
606 !> Slope limiters
607 !===============================================================================
608 subroutine uslope(self,bf,q,dq,dbf,rslope_type)
609  class(mhd_t):: self
610  real, dimension(:,:,:,:):: q
611  real, dimension(:,:,:,:):: bf
612  real, dimension(:,:,:,:,:):: dq
613  real, dimension(:,:,:,:,:):: dbf
614  real:: rslope_type
615  integer::nv, islope_type
616  integer::i, j, k, n, l(3), u(3)
617  real(dp)::dsgn, dlim, dcen, dlft, drgt, slop, xslope_type
618  real(dp)::dfll,dflm,dflr,dfml,dfmm,dfmr,dfrl,dfrm,dfrr
619  real(dp)::dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl
620  real(dp)::dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm
621  real(dp)::dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr
622  real(dp)::vmin,vmax,dfx,dfy,dfz,dff,idff
623  integer,save:: itimer=0
624  !-----------------------------------------------------------------------------
625  call trace%begin ('mhd_t%uslope', itimer=itimer, detailed_timer=detailed_timer)
626  !
627  islope_type = int(rslope_type)
628  l = self%mesh%lb+1
629  u = self%mesh%ub-1
630  nv=self%nv
631  if(slope_type==-1)then
632  do n = 1, nv
633  do k = l(3),u(3)
634  do j = l(2),u(2)
635  do i = l(1),u(1)
636  dq(i,j,k,n,1) = half*(q(i+1,j,k,n) - q(i-1,j,k,n))
637  dq(i,j,k,n,2) = half*(q(i,j+1,k,n) - q(i,j-1,k,n))
638  dq(i,j,k,n,3) = half*(q(i,j,k+1,n) - q(i,j,k-1,n))
639  end do
640  end do
641  end do
642  end do
643  else if(islope_type==0)then
644  dq=0.0_dp
645  return
646  else if (islope_type==1) then ! minmod
647  do n = 1, nv
648  do k = l(3),u(3)
649  do j = l(2),u(2)
650  do i = l(1),u(1)
651  ! slopes in first coordinate direction
652  dlft = q(i ,j,k,n) - q(i-1,j,k,n)
653  drgt = q(i+1,j,k,n) - q(i ,j,k,n)
654  dq(i,j,k,n,1) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
655  ! slopes in second coordinate direction
656  dlft = q(i,j ,k,n) - q(i,j-1,k,n)
657  drgt = q(i,j+1,k,n) - q(i,j ,k,n)
658  dq(i,j,k,n,2) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
659  ! slopes in third coordinate direction
660  dlft = q(i,j,k ,n) - q(i,j,k-1,n)
661  drgt = q(i,j,k+1,n) - q(i,j,k ,n)
662  dq(i,j,k,n,3) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
663  end do
664  end do
665  end do
666  end do
667  else if (islope_type==2) then ! moncen
668  xslope_type= 2 !MIN(slope_type,2)
669  do n = 1, nv
670  do k = l(3),u(3)
671  do j = l(2),u(2)
672  do i = l(1),u(1)
673  ! slopes in first coordinate direction
674  dlft = (q(i ,j,k,n) - q(i-1,j,k,n)) ! = dq
675  drgt = (q(i+1,j,k,n) - q(i ,j,k,n)) ! = dq
676  dcen = half*(dlft+drgt) ! = dq
677  dsgn = sign(one, dcen) ! sign(slope)
678  slop = xslope_type*min(abs(dlft),abs(drgt)) ! slope_type=1: slope=dq, slope_type=2: slope=2*dq
679  slop = merge(slop,zero,(dlft*drgt)>zero) ! if both are non-zero, use the smaller one, or times 2
680  dq(i,j,k,n,1) = dsgn*min(slop,abs(dcen)) ! the smaller of the twice the smaller slope, and the central one ! slopes in second coordinate direction
681  dlft = (q(i,j ,k,n) - q(i,j-1,k,n))
682  drgt = (q(i,j+1,k,n) - q(i,j ,k,n))
683  dcen = half*(dlft+drgt)
684  dsgn = sign(one, dcen)
685  slop = xslope_type*min(abs(dlft),abs(drgt))
686  slop = merge(slop,zero,(dlft*drgt)>zero)
687  dq(i,j,k,n,2) = dsgn*min(slop,abs(dcen))
688  ! slopes in third coordinate direction
689  dlft = (q(i,j,k ,n) - q(i,j,k-1,n))
690  drgt = (q(i,j,k+1,n) - q(i,j,k ,n))
691  dcen = half*(dlft+drgt)
692  dsgn = sign(one, dcen)
693  slop = xslope_type*min(abs(dlft),abs(drgt))
694  slop = merge(slop,zero,(dlft*drgt)>zero)
695  dq(i,j,k,n,3) = dsgn*min(slop,abs(dcen))
696  end do
697  end do
698  end do
699  end do
700 
701  do k = l(3),u(3)
702  do j = l(2),u(2)
703  do i = l(1),u(1)
704  ! Bx
705  dlft = (bf(i,j ,k,1) - bf(i,j-1,k,1))
706  drgt = (bf(i,j+1,k,1) - bf(i,j ,k,1))
707  dcen = half*(dlft+drgt)
708  dsgn = sign(one, dcen)
709  slop = xslope_type*min(abs(dlft),abs(drgt))
710  slop = merge(slop,zero,(dlft*drgt)>zero)
711  dbf(i,j,k,1,1) = dsgn*min(slop,abs(dcen))
712  ! slopes in second coordinate direction
713  dlft = (bf(i,j,k ,1) - bf(i,j,k-1,1))
714  drgt = (bf(i,j,k+1,1) - bf(i,j,k ,1))
715  dcen = half*(dlft+drgt)
716  dsgn = sign(one, dcen)
717  slop = xslope_type*min(abs(dlft),abs(drgt))
718  slop = merge(slop,zero,(dlft*drgt)>zero)
719  dbf(i,j,k,1,2) = dsgn*min(slop,abs(dcen))
720  ! By
721  ! slopes in first coordinate direction
722  dlft = (bf(i ,j,k,2) - bf(i-1,j,k,2))
723  drgt = (bf(i+1,j,k,2) - bf(i ,j,k,2))
724  dcen = half*(dlft+drgt)
725  dsgn = sign(one, dcen)
726  slop = xslope_type*min(abs(dlft),abs(drgt))
727  slop = merge(slop,zero,(dlft*drgt)>zero)
728  dbf(i,j,k,2,1) = dsgn*min(slop,abs(dcen))
729  ! slopes in second coordinate direction
730  dlft = (bf(i,j,k ,2) - bf(i,j,k-1,2))
731  drgt = (bf(i,j,k+1,2) - bf(i,j,k ,2))
732  dcen = half*(dlft+drgt)
733  dsgn = sign(one, dcen)
734  slop = xslope_type*min(abs(dlft),abs(drgt))
735  slop = merge(slop,zero,(dlft*drgt)>zero)
736  dbf(i,j,k,2,2) = dsgn*min(slop,abs(dcen))
737  ! Bz
738  ! slopes in first coordinate direction
739  dlft = (bf(i ,j,k,3) - bf(i-1,j,k,3))
740  drgt = (bf(i+1,j,k,3) - bf(i ,j,k,3))
741  dcen = half*(dlft+drgt)
742  dsgn = sign(one, dcen)
743  slop = xslope_type*min(abs(dlft),abs(drgt))
744  slop = merge(slop,zero,(dlft*drgt)>zero)
745  dbf(i,j,k,3,1) = dsgn*min(slop,abs(dcen))
746  ! slopes in second coordinate direction
747  dlft = (bf(i,j ,k,3) - bf(i,j-1,k,3))
748  drgt = (bf(i,j+1,k,3) - bf(i,j ,k,3))
749  dcen = half*(dlft+drgt)
750  dsgn = sign(one, dcen)
751  slop = xslope_type*min(abs(dlft),abs(drgt))
752  slop = merge(slop,zero,(dlft*drgt)>zero)
753  dbf(i,j,k,3,2) = dsgn*min(slop,abs(dcen))
754  end do
755  end do
756  end do
757 
758  else if (islope_type==3) then ! positivity preserving 3d unsplit slope
759  xslope_type = real(slope_type - 2_dp,kind=dp)
760  do n = 1, nv
761  do k = l(3),u(3)
762  do j = l(2),u(2)
763  do i = l(1),u(1)
764  dflll = q(i-1,j-1,k-1,n)-q(i,j,k,n)
765  dflml = q(i-1,j ,k-1,n)-q(i,j,k,n)
766  dflrl = q(i-1,j+1,k-1,n)-q(i,j,k,n)
767  dfmll = q(i ,j-1,k-1,n)-q(i,j,k,n)
768  dfmml = q(i ,j ,k-1,n)-q(i,j,k,n)
769  dfmrl = q(i ,j+1,k-1,n)-q(i,j,k,n)
770  dfrll = q(i+1,j-1,k-1,n)-q(i,j,k,n)
771  dfrml = q(i+1,j ,k-1,n)-q(i,j,k,n)
772  dfrrl = q(i+1,j+1,k-1,n)-q(i,j,k,n)
773  !
774  dfllm = q(i-1,j-1,k ,n)-q(i,j,k,n)
775  dflmm = q(i-1,j ,k ,n)-q(i,j,k,n)
776  dflrm = q(i-1,j+1,k ,n)-q(i,j,k,n)
777  dfmlm = q(i ,j-1,k ,n)-q(i,j,k,n)
778  dfmmm = q(i ,j ,k ,n)-q(i,j,k,n)
779  dfmrm = q(i ,j+1,k ,n)-q(i,j,k,n)
780  dfrlm = q(i+1,j-1,k ,n)-q(i,j,k,n)
781  dfrmm = q(i+1,j ,k ,n)-q(i,j,k,n)
782  dfrrm = q(i+1,j+1,k ,n)-q(i,j,k,n)
783  !
784  dfllr = q(i-1,j-1,k+1,n)-q(i,j,k,n)
785  dflmr = q(i-1,j ,k+1,n)-q(i,j,k,n)
786  dflrr = q(i-1,j+1,k+1,n)-q(i,j,k,n)
787  dfmlr = q(i ,j-1,k+1,n)-q(i,j,k,n)
788  dfmmr = q(i ,j ,k+1,n)-q(i,j,k,n)
789  dfmrr = q(i ,j+1,k+1,n)-q(i,j,k,n)
790  dfrlr = q(i+1,j-1,k+1,n)-q(i,j,k,n)
791  dfrmr = q(i+1,j ,k+1,n)-q(i,j,k,n)
792  dfrrr = q(i+1,j+1,k+1,n)-q(i,j,k,n)
793  !
794  vmin = min(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
795  & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
796  & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
797  vmax = max(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
798  & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
799  & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
800  !
801  dfx = half*(q(i+1,j,k,n)-q(i-1,j,k,n))
802  dfy = half*(q(i,j+1,k,n)-q(i,j-1,k,n))
803  dfz = half*(q(i,j,k+1,n)-q(i,j,k-1,n))
804  idff = xslope_type / sqrt(dfx**2+dfy**2+dfz**2+tiny(1.0_dp)) ! xslope_type = 1 or 2 dep. on slope_type
805  !
806  ! vmin and vmax are the lower and upper limit of one-cell differences
807  ! so a constant slope with step dq will give vmin=-dq and vmax=+dq,
808  ! while dfx (dff for any direction) will also be = dq, unless we keep
809  ! the 'half' in front of the sqrt. With that kept, the dlim factor
810  ! will remain unity until the slopes differ with a factor of 2. Then,
811  ! if any of the 26 slopes is steeper than twice the gradient slope,
812  ! the gradient will be diminished.
813  !
814  ! idff > 0, so 1/dff = idff < 0.001 * huge is equivalent to dff > 0
815  dlim = merge(min(one,min(abs(vmin),abs(vmax))*idff), one, idff<0.001_dp * huge(1.0_dp))
816  !
817  dq(i,j,k,n,1) = dlim*dfx
818  dq(i,j,k,n,2) = dlim*dfy
819  dq(i,j,k,n,3) = dlim*dfz
820  end do
821  end do
822  end do
823  end do
824  do k = l(3),u(3)
825  do j = l(2),u(2)
826  do i = l(1),u(1)
827  ! Bx
828  dlft = (bf(i,j ,k,1) - bf(i,j-1,k,1))
829  drgt = (bf(i,j+1,k,1) - bf(i,j ,k,1))
830  dcen = half*(dlft+drgt)
831  dsgn = sign(one, dcen)
832  slop = xslope_type*min(abs(dlft),abs(drgt))
833  slop = merge(slop,zero,(dlft*drgt)>zero)
834  dbf(i,j,k,1,1) = dsgn*min(slop,abs(dcen))
835  ! slopes in second coordinate direction
836  dlft = (bf(i,j,k ,1) - bf(i,j,k-1,1))
837  drgt = (bf(i,j,k+1,1) - bf(i,j,k ,1))
838  dcen = half*(dlft+drgt)
839  dsgn = sign(one, dcen)
840  slop = xslope_type*min(abs(dlft),abs(drgt))
841  slop = merge(slop,zero,(dlft*drgt)>zero)
842  dbf(i,j,k,1,2) = dsgn*min(slop,abs(dcen))
843  !By
844  ! slopes in first coordinate direction
845  dlft = (bf(i ,j,k,2) - bf(i-1,j,k,2))
846  drgt = (bf(i+1,j,k,2) - bf(i ,j,k,2))
847  dcen = half*(dlft+drgt)
848  dsgn = sign(one, dcen)
849  slop = xslope_type*min(abs(dlft),abs(drgt))
850  slop = merge(slop,zero,(dlft*drgt)>zero)
851  dbf(i,j,k,2,1) = dsgn*min(slop,abs(dcen))
852  ! slopes in second coordinate direction
853  dlft = (bf(i,j,k ,2) - bf(i,j,k-1,2))
854  drgt = (bf(i,j,k+1,2) - bf(i,j,k ,2))
855  dcen = half*(dlft+drgt)
856  dsgn = sign(one, dcen)
857  slop = xslope_type*min(abs(dlft),abs(drgt))
858  slop = merge(slop,zero,(dlft*drgt)>zero)
859  dbf(i,j,k,2,2) = dsgn*min(slop,abs(dcen))
860  ! Bz
861  ! slopes in first coordinate direction
862  dlft = (bf(i ,j,k,3) - bf(i-1,j,k,3))
863  drgt = (bf(i+1,j,k,3) - bf(i ,j,k,3))
864  dcen = half*(dlft+drgt)
865  dsgn = sign(one, dcen)
866  slop = xslope_type*min(abs(dlft),abs(drgt))
867  slop = merge(slop,zero,(dlft*drgt)>zero)
868  dbf(i,j,k,3,1) = dsgn*min(slop,abs(dcen))
869  ! slopes in second coordinate direction
870  dlft = (bf(i,j ,k,3) - bf(i,j-1,k,3))
871  drgt = (bf(i,j+1,k,3) - bf(i,j ,k,3))
872  dcen = half*(dlft+drgt)
873  dsgn = sign(one, dcen)
874  slop = xslope_type*min(abs(dlft),abs(drgt))
875  slop = merge(slop,zero,(dlft*drgt)>zero)
876  dbf(i,j,k,3,2) = dsgn*min(slop,abs(dcen))
877  end do
878  end do
879  end do
880 
881  else
882  write(*,*)'Unknown slope type'
883  stop
884  endif
885  call trace%end (itimer, detailed_timer=detailed_timer)
886 end subroutine uslope
887 
888 !===============================================================================
889 !> Predictor step
890 !===============================================================================
891 SUBROUTINE trace3d(self)
892  IMPLICIT NONE
893  class(mhd_t):: self
894  real(dp):: dx, dy, dz, dt
895  integer :: i, j, k, n
896  integer, parameter ::ir=1, iu=2, iv=3, iw=4, ip=5, ia=6, ib=7, ic=8
897  real(dp)::dtdx, dtdy, dtdz, smallp
898  real(dp) :: invr
899  real(dp) :: r, u, v, w, p, a, b, c
900  real(dp) :: ell, elr, erl, err
901  real(dp) :: fll, flr, frl, frr
902  real(dp) :: gll, glr, grl, grr
903  real(dp) :: drx, dux, dvx, dwx, dpx, dax, dbx, dcx
904  real(dp) :: dry, duy, dvy, dwy, dpy, day, dby, dcy
905  real(dp) :: drz, duz, dvz, dwz, dpz, daz, dbz, dcz
906  real(dp) :: drx2,dpx2,dry2,dpy2,drz2,dpz2
907  real(dp) :: sr0, su0, sv0, sw0, sp0, sa0, sb0, sc0
908  real(dp) :: al, ar, bl, br, cl, cr
909  real(dp) :: daly, dary, dalz, darz
910  real(dp) :: dblx, dbrx, dblz, dbrz
911  real(dp) :: dclx, dcrx, dcly, dcry
912  real(dp) :: sal0, sar0, sbl0, sbr0, scl0, scr0
913  real(dp), dimension(1:self%gn(1),1:self%gn(2),1:self%gn(3)):: rho_t, p_t
914  real(dp) :: eff_gamma, off
915  integer, save:: itimer(4)=0
916  real:: ff(self%gn(1),3)
917  !-----------------------------------------------------------------------------
918  associate(q=>self%q,bf=>self%bf,dq=>self%dq,dq2=>self%dq2,qm=>self%qm,qp=>self%qp, &
919  dbf=>self%dbf,qrt=>self%qRT,qrb=>self%qRB,qlt=>self%qLT,qlb=>self%qLB, &
920  bx=>self%bf(:,:,:,1),by=>self%bf(:,:,:,2),bz=>self%bf(:,:,:,3), &
921  ux=>self%q(:,:,:,2),uy=>self%q(:,:,:,3),uz=>self%q(:,:,:,4), &
922  ex=>self%Ex,ey=>self%Ey,ez=>self%Ez)
923  !
924  call trace%begin ('mhd_t%trace3d', itimer=itimer(1), detailed_timer=detailed_timer)
925  su0=0.0_dp; sv0=0.0_dp; sw0=0.0_dp
926  !
927  dtdx = self%dtime/self%ds(1)
928  dtdy = self%dtime/self%ds(2)
929  dtdz = self%dtime/self%ds(3)
930  smallp = smallr*smallc**2/gamma
931  !
932  ! E = - UxB = BxU
933  ex = ydn(zdn(uy))*ydn(bz) - zdn(ydn(uz))*zdn(by)
934  ey = zdn(xdn(uz))*zdn(bx) - xdn(zdn(ux))*xdn(bz)
935  ez = xdn(ydn(ux))*xdn(by) - ydn(xdn(uy))*ydn(bx)
936 
937  ! fill positive definite variables with reasonable values in "ghost-zones"
938  rho_t(:,:,:) = q(:,:,:,ir)
939  p_t(:,:,:) = q(:,:,:,ip)
940 
941  ff = 0.0
942 
943  !
944  !call trace%end (itimer(1))
945  !call trace%begin ('mhd_t%trace3d(2)', itimer=itimer(2))
946  eff_gamma = gamma
947  do k=2,self%gn(3)-1
948  do j=2,self%gn(2)-1
949  if (allocated(self%force_per_unit_mass)) then
950  ff(:,1) = self%force_per_unit_mass(:,j,k,1)
951  ff(:,2) = self%force_per_unit_mass(:,j,k,2)
952  ff(:,3) = self%force_per_unit_mass(:,j,k,3)
953  end if
954  if (allocated(self%force_per_unit_volume)) then
955  ff(:,1) = self%force_per_unit_mass(:,j,k,1)/q(:,j,k,ir)
956  ff(:,2) = self%force_per_unit_mass(:,j,k,2)/q(:,j,k,ir)
957  ff(:,3) = self%force_per_unit_mass(:,j,k,3)/q(:,j,k,ir)
958  end if
959  !!dir! simd
960  do i=2,self%gn(1)-1
961 ! associate ( &
962  r=q(i,j,k,ir)
963  u=q(i,j,k,iu)
964  v=q(i,j,k,iv)
965  w=q(i,j,k,iw)
966  p=q(i,j,k,ip)
967  a=q(i,j,k,ia)
968  b=q(i,j,k,ib)
969  c=q(i,j,k,ic)
970  ! Face centered variables
971  al = bf(i,j,k,1)
972  bl = bf(i,j,k,2)
973  cl = bf(i,j,k,3)
974  ar = bf(i+1,j ,k ,1)
975  br = bf(i ,j+1,k ,2)
976  cr = bf(i ,j ,k+1,3)
977  ! Edge centered electric field in X, Y and Z directions
978  ell = ex(i,j ,k )
979  elr = ex(i,j ,k+1)
980  erl = ex(i,j+1,k )
981  err = ex(i,j+1,k+1)
982  fll = ey(i ,j,k )
983  flr = ey(i ,j,k+1)
984  frl = ey(i+1,j,k )
985  frr = ey(i+1,j,k+1)
986  gll = ez(i ,j ,k)
987  glr = ez(i ,j+1,k)
988  grl = ez(i+1,j ,k)
989  grr = ez(i+1,j+1,k)
990  !
991  ! Cell centered TVD slopes in X, Y and Z directions
992  drx = half * dq(i,j,k,ir,1)
993  drx2= half *dq2(i,j,k,ir,1)
994  dux = half * dq(i,j,k,iu,1)
995  dvx = half * dq(i,j,k,iv,1)
996  dwx = half * dq(i,j,k,iw,1)
997  dpx = half * dq(i,j,k,ip,1)
998  dpx2= half *dq2(i,j,k,ip,1)
999  dbx = half * dq(i,j,k,ib,1)
1000  dcx = half * dq(i,j,k,ic,1)
1001  !
1002  dry = half * dq(i,j,k,ir,2)
1003  dry2= half *dq2(i,j,k,ir,2)
1004  duy = half * dq(i,j,k,iu,2)
1005  dvy = half * dq(i,j,k,iv,2)
1006  dwy = half * dq(i,j,k,iw,2)
1007  dpy = half * dq(i,j,k,ip,2)
1008  dpy2= half *dq2(i,j,k,ip,2)
1009  duy = half * dq(i,j,k,iu,2)
1010  day = half * dq(i,j,k,ia,2)
1011  dcy = half * dq(i,j,k,ic,2)
1012  !
1013  drz = half * dq(i,j,k,ir,3)
1014  drz2= half *dq2(i,j,k,ir,3)
1015  duz = half * dq(i,j,k,iu,3)
1016  dvz = half * dq(i,j,k,iv,3)
1017  dwz = half * dq(i,j,k,iw,3)
1018  dpz = half * dq(i,j,k,ip,3)
1019  dpz2= half *dq2(i,j,k,ip,3)
1020  daz = half * dq(i,j,k,ia,3)
1021  dbz = half * dq(i,j,k,ib,3)
1022  !
1023  ! Face centered TVD slopes in transverse directions
1024  daly = half * dbf(i ,j ,k ,1,1)
1025  dary = half * dbf(i+1,j ,k ,1,1)
1026  dalz = half * dbf(i ,j ,k ,1,2)
1027  darz = half * dbf(i+1,j ,k ,1,2)
1028 
1029  dblx = half * dbf(i ,j ,k ,2,1)
1030  dbrx = half * dbf(i ,j+1,k ,2,1)
1031  dblz = half * dbf(i ,j ,k ,2,2)
1032  dbrz = half * dbf(i ,j+1,k ,2,2)
1033 
1034  dclx = half * dbf(i ,j ,k ,3,1)
1035  dcrx = half * dbf(i ,j ,k+1,3,1)
1036  dcly = half * dbf(i ,j ,k ,3,2)
1037  dcry = half * dbf(i ,j ,k+1,3,2)
1038 
1039  !
1040  ! Face-centered predicted states
1041  sal0 = +(glr-gll)*dtdy*half -(flr-fll)*dtdz*half
1042  sar0 = +(grr-grl)*dtdy*half -(frr-frl)*dtdz*half
1043  sbl0 = -(grl-gll)*dtdx*half +(elr-ell)*dtdz*half
1044  sbr0 = -(grr-glr)*dtdx*half +(err-erl)*dtdz*half
1045  scl0 = +(frl-fll)*dtdx*half -(erl-ell)*dtdy*half
1046  scr0 = +(frr-flr)*dtdx*half -(err-elr)*dtdy*half
1047  !
1048  al = al + sal0
1049  ar = ar + sar0
1050  bl = bl + sbl0
1051  br = br + sbr0
1052  cl = cl + scl0
1053  cr = cr + scr0
1054  !
1055  ! Source terms (including transverse derivatives)
1056  invr = one / r
1057  sr0 = (-u*drx-dux*r)*dtdx + (-v*dry-dvy*r)*dtdy + (-w*drz-dwz*r)*dtdz
1058  su0 = (-u*dux-(dpx+b*dbx+c*dcx)*invr)*dtdx + (-v*duy+b*day*invr)*dtdy + (-w*duz+c*daz*invr)*dtdz
1059  sv0 = (-u*dvx+a*dbx*invr)*dtdx + (-v*dvy-(dpy+a*day+c*dcy)*invr)*dtdy + (-w*dvz+c*dbz*invr)*dtdz
1060  sw0 = (-u*dwx+a*dcx*invr)*dtdx + (-v*dwy+b*dcy*invr)*dtdy + (-w*dwz-(dpz+a*daz+b*dbz)*invr)*dtdz
1061  sp0 = (-u*dpx-dux*eff_gamma*p)*dtdx + (-v*dpy-dvy*eff_gamma*p)*dtdy + (-w*dpz-dwz*eff_gamma*p)*dtdz
1062  !
1063  ! Add external force
1064  su0 = su0 + ff(i,1)*half*self%dtime
1065  sv0 = sv0 + ff(i,2)*half*self%dtime
1066  sw0 = sw0 + ff(i,3)*half*self%dtime
1067  !
1068  ! Evolve entropy instead of pressure
1069  !sS0 = (-u*dSx)*dtdx + (-v*dSy)*dtdy + (-w*dSz)*dtdz
1070  !
1071 ! ! Save predicted cell-centered pseudo-entropy; we are borrowing space in q(:,:,:,:,ip),
1072 ! ! which is not used below this point. In the calling routine this is again copied over
1073 ! ! into uin, and is thus made available to godfine1, through that array.
1074 ! if (entropy_fix) then
1075 ! if ((i==1.or.i==2).and.(j==1.or.j==2).and.(k==1.or.k==2)) then
1076 ! q(:,i,j,k,ip) = S + 2.0*sS0
1077 ! if (verbose>0) then
1078 ! do ig=1,ngrid
1079 ! if (in_debug_region(ind_grid(ig))) then
1080 ! print'(3hMKE,i8,i3,3f10.6,1p,12e11.3)', &
1081 ! ind_grid(ig),(i-1)+(j-1)*2+(k-1)*4, &
1082 ! xg(ind_grid(ig),1)+(i-1.5)/32.0, &
1083 ! xg(ind_grid(ig),2)+(j-1.5)/32.0, &
1084 ! xg(ind_grid(ig),3)+(k-1.5)/32.0, &
1085 ! log(p(ig)-dpx(ig)) - gamma*log(r(ig)-drx(ig)), &
1086 ! log(p(ig))-gamma*log(r(ig)), &
1087 ! log(p(ig)+dpx(ig)) - gamma*log(r(ig)+drx(ig)), &
1088 ! q(ig,i,j,k,ip), &
1089 ! dSx(ig),dSy(ig),dSz(ig), &
1090 ! u(ig),v(ig),w(ig),dtdx
1091 ! end if
1092 ! end do
1093 ! end if
1094 ! end if
1095 ! end if
1096  !
1097 
1098  ! Cell-centered predicted states
1099  r = r + sr0
1100  r = max(r,smallr)
1101  rho_t(i,j,k) = r
1102  u = u + su0
1103  v = v + sv0
1104  w = w + sw0
1105  p = p + sp0
1106  p = max(p,smallp)
1107  p_t(i,j,k) = p
1108  a = half*(al+ar)
1109  b = half*(bl+br)
1110  c = half*(cl+cr)
1111  !
1112  ! Face averaged right state at left interface
1113  qp(i,j,k,ir,1) = r - drx
1114  qp(i,j,k,iu,1) = u - dux
1115  qp(i,j,k,iv,1) = v - dvx
1116  qp(i,j,k,iw,1) = w - dwx
1117  qp(i,j,k,ip,1) = p - dpx
1118  qp(i,j,k,ia,1) = al
1119  qp(i,j,k,ib,1) = b - dbx
1120  qp(i,j,k,ic,1) = c - dcx
1121  !
1122  ! Face averaged left state at right interface
1123  qm(i,j,k,ir,1) = r + drx
1124  qm(i,j,k,iu,1) = u + dux
1125  qm(i,j,k,iv,1) = v + dvx
1126  qm(i,j,k,iw,1) = w + dwx
1127  qm(i,j,k,ip,1) = p + dpx
1128  qm(i,j,k,ia,1) = ar
1129  qm(i,j,k,ib,1) = b + dbx
1130  qm(i,j,k,ic,1) = c + dcx
1131  !
1132  ! Face averaged top state at bottom interface
1133  qp(i,j,k,ir,2) = r - dry
1134  qp(i,j,k,iu,2) = u - duy
1135  qp(i,j,k,iv,2) = v - dvy
1136  qp(i,j,k,iw,2) = w - dwy
1137  qp(i,j,k,ip,2) = p - dpy
1138  qp(i,j,k,ia,2) = a - day
1139  qp(i,j,k,ib,2) = bl
1140  qp(i,j,k,ic,2) = c - dcy
1141  !
1142  ! Face averaged bottom state at top interface
1143  qm(i,j,k,ir,2) = r + dry
1144  qm(i,j,k,iu,2) = u + duy
1145  qm(i,j,k,iv,2) = v + dvy
1146  qm(i,j,k,iw,2) = w + dwy
1147  qm(i,j,k,ip,2) = p + dpy
1148  qm(i,j,k,ia,2) = a + day
1149  qm(i,j,k,ib,2) = br
1150  qm(i,j,k,ic,2) = c + dcy
1151  !
1152  ! Face averaged front state at back interface
1153  qp(i,j,k,ir,3) = r - drz
1154  qp(i,j,k,iu,3) = u - duz
1155  qp(i,j,k,iv,3) = v - dvz
1156  qp(i,j,k,iw,3) = w - dwz
1157  qp(i,j,k,ip,3) = p - dpz
1158  qp(i,j,k,ia,3) = a - daz
1159  qp(i,j,k,ib,3) = b - dbz
1160  qp(i,j,k,ic,3) = cl
1161  !
1162  ! Face averaged back state at front interface
1163  qm(i,j,k,ir,3) = r + drz
1164  qm(i,j,k,iu,3) = u + duz
1165  qm(i,j,k,iv,3) = v + dvz
1166  qm(i,j,k,iw,3) = w + dwz
1167  qm(i,j,k,ip,3) = p + dpz
1168  qm(i,j,k,ia,3) = a + daz
1169  qm(i,j,k,ib,3) = b + dbz
1170  qm(i,j,k,ic,3) = cr
1171  !
1172  ! X-edge averaged right-top corner state (RT->LL)
1173  qrt(i,j,k,ir,1) = r + (+dry2+drz2)
1174  qrt(i,j,k,iu,1) = u + (+duy+duz)
1175  qrt(i,j,k,iv,1) = v + (+dvy+dvz)
1176  qrt(i,j,k,iw,1) = w + (+dwy+dwz)
1177  qrt(i,j,k,ip,1) = p + (+dpy2+dpz2)
1178  qrt(i,j,k,ia,1) = a + (+day+daz)
1179  qrt(i,j,k,ib,1) = br+ ( +dbrz)
1180  qrt(i,j,k,ic,1) = cr+ (+dcry )
1181  !
1182  ! X-edge averaged right-bottom corner state (RB->LR)
1183  qrb(i,j,k,ir,1) = r + (+dry2-drz2)
1184  qrb(i,j,k,iu,1) = u + (+duy-duz)
1185  qrb(i,j,k,iv,1) = v + (+dvy-dvz)
1186  qrb(i,j,k,iw,1) = w + (+dwy-dwz)
1187  qrb(i,j,k,ip,1) = p + (+dpy2-dpz2)
1188  qrb(i,j,k,ia,1) = a + (+day-daz)
1189  qrb(i,j,k,ib,1) = br+ ( -dbrz)
1190  qrb(i,j,k,ic,1) = cl+ (+dcly )
1191  !
1192  ! X-edge averaged left-top corner state (LT->RL)
1193  qlt(i,j,k,ir,1) = r + (-dry2+drz2)
1194  qlt(i,j,k,iu,1) = u + (-duy+duz)
1195  qlt(i,j,k,iv,1) = v + (-dvy+dvz)
1196  qlt(i,j,k,iw,1) = w + (-dwy+dwz)
1197  qlt(i,j,k,ip,1) = p + (-dpy2+dpz2)
1198  qlt(i,j,k,ia,1) = a + (-day+daz)
1199  qlt(i,j,k,ib,1) = bl+ ( +dblz)
1200  qlt(i,j,k,ic,1) = cr+ (-dcry )
1201  !
1202  ! X-edge averaged left-bottom corner state (LB->RR)
1203  qlb(i,j,k,ir,1) = r + (-dry2-drz2)
1204  qlb(i,j,k,iu,1) = u + (-duy-duz)
1205  qlb(i,j,k,iv,1) = v + (-dvy-dvz)
1206  qlb(i,j,k,iw,1) = w + (-dwy-dwz)
1207  qlb(i,j,k,ip,1) = p + (-dpy2-dpz2)
1208  qlb(i,j,k,ia,1) = a + (-day-daz)
1209  qlb(i,j,k,ib,1) = bl+ ( -dblz)
1210  qlb(i,j,k,ic,1) = cl+ (-dcly )
1211  !
1212  ! Y-edge averaged right-top corner state (RT->LL)
1213  qrt(i,j,k,ir,2) = r + (+drx2+drz2)
1214  qrt(i,j,k,iu,2) = u + (+dux+duz)
1215  qrt(i,j,k,iv,2) = v + (+dvx+dvz)
1216  qrt(i,j,k,iw,2) = w + (+dwx+dwz)
1217  qrt(i,j,k,ip,2) = p + (+dpx2+dpz2)
1218  qrt(i,j,k,ia,2) = ar+ ( +darz)
1219  qrt(i,j,k,ib,2) = b + (+dbx+dbz)
1220  qrt(i,j,k,ic,2) = cr+ (+dcrx )
1221  !
1222  ! Y-edge averaged right-bottom corner state (RB->LR)
1223  qrb(i,j,k,ir,2) = r + (+drx2-drz2)
1224  qrb(i,j,k,iu,2) = u + (+dux-duz)
1225  qrb(i,j,k,iv,2) = v + (+dvx-dvz)
1226  qrb(i,j,k,iw,2) = w + (+dwx-dwz)
1227  qrb(i,j,k,ip,2) = p + (+dpx2-dpz2)
1228  qrb(i,j,k,ia,2) = ar+ ( -darz)
1229  qrb(i,j,k,ib,2) = b + (+dbx-dbz)
1230  qrb(i,j,k,ic,2) = cl+ (+dclx )
1231  !
1232  ! Y-edge averaged left-top corner state (LT->RL)
1233  qlt(i,j,k,ir,2) = r + (-drx2+drz2)
1234  qlt(i,j,k,iu,2) = u + (-dux+duz)
1235  qlt(i,j,k,iv,2) = v + (-dvx+dvz)
1236  qlt(i,j,k,iw,2) = w + (-dwx+dwz)
1237  qlt(i,j,k,ip,2) = p + (-dpx2+dpz2)
1238  qlt(i,j,k,ia,2) = al+ ( +dalz)
1239  qlt(i,j,k,ib,2) = b + (-dbx+dbz)
1240  qlt(i,j,k,ic,2) = cr+ (-dcrx )
1241  !
1242  ! Y-edge averaged left-bottom corner state (LB->RR)
1243  qlb(i,j,k,ir,2) = r + (-drx2-drz2)
1244  qlb(i,j,k,iu,2) = u + (-dux-duz)
1245  qlb(i,j,k,iv,2) = v + (-dvx-dvz)
1246  qlb(i,j,k,iw,2) = w + (-dwx-dwz)
1247  qlb(i,j,k,ip,2) = p + (-dpx2-dpz2)
1248  qlb(i,j,k,ia,2) = al+ ( -dalz)
1249  qlb(i,j,k,ib,2) = b + (-dbx-dbz)
1250  qlb(i,j,k,ic,2) = cl+ (-dclx )
1251  !
1252  ! Z-edge averaged right-top corner state (RT->LL)
1253  qrt(i,j,k,ir,3) = r + (+drx2+dry2)
1254  qrt(i,j,k,iu,3) = u + (+dux+duy)
1255  qrt(i,j,k,iv,3) = v + (+dvx+dvy)
1256  qrt(i,j,k,iw,3) = w + (+dwx+dwy)
1257  qrt(i,j,k,ip,3) = p + (+dpx2+dpy2)
1258  qrt(i,j,k,ia,3) = ar+ ( +dary)
1259  qrt(i,j,k,ib,3) = br+ (+dbrx )
1260  qrt(i,j,k,ic,3) = c + (+dcx+dcy)
1261  !
1262  ! Z-edge averaged right-bottom corner state (RB->LR)
1263  qrb(i,j,k,ir,3) = r + (+drx2-dry2)
1264  qrb(i,j,k,iu,3) = u + (+dux-duy)
1265  qrb(i,j,k,iv,3) = v + (+dvx-dvy)
1266  qrb(i,j,k,iw,3) = w + (+dwx-dwy)
1267  qrb(i,j,k,ip,3) = p + (+dpx2-dpy2)
1268  qrb(i,j,k,ia,3) = ar+ ( -dary)
1269  qrb(i,j,k,ib,3) = bl+ (+dblx )
1270  qrb(i,j,k,ic,3) = c + (+dcx-dcy)
1271  !
1272  ! Z-edge averaged left-top corner state (LT->RL)
1273  qlt(i,j,k,ir,3) = r + (-drx2+dry2)
1274  qlt(i,j,k,iu,3) = u + (-dux+duy)
1275  qlt(i,j,k,iv,3) = v + (-dvx+dvy)
1276  qlt(i,j,k,iw,3) = w + (-dwx+dwy)
1277  qlt(i,j,k,ip,3) = p + (-dpx2+dpy2)
1278  qlt(i,j,k,ia,3) = al+ ( +daly)
1279  qlt(i,j,k,ib,3) = br+ (-dbrx )
1280  qlt(i,j,k,ic,3) = c + (-dcx+dcy)
1281  !
1282  ! Z-edge averaged left-bottom corner state (LB->RR)
1283  qlb(i,j,k,ir,3) = r + (-drx2-dry2)
1284  qlb(i,j,k,iu,3) = u + (-dux-duy)
1285  qlb(i,j,k,iv,3) = v + (-dvx-dvy)
1286  qlb(i,j,k,iw,3) = w + (-dwx-dwy)
1287  qlb(i,j,k,ip,3) = p + (-dpx2-dpy2)
1288  qlb(i,j,k,ia,3) = al+ ( -daly)
1289  qlb(i,j,k,ib,3) = bl+ (-dblx )
1290  qlb(i,j,k,ic,3) = c + (-dcx-dcy)
1291  end do
1292  end do
1293  end do
1294  !
1295  !call trace%end (itimer(2))
1296  !call trace%begin ('mhd_t%trace3d(3)', itimer=itimer(3))
1297  do k=2,self%gn(3)-1
1298  do j=2,self%gn(2)-1
1299  do i=2,self%gn(1)-1
1300  r = rho_t(i,j,k)
1301  p = p_t(i,j,k)
1302 
1303  ! X-edge averaged right-top corner state (RT->LL)
1304  !qRT(i,j,k,ir,1) = r + (+dry+drz)
1305  drx = sqrt(sqrt(r*rho_t(i,j+1,k)*rho_t(i,j,k+1)*rho_t(i,j+1,k+1))) - r
1306  dpx = sqrt(sqrt(p* p_t(i,j+1,k)* p_t(i,j,k+1)* p_t(i,j+1,k+1))) - p
1307  dry = qrt(i,j,k,ir,1) - r
1308  dpy = qrt(i,j,k,ip,1) - p
1309  qrt(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1310  qrt(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1311  ! X-edge averaged right-bottom corner state (RB->LR)
1312  !qRB(i,j,k,ir,1) = r + (+dry-drz)
1313  drx = sqrt(sqrt(r*rho_t(i,j+1,k)*rho_t(i,j,k-1)*rho_t(i,j+1,k-1))) - r
1314  dpx = sqrt(sqrt(p* p_t(i,j+1,k)* p_t(i,j,k-1)* p_t(i,j+1,k-1))) - p
1315  dry = qrb(i,j,k,ir,1) - r
1316  dpy = qrb(i,j,k,ip,1) - p
1317  qrb(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1318  qrb(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1319  ! X-edge averaged left-top corner state (LT->RL)
1320  !qLT(i,j,k,ir,1) = r + (-dry+drz)
1321  drx = sqrt(sqrt(r*rho_t(i,j-1,k)*rho_t(i,j,k+1)*rho_t(i,j-1,k+1))) - r
1322  dpx = sqrt(sqrt(p* p_t(i,j-1,k)* p_t(i,j,k+1)* p_t(i,j-1,k+1))) - p
1323  dry = qlt(i,j,k,ir,1) - r
1324  dpy = qlt(i,j,k,ip,1) - p
1325  qlt(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1326  qlt(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1327  ! X-edge averaged left-bottom corner state (LB->RR)
1328  !qLB(i,j,k,ir,1) = r + (-dry-drz)
1329  drx = sqrt(sqrt(r*rho_t(i,j-1,k)*rho_t(i,j,k-1)*rho_t(i,j-1,k-1))) - r
1330  dpx = sqrt(sqrt(p* p_t(i,j-1,k)* p_t(i,j,k-1)* p_t(i,j-1,k-1))) - p
1331  dry = qlb(i,j,k,ir,1) - r
1332  dpy = qlb(i,j,k,ip,1) - p
1333  qlb(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1334  qlb(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1335  ! Y-edge averaged right-top corner state (RT->LL)
1336  !qRT(i,j,k,ir,2) = r + (+drx+drz)
1337  drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j,k+1)*rho_t(i+1,j,k+1))) - r
1338  dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j,k+1)* p_t(i+1,j,k+1))) - p
1339  dry = qrt(i,j,k,ir,2) - r
1340  dpy = qrt(i,j,k,ip,2) - p
1341  qrt(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1342  qrt(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1343  ! Y-edge averaged right-bottom corner state (RB->LR)
1344  !qRB(i,j,k,ir,2) = r + (+drx-drz)
1345  drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j,k-1)*rho_t(i+1,j,k-1))) - r
1346  dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j,k-1)* p_t(i+1,j,k-1))) - p
1347  dry = qrb(i,j,k,ir,2) - r
1348  dpy = qrb(i,j,k,ip,2) - p
1349  qrb(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1350  qrb(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1351 
1352  ! Y-edge averaged left-top corner state (LT->RL)
1353  !qLT(i,j,k,ir,2) = r + (-drx+drz)
1354  drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j,k+1)*rho_t(i-1,j,k+1))) - r
1355  dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j,k+1)* p_t(i-1,j,k+1))) - p
1356  dry = qlt(i,j,k,ir,2) - r
1357  dpy = qlt(i,j,k,ip,2) - p
1358  qlt(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1359  qlt(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1360  ! Y-edge averaged left-bottom corner state (LB->RR)
1361  !qLB(i,j,k,ir,2) = r + (-drx-drz)
1362  drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j,k-1)*rho_t(i-1,j,k-1))) - r
1363  dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j,k-1)* p_t(i-1,j,k-1))) - p
1364  dry = qlb(i,j,k,ir,2) - r
1365  dpy = qlb(i,j,k,ip,2) - p
1366  qlb(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1367  qlb(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1368  ! Z-edge averaged right-top corner state (RT->LL)
1369  !qRT(i,j,k,ir,3) = r + (+drx+dry)
1370  drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j+1,k)*rho_t(i+1,j+1,k))) - r
1371  dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j+1,k)* p_t(i+1,j+1,k))) - p
1372  dry = qrt(i,j,k,ir,3) - r
1373  dpy = qrt(i,j,k,ip,3) - p
1374  qrt(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1375  qrt(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1376  ! Z-edge averaged right-bottom corner state (RB->LR)
1377  !qRB(i,j,k,ir,3) = r + (+drx-dry)
1378  drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j-1,k)*rho_t(i+1,j-1,k))) - r
1379  dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j-1,k)* p_t(i+1,j-1,k))) - p
1380  dry = qrb(i,j,k,ir,3) - r
1381  dpy = qrb(i,j,k,ip,3) - p
1382  qrb(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1383  qrb(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1384  ! Z-edge averaged left-top corner state (LT->RL)
1385  !qLT(i,j,k,ir,3) = r + (-drx+dry)
1386  drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j+1,k)*rho_t(i-1,j+1,k))) - r
1387  dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j+1,k)* p_t(i-1,j+1,k))) - p
1388  dry = qlt(i,j,k,ir,3) - r
1389  dpy = qlt(i,j,k,ip,3) - p
1390  qlt(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1391  qlt(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1392  ! Z-edge averaged left-bottom corner state (LB->RR)
1393  !qLB(i,j,k,ir,3) = r + (-drx-dry)
1394  drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j-1,k)*rho_t(i-1,j-1,k))) - r
1395  dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j-1,k)* p_t(i-1,j-1,k))) - p
1396  dry = qlb(i,j,k,ir,3) - r
1397  dpy = qlb(i,j,k,ip,3) - p
1398  qlb(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1399  qlb(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1400 
1401 ! end associate
1402  end do
1403  end do
1404  end do
1405  !call trace%end (itimer(3))
1406  call trace%end (itimer(1), detailed_timer=detailed_timer)
1407  end associate
1408 contains
1409  FUNCTION xdn(a) RESULT (b)
1410  real, dimension(:,:,:):: a
1411  real, dimension(size(a,1),size(a,2),size(a,3)):: b
1412  do k=1,size(a,3)
1413  do j=1,size(a,2)
1414  do i=2,size(a,1)
1415  b(i,j,k)=half*(a(i-1,j,k)+a(i,j,k))
1416  end do
1417  end do
1418  end do
1419  b(1,:,:) = a(1,:,:)
1420  END FUNCTION
1421  FUNCTION ydn(a) RESULT (b)
1422  real, dimension(:,:,:):: a
1423  real, dimension(size(a,1),size(a,2),size(a,3)):: b
1424  do k=1,size(a,3)
1425  do j=2,size(a,2)
1426  do i=1,size(a,1)
1427  b(i,j,k)=half*(a(i,j-1,k)+a(i,j,k))
1428  end do
1429  end do
1430  end do
1431  b(:,1,:) = a(:,1,:)
1432  END FUNCTION
1433  FUNCTION zdn(a) RESULT (b)
1434  real, dimension(:,:,:):: a
1435  real, dimension(size(a,1),size(a,2),size(a,3)):: b
1436  do k=2,size(a,3)
1437  do j=1,size(a,2)
1438  do i=1,size(a,1)
1439  b(i,j,k)=half*(a(i,j,k-1)+a(i,j,k))
1440  end do
1441  end do
1442  end do
1443  b(:,:,1) = a(:,:,1)
1444  END FUNCTION
1445 END SUBROUTINE trace3d
1446 
1447 !===============================================================================
1448 !> Compute fluxes
1449 !===============================================================================
1450 subroutine cmpflxm(self, nx, ln, lt1, lt2, bn, bt1, bt2, fl)
1451  class(mhd_t):: self
1452  real(dp), dimension(:,:,:,:):: fl
1453  integer:: nx, ln, lt1, lt2, bn, bt1, bt2
1454  ! local variables
1455  integer:: i, j, k, lb(3), ub(3), xdim, l(3), idim
1456  real(dp), dimension(nx,self%nv):: qleft, qright, fgdnv
1457  real(dp), dimension(nx):: bn_mean
1458  real(dp):: dtds(3)
1459  integer, save:: itimer=0
1460  integer:: n
1461  associate(qp=>self%qp, qm=>self%qm)
1462  !-----------------------------------------------------------------------------
1463  call trace%begin ('mhd_t%cmpflxm', itimer=itimer, detailed_timer=detailed_timer)
1464  lb = self%mesh%lb+2
1465  ub = self%mesh%ub-1
1466  xdim = ln-1
1467  if (verbose>0) print '(a,2(2x,3i4))', 'mhd_t%cmpflxm: lb, ub=', lb, ub
1468  !
1469  ! Avoid using spurious values from allocated, but non-defined entries
1470  bn_mean=0d0
1471 
1472  if (xdim == 1) then
1473  do k = lb(3),ub(3)
1474  do j = lb(2),ub(2)
1475  ! Enforce continuity for normal magnetic field
1476  n = ub(1)-lb(1)+1
1477  bn_mean(1:n) = half*(qm(lb(1)-1:ub(1)-1,j,k,bn,xdim)+qp(lb(1):ub(1),j,k,bn,xdim))
1478  ! Left state
1479  qleft(1:n,1) = qm(lb(1)-1:ub(1)-1,j,k,1 ,xdim) ! Mass density
1480  qleft(1:n,2) = qm(lb(1)-1:ub(1)-1,j,k,5 ,xdim) ! Pressure
1481  qleft(1:n,3) = qm(lb(1)-1:ub(1)-1,j,k,ln ,xdim) ! Normal velocity
1482  qleft(1:n,4) = bn_mean(1:n) ! Normal magnetic field
1483  qleft(1:n,5) = qm(lb(1)-1:ub(1)-1,j,k,lt1,xdim) ! Tangential velocity 1
1484  qleft(1:n,6) = qm(lb(1)-1:ub(1)-1,j,k,bt1,xdim) ! Tangential magnetic field 1
1485  qleft(1:n,7) = qm(lb(1)-1:ub(1)-1,j,k,lt2,xdim) ! Tangential velocity 2
1486  qleft(1:n,8) = qm(lb(1)-1:ub(1)-1,j,k,bt2,xdim) ! Tangential magnetic field 2
1487  ! Right state
1488  qright(1:n,1) = qp(lb(1):ub(1),j,k,1 ,xdim) ! Mass density
1489  qright(1:n,2) = qp(lb(1):ub(1),j,k,5 ,xdim) ! Pressure
1490  qright(1:n,3) = qp(lb(1):ub(1),j,k,ln ,xdim) ! Normal velocity
1491  qright(1:n,4) = bn_mean(1:n) ! Normal magnetic field
1492  qright(1:n,5) = qp(lb(1):ub(1),j,k,lt1,xdim) ! Tangential velocity 1
1493  qright(1:n,6) = qp(lb(1):ub(1),j,k,bt1,xdim) ! Tangential magnetic field 1
1494  qright(1:n,7) = qp(lb(1):ub(1),j,k,lt2,xdim) ! Tangential velocity 2
1495  qright(1:n,8) = qp(lb(1):ub(1),j,k,bt2,xdim) ! Tangential magnetic field 2
1496  call solver (self,qleft,qright,fgdnv,n)
1497  ! Output fluxes
1498  fl(lb(1):ub(1),j,k,1 ) = fgdnv(1:n,1) ! Mass density
1499  fl(lb(1):ub(1),j,k,5 ) = fgdnv(1:n,2) ! Total energy
1500  fl(lb(1):ub(1),j,k,ln ) = fgdnv(1:n,3) ! Normal momentum
1501  fl(lb(1):ub(1),j,k,bn ) = fgdnv(1:n,4) ! Normal magnetic field
1502  fl(lb(1):ub(1),j,k,lt1) = fgdnv(1:n,5) ! Transverse momentum 1
1503  fl(lb(1):ub(1),j,k,bt1) = fgdnv(1:n,6) ! Transverse magnetic field 1
1504  fl(lb(1):ub(1),j,k,lt2) = fgdnv(1:n,7) ! Transverse momentum 2
1505  fl(lb(1):ub(1),j,k,bt2) = fgdnv(1:n,8) ! Transverse magnetic field 2
1506  end do
1507  end do
1508  else if (xdim == 2) then
1509  do i = lb(1),ub(1)
1510  do k = lb(3),ub(3)
1511  ! Enforce continuity for normal magnetic field
1512  n = ub(2)-lb(2)+1
1513  bn_mean(1:n) = half*(qm(i,lb(2)-1:ub(2)-1,k,bn,xdim)+qp(i,lb(2):ub(2),k,bn,xdim))
1514  ! Left state
1515  qleft(1:n,1) = qm(i,lb(2)-1:ub(2)-1,k,1 ,xdim) ! Mass density
1516  qleft(1:n,2) = qm(i,lb(2)-1:ub(2)-1,k,5 ,xdim) ! Pressure
1517  qleft(1:n,3) = qm(i,lb(2)-1:ub(2)-1,k,ln ,xdim) ! Normal velocity
1518  qleft(1:n,4) = bn_mean(1:n) ! Normal magnetic field
1519  qleft(1:n,5) = qm(i,lb(2)-1:ub(2)-1,k,lt1,xdim) ! Tangential velocity 1
1520  qleft(1:n,6) = qm(i,lb(2)-1:ub(2)-1,k,bt1,xdim) ! Tangential magnetic field 1
1521  qleft(1:n,7) = qm(i,lb(2)-1:ub(2)-1,k,lt2,xdim) ! Tangential velocity 2
1522  qleft(1:n,8) = qm(i,lb(2)-1:ub(2)-1,k,bt2,xdim) ! Tangential magnetic field 2
1523  ! Right state
1524  qright(1:n,1) = qp(i,lb(2):ub(2),k,1 ,xdim) ! Mass density
1525  qright(1:n,2) = qp(i,lb(2):ub(2),k,5 ,xdim) ! Pressure
1526  qright(1:n,3) = qp(i,lb(2):ub(2),k,ln ,xdim) ! Normal velocity
1527  qright(1:n,4) = bn_mean(1:n) ! Normal magnetic field
1528  qright(1:n,5) = qp(i,lb(2):ub(2),k,lt1,xdim) ! Tangential velocity 1
1529  qright(1:n,6) = qp(i,lb(2):ub(2),k,bt1,xdim) ! Tangential magnetic field 1
1530  qright(1:n,7) = qp(i,lb(2):ub(2),k,lt2,xdim) ! Tangential velocity 2
1531  qright(1:n,8) = qp(i,lb(2):ub(2),k,bt2,xdim) ! Tangential magnetic field 2
1532  call solver (self,qleft,qright,fgdnv,n)
1533  ! Output fluxes
1534  fl(i,lb(2):ub(2),k,1 ) = fgdnv(1:n,1) ! Mass density
1535  fl(i,lb(2):ub(2),k,5 ) = fgdnv(1:n,2) ! Total energy
1536  fl(i,lb(2):ub(2),k,ln ) = fgdnv(1:n,3) ! Normal momentum
1537  fl(i,lb(2):ub(2),k,bn ) = fgdnv(1:n,4) ! Normal magnetic field
1538  fl(i,lb(2):ub(2),k,lt1) = fgdnv(1:n,5) ! Transverse momentum 1
1539  fl(i,lb(2):ub(2),k,bt1) = fgdnv(1:n,6) ! Transverse magnetic field 1
1540  fl(i,lb(2):ub(2),k,lt2) = fgdnv(1:n,7) ! Transverse momentum 2
1541  fl(i,lb(2):ub(2),k,bt2) = fgdnv(1:n,8) ! Transverse magnetic field 2
1542  end do
1543  end do
1544  else if (xdim == 3) then
1545  do j = lb(2),ub(2)
1546  do i = lb(1),ub(1)
1547  ! Enforce continuity for normal magnetic field
1548  n = ub(3)-lb(3)+1
1549  bn_mean(1:n) = half*(qm(i,j,lb(3)-1:ub(3)-1,bn,xdim)+qp(i,j,lb(3):ub(3),bn,xdim))
1550  ! Left state
1551  qleft(1:n,1) = qm(i,j,lb(3)-1:ub(3)-1,1 ,xdim) ! Mass density
1552  qleft(1:n,2) = qm(i,j,lb(3)-1:ub(3)-1,5 ,xdim) ! Pressure
1553  qleft(1:n,3) = qm(i,j,lb(3)-1:ub(3)-1,ln ,xdim) ! Normal velocity
1554  qleft(1:n,4) = bn_mean(1:n) ! Normal magnetic field
1555  qleft(1:n,5) = qm(i,j,lb(3)-1:ub(3)-1,lt1,xdim) ! Tangential velocity 1
1556  qleft(1:n,6) = qm(i,j,lb(3)-1:ub(3)-1,bt1,xdim) ! Tangential magnetic field 1
1557  qleft(1:n,7) = qm(i,j,lb(3)-1:ub(3)-1,lt2,xdim) ! Tangential velocity 2
1558  qleft(1:n,8) = qm(i,j,lb(3)-1:ub(3)-1,bt2,xdim) ! Tangential magnetic field 2
1559  ! Right state
1560  qright(1:n,1) = qp(i,j,lb(3):ub(3),1 ,xdim) ! Mass density
1561  qright(1:n,2) = qp(i,j,lb(3):ub(3),5 ,xdim) ! Pressure
1562  qright(1:n,3) = qp(i,j,lb(3):ub(3),ln ,xdim) ! Normal velocity
1563  qright(1:n,4) = bn_mean(1:n) ! Normal magnetic field
1564  qright(1:n,5) = qp(i,j,lb(3):ub(3),lt1,xdim) ! Tangential velocity 1
1565  qright(1:n,6) = qp(i,j,lb(3):ub(3),bt1,xdim) ! Tangential magnetic field 1
1566  qright(1:n,7) = qp(i,j,lb(3):ub(3),lt2,xdim) ! Tangential velocity 2
1567  qright(1:n,8) = qp(i,j,lb(3):ub(3),bt2,xdim) ! Tangential magnetic field 2
1568  call solver (self,qleft,qright,fgdnv,n)
1569  ! Output fluxes
1570  fl(i,j,lb(3):ub(3),1 ) = fgdnv(1:n,1) ! Mass density
1571  fl(i,j,lb(3):ub(3),5 ) = fgdnv(1:n,2) ! Total energy
1572  fl(i,j,lb(3):ub(3),ln ) = fgdnv(1:n,3) ! Normal momentum
1573  fl(i,j,lb(3):ub(3),bn ) = fgdnv(1:n,4) ! Normal magnetic field
1574  fl(i,j,lb(3):ub(3),lt1) = fgdnv(1:n,5) ! Transverse momentum 1
1575  fl(i,j,lb(3):ub(3),bt1) = fgdnv(1:n,6) ! Transverse magnetic field 1
1576  fl(i,j,lb(3):ub(3),lt2) = fgdnv(1:n,7) ! Transverse momentum 2
1577  fl(i,j,lb(3):ub(3),bt2) = fgdnv(1:n,8) ! Transverse magnetic field 2
1578  end do
1579  end do
1580 
1581  endif
1582 
1583  end associate
1584  call trace%end (itimer, detailed_timer=detailed_timer)
1585 end subroutine cmpflxm
1586 
1587 !===============================================================================
1588 !> 2D solver for EMF
1589 !===============================================================================
1590 SUBROUTINE cmp_mag_flx (self, qRT, qRB, qLT, qLB, lp1 ,lp2 ,lor ,bp1 , bp2, bor, emf, ngrid, nvar)
1591  class(mhd_t):: self
1592  ! 2D Riemann solver to compute EMF at cell edges
1593  ! indices of the 2 planar velocity lp1 and lp2 and the orthogonal one,
1594  ! lor and idem for the magnetic field
1595  integer:: lp1, lp2, lor, bp1, bp2, bor, ngrid, nvar
1596  integer:: irt, jrt, krt, irb, jrb, krb, ilt, jlt, klt, ilb, jlb, klb
1597  real(dp),DIMENSION(:,:,:):: emf
1598  real(dp),DIMENSION(:,:,:,:,:):: qrt, qlt, qrb, qlb
1599  ! local variables
1600  integer:: i, j, k, n, xdim, l, u
1601  real(dp),DIMENSION(1:ngrid,nvar)::qll,qrl,qlr,qrr
1602  real(dp):: smallp
1603  integer, save:: itimer=0
1604  !associate (qRT=>self%qRT,qLT=>self%qLT,qRB=>self%qRB,qLB=>self%qLB)
1605  call trace%begin ('mhd_t%cmp_mag_flx', itimer=itimer, detailed_timer=detailed_timer)
1606  !-----------------------------------------------------------------------------
1607  xdim = lor - 1
1608 ! if (xdim == 3 .or. xdim == 1) then
1609 ! associate (qRT=>self%qRT,qLT=>self%qLT,qRB=>self%qRB,qLB=>self%qLB)
1610 ! else if (xdim == 2) then
1611 ! associate (qRT=>self%qRT,qLT=>self%qRB,qRB=>self%qLT,qLB=>self%qLB)
1612 ! endif
1613  !
1614  smallp = smallr*smallc**2
1615  !
1616  qll = 0d0
1617  qrl = 0d0
1618  qlr = 0d0
1619  qrr = 0d0
1620  !
1621  if (xdim == 3) then
1622  irt = 1; jrt = 1; krt = 0
1623  irb = 1; jrb = 0; krb = 0
1624  ilt = 0; jlt = 1; klt = 0
1625  ilb = 0; jlb = 0; klb = 0
1626  else if (xdim == 2) then
1627  irt = 1; jrt = 0; krt = 1
1628  irb = 0; jrb = 0; krb = 1
1629  ilt = 1; jlt = 0; klt = 0
1630  ilb = 0; jlb = 0; klb = 0
1631  else if (xdim == 1) then
1632  irt = 0; jrt = 1; krt = 1
1633  irb = 0; jrb = 1; krb = 0
1634  ilt = 0; jlt = 0; klt = 1
1635  ilb = 0; jlb = 0; klb = 0
1636  endif
1637  !
1638  if (xdim == 1) then
1639  do k = self%mesh(3)%li,self%mesh(3)%uo
1640  do j = self%mesh(2)%li,self%mesh(2)%uo
1641  l = self%mesh(1)%li
1642  u = self%mesh(1)%ui
1643  do i = l,u
1644  qll(i,1) = max(qrt(i-irt, j-jrt, k-krt, 1, xdim), smallr)
1645  qrl(i,1) = max(qlt(i-ilt, j-jlt, k-klt, 1, xdim), smallr)
1646  qlr(i,1) = max(qrb(i-irb, j-jrb, k-krb, 1, xdim), smallr)
1647  qrr(i,1) = max(qlb(i-ilb, j-jlb, k-klb, 1, xdim), smallr)
1648  qll(i,2) = max(qrt(i-irt, j-jrt, k-krt, 5, xdim), smallp)
1649  qrl(i,2) = max(qlt(i-ilt, j-jlt, k-klt, 5, xdim), smallp)
1650  qlr(i,2) = max(qrb(i-irb, j-jrb, k-krb, 5, xdim), smallp)
1651  qrr(i,2) = max(qlb(i-ilb, j-jlb, k-klb, 5, xdim), smallp)
1652  qll(i,3) = qrt(i-irt, j-jrt, k-krt, lp1, xdim)
1653  qrl(i,3) = qlt(i-ilt, j-jlt, k-klt, lp1, xdim)
1654  qlr(i,3) = qrb(i-irb, j-jrb, k-krb, lp1, xdim)
1655  qrr(i,3) = qlb(i-ilb, j-jlb, k-klb, lp1, xdim)
1656  qll(i,4) = qrt(i-irt, j-jrt, k-krt, lp2, xdim)
1657  qrl(i,4) = qlt(i-ilt, j-jlt, k-klt, lp2, xdim)
1658  qlr(i,4) = qrb(i-irb, j-jrb, k-krb, lp2, xdim)
1659  qrr(i,4) = qlb(i-ilb, j-jlb, k-klb, lp2, xdim)
1660  qll(i,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1661  qrl(i,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1662  qlr(i,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1663  qrr(i,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1664  qll(i,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1665  qrl(i,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1666  qlr(i,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1667  qrr(i,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1668  qll(i,5) = qrt(i-irt, j-jrt, k-krt, lor, xdim)
1669  qrl(i,5) = qlt(i-ilt, j-jlt, k-klt, lor, xdim)
1670  qlr(i,5) = qrb(i-irb, j-jrb, k-krb, lor, xdim)
1671  qrr(i,5) = qlb(i-ilb, j-jlb, k-klb, lor, xdim)
1672  qll(i,8) = qrt(i-irt, j-jrt, k-krt, bor, xdim)
1673  qrl(i,8) = qlt(i-ilt, j-jlt, k-klt, bor, xdim)
1674  qlr(i,8) = qrb(i-irb, j-jrb, k-krb, bor, xdim)
1675  qrr(i,8) = qlb(i-ilb, j-jlb, k-klb, bor, xdim)
1676  end do
1677  call solver_2d (self,qll,qrl,qlr,qrr,emf,j,k,l,u,nvar,xdim,detailed_timer=detailed_timer)
1678  end do
1679  end do
1680  else if (xdim == 2) then
1681  do i = self%mesh(1)%li,self%mesh(1)%uo
1682  do k = self%mesh(3)%li,self%mesh(3)%uo
1683  l = self%mesh(2)%li
1684  u = self%mesh(2)%ui
1685  do j = l,u
1686  qll(j,1) = max(qrt(i-irt, j-jrt, k-krt, 1, xdim), smallr)
1687  qrl(j,1) = max(qlt(i-ilt, j-jlt, k-klt, 1, xdim), smallr)
1688  qlr(j,1) = max(qrb(i-irb, j-jrb, k-krb, 1, xdim), smallr)
1689  qrr(j,1) = max(qlb(i-ilb, j-jlb, k-klb, 1, xdim), smallr)
1690  qll(j,2) = max(qrt(i-irt, j-jrt, k-krt, 5, xdim), smallp)
1691  qrl(j,2) = max(qlt(i-ilt, j-jlt, k-klt, 5, xdim), smallp)
1692  qlr(j,2) = max(qrb(i-irb, j-jrb, k-krb, 5, xdim), smallp)
1693  qrr(j,2) = max(qlb(i-ilb, j-jlb, k-klb, 5, xdim), smallp)
1694  qll(j,3) = qrt(i-irt, j-jrt, k-krt, lp1, xdim)
1695  qrl(j,3) = qlt(i-ilt, j-jlt, k-klt, lp1, xdim)
1696  qlr(j,3) = qrb(i-irb, j-jrb, k-krb, lp1, xdim)
1697  qrr(j,3) = qlb(i-ilb, j-jlb, k-klb, lp1, xdim)
1698  qll(j,4) = qrt(i-irt, j-jrt, k-krt, lp2, xdim)
1699  qrl(j,4) = qlt(i-ilt, j-jlt, k-klt, lp2, xdim)
1700  qlr(j,4) = qrb(i-irb, j-jrb, k-krb, lp2, xdim)
1701  qrr(j,4) = qlb(i-ilb, j-jlb, k-klb, lp2, xdim)
1702  qll(j,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1703  qrl(j,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1704  qlr(j,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1705  qrr(j,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1706  qll(j,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1707  qrl(j,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1708  qlr(j,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1709  qrr(j,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1710  qll(j,5) = qrt(i-irt, j-jrt, k-krt, lor, xdim)
1711  qrl(j,5) = qlt(i-ilt, j-jlt, k-klt, lor, xdim)
1712  qlr(j,5) = qrb(i-irb, j-jrb, k-krb, lor, xdim)
1713  qrr(j,5) = qlb(i-ilb, j-jlb, k-klb, lor, xdim)
1714  qll(j,8) = qrt(i-irt, j-jrt, k-krt, bor, xdim)
1715  qrl(j,8) = qlt(i-ilt, j-jlt, k-klt, bor, xdim)
1716  qlr(j,8) = qrb(i-irb, j-jrb, k-krb, bor, xdim)
1717  qrr(j,8) = qlb(i-ilb, j-jlb, k-klb, bor, xdim)
1718  end do
1719  call solver_2d (self,qll,qrl,qlr,qrr,emf,k,i,l,u,nvar,xdim,detailed_timer=detailed_timer)
1720  end do
1721  end do
1722  else if (xdim == 3) then
1723  do j = self%mesh(2)%li,self%mesh(2)%uo
1724  do i = self%mesh(1)%li,self%mesh(1)%uo
1725  l = self%mesh(3)%li
1726  u = self%mesh(3)%ui
1727  do k = l,u
1728  qll(k,1) = max(qrt(i-irt, j-jrt, k-krt, 1, xdim), smallr)
1729  qrl(k,1) = max(qlt(i-ilt, j-jlt, k-klt, 1, xdim), smallr)
1730  qlr(k,1) = max(qrb(i-irb, j-jrb, k-krb, 1, xdim), smallr)
1731  qrr(k,1) = max(qlb(i-ilb, j-jlb, k-klb, 1, xdim), smallr)
1732  qll(k,2) = max(qrt(i-irt, j-jrt, k-krt, 5, xdim), smallp)
1733  qrl(k,2) = max(qlt(i-ilt, j-jlt, k-klt, 5, xdim), smallp)
1734  qlr(k,2) = max(qrb(i-irb, j-jrb, k-krb, 5, xdim), smallp)
1735  qrr(k,2) = max(qlb(i-ilb, j-jlb, k-klb, 5, xdim), smallp)
1736  qll(k,3) = qrt(i-irt, j-jrt, k-krt, lp1, xdim)
1737  qrl(k,3) = qlt(i-ilt, j-jlt, k-klt, lp1, xdim)
1738  qlr(k,3) = qrb(i-irb, j-jrb, k-krb, lp1, xdim)
1739  qrr(k,3) = qlb(i-ilb, j-jlb, k-klb, lp1, xdim)
1740  qll(k,4) = qrt(i-irt, j-jrt, k-krt, lp2, xdim)
1741  qrl(k,4) = qlt(i-ilt, j-jlt, k-klt, lp2, xdim)
1742  qlr(k,4) = qrb(i-irb, j-jrb, k-krb, lp2, xdim)
1743  qrr(k,4) = qlb(i-ilb, j-jlb, k-klb, lp2, xdim)
1744  qll(k,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1745  qrl(k,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1746  qlr(k,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1747  qrr(k,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1748  qll(k,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1749  qrl(k,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1750  qlr(k,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1751  qrr(k,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1752  qll(k,5) = qrt(i-irt, j-jrt, k-krt, lor, xdim)
1753  qrl(k,5) = qlt(i-ilt, j-jlt, k-klt, lor, xdim)
1754  qlr(k,5) = qrb(i-irb, j-jrb, k-krb, lor, xdim)
1755  qrr(k,5) = qlb(i-ilb, j-jlb, k-klb, lor, xdim)
1756  qll(k,8) = qrt(i-irt, j-jrt, k-krt, bor, xdim)
1757  qrl(k,8) = qlt(i-ilt, j-jlt, k-klt, bor, xdim)
1758  qlr(k,8) = qrb(i-irb, j-jrb, k-krb, bor, xdim)
1759  qrr(k,8) = qlb(i-ilb, j-jlb, k-klb, bor, xdim)
1760  end do
1761  call solver_2d (self,qll,qrl,qlr,qrr,emf,i,j,l,u,nvar,xdim,detailed_timer=detailed_timer)
1762  end do
1763  end do
1764  endif
1765  call trace%end (itimer, detailed_timer=detailed_timer)
1766 ! end associate
1767 END SUBROUTINE cmp_mag_flx
1768 
1769 !===============================================================================
1770 !> Selection of 1-D solvers
1771 !===============================================================================
1772 SUBROUTINE solver (self, qleft, qright, fgdnv, ub)
1773  class(mhd_t):: self
1774  real(dp),dimension(:,:):: qleft,qright,fgdnv
1775  integer :: lb, ub
1776  !.............................................................................
1777  if (self%riemann == 'hlld') then
1778  call rieman%hlld (qleft, qright, fgdnv, ub)
1779  else if (self%riemann == 'llf')then
1780  call rieman%llf (qleft, qright, fgdnv, ub)
1781  else
1782  write(*,*)'unknown Riemann solver'
1783  stop
1784  end if
1785 END SUBROUTINE solver
1786 
1787 !===============================================================================
1788 !> Selection of 2-D solvers
1789 !===============================================================================
1790 SUBROUTINE solver_2d (self,qLL,qRL,qLR,qRR,emf,j,k,l,u,nvar,xdim,detailed_timer)
1791  class(mhd_t):: self
1792  real(DP),dimension(:,:),intent(in):: qll,qrl,qlr,qrr
1793  real(DP),dimension(:,:,:):: emf
1794  integer, intent(in):: j,k,l,u,nvar,xdim
1795  logical, optional:: detailed_timer
1796  !.............................................................................
1797  if (self%riemann == 'hlld') then
1798  call rieman%hlld_2d (qll,qrl,qlr,qrr,emf,j,k,l,u,nvar,xdim,detailed_timer=detailed_timer)
1799  else if (self%riemann == 'llf')then
1800  call rieman%llf_2d (qll,qrl,qlr,qrr,emf,j,k,l,u,nvar,xdim,detailed_timer=detailed_timer)
1801  else
1802  write(*,*)'unknown Riemann 2D solver'
1803  stop
1804  end if
1805 END SUBROUTINE solver_2d
1806 
1807 !===============================================================================
1808 FUNCTION up(f,i) RESULT (g)
1809  real, dimension(:,:,:), intent(in):: f
1810  real, dimension(size(f,1),size(f,2),size(f,3)):: g
1811  integer:: i
1812  !.............................................................................
1813  g = 0.5*(cshift(f,1,i)+f)
1814 END FUNCTION
1815 
1816 !===============================================================================
1817 FUNCTION gas_pressure (self) RESULT (pg)
1818  class(mhd_t):: self
1819  real, dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
1820  integer :: l(3), u(3)
1821  !-----------------------------------------------------------------------------
1822  associate(d => self%mem(:,:,:,self%idx%d ,self%it,1), &
1823  px => self%mem(:,:,:,self%idx%px,self%it,1), &
1824  py => self%mem(:,:,:,self%idx%py,self%it,1), &
1825  pz => self%mem(:,:,:,self%idx%pz,self%it,1), &
1826  bx => self%mem(:,:,:,self%idx%bx,self%it,1), &
1827  by => self%mem(:,:,:,self%idx%by,self%it,1), &
1828  bz => self%mem(:,:,:,self%idx%bz,self%it,1), &
1829  e => self%mem(:,:,:,self%idx%e ,self%it,1))
1830  if (isothermal) then
1831  pg = d*self%csound**2
1832  else
1833  pg = (gamma-1d0)* &
1834  (e-0.5*(up(bx,1)**2+up(by,2)**2+up(bz,3)**2+(px**2+py**2+pz**2)/d))
1835  end if
1836  end associate
1837 END FUNCTION gas_pressure
1838 
1839 !===============================================================================
1840 !> Add gradient of phi gravitational force. FIXME: This should be in extras
1841 !===============================================================================
1842 SUBROUTINE gravitational_force (self)
1843  class(mhd_t):: self
1844  real, dimension(:,:,:), pointer:: phi, d, et, gradphi1, gradphi2, gradphi3
1845  real, dimension(:,:,:,:), pointer:: p
1846  integer:: i,j,k
1847  real, dimension(3):: wtp1, wtm1
1848  real:: etemp
1849  !.............................................................................
1850  associate(l => self%li, u => self%ui)
1851  wtp1 = 0.50 * 4.0 / 3.0 / self%ds
1852  wtm1 = 0.25 * 1.0 / 3.0 / self%ds
1853  phi => self%mem(:,:,:,self%idx%phi,self%it,1)
1854  d => self%mem(:,:,:,self%idx%d,self%it,1)
1855  et => self%mem(:,:,:,self%idx%e,self%new,1)
1856  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
1857  allocate (gradphi1(self%gn(1),self%gn(2),self%gn(3)), &
1858  gradphi2(self%gn(1),self%gn(2),self%gn(3)), &
1859  gradphi3(self%gn(1),self%gn(2),self%gn(3)))
1860 
1861  do k=l(3),u(3)
1862  do j=l(2),u(2)
1863  do i=l(1),u(1)
1864  ! calculate energy minus the kinetic energy
1865  etemp = et(i,j,k) - 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1866  - 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1867  - 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1868  ! calculate gradient of phi
1869  gradphi1(i,j,k) = wtp1(1) * (phi(i+1,j,k) - phi(i-1,j,k)) &
1870  - wtm1(1) * (phi(i+2,j,k) - phi(i-2,j,k))
1871  gradphi2(i,j,k) = wtp1(2) * (phi(i,j+1,k) - phi(i,j-1,k)) &
1872  - wtm1(2) * (phi(i,j+2,k) - phi(i,j-2,k))
1873  gradphi3(i,j,k) = wtp1(3) * (phi(i,j,k+1) - phi(i,j,k-1)) &
1874  - wtm1(3) * (phi(i,j,k+2) - phi(i,j,k-2))
1875  ! update momentum
1876  p(i,j,k,1) = p(i,j,k,1) - self%dtime * d(i,j,k) * gradphi1(i,j,k)
1877  p(i,j,k,2) = p(i,j,k,2) - self%dtime * d(i,j,k) * gradphi2(i,j,k)
1878  p(i,j,k,3) = p(i,j,k,3) - self%dtime * d(i,j,k) * gradphi3(i,j,k)
1879  ! update total energy with new momentum
1880  et(i,j,k) = etemp + 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1881  + 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1882  + 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1883  end do
1884  end do
1885  end do
1886 
1887  deallocate (gradphi1, gradphi2, gradphi3)
1888  end associate
1889 END SUBROUTINE gravitational_force
1890 
1891 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
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
non-ideal MHD module
RAMSES Godunov solvers, use of guard zones; specifically in HLLD.
Definition: mhd_mod.f90:16
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...