DISPATCH
mhd_mod.f90
1 !===============================================================================
2 !> Staggered mesh solver using entropy per unit volume as the energy variable.
3 !>
4 !> Flop counting is enable via annotation at the end of lines. To count flops,
5 !> use a Linux pipeline such as
6 !>
7 !> % grep ops: ../../../solvers/stagger2/mhd_mod.f90 | sed 's/.*ops://' | sumcol 2
8 !>
9 !> where 'sumcol' is any script that sums up values in column 2
10 !===============================================================================
11 MODULE mhd_mod
12  USE io_mod
13  USE bits_mod
14  USE omp_mod
15  USE omp_timer_mod
16  USE mpi_mod
17  USE trace_mod
18  USE kinds_mod
19  USE units_mod
20  USE scaling_mod
21  USE scalar_mod
22  USE vector_mod
23  USE vector_ops
24  USE stagger_mod
25  USE extras_mod
26  USE index_mod
27  USE timestep_mod
28  USE eos_mod
29  USE non_ideal_mod
30  USE initial_mod
31  implicit none
32  private
33  type, extends(extras_t), public:: mhd_t
34  real, dimension(:,:,:), pointer:: d, dddt, s, dsdt, pg
35  real, dimension(:,:,:,:), pointer:: p, dpdt, b, dbdt
36  real:: nu(6)
37  real:: cdtd
38  logical:: mhd=.false.
39  logical:: first_time=.true.
40  logical:: do_force=.true.
41  type(extras_t):: extras
42  contains
43  procedure:: init
44  procedure:: test
45  procedure:: pde
46  procedure:: update
47  procedure:: output
48  procedure:: gas_pressure
49  procedure:: gas_pressure_sub
50  end type
51  integer:: verbose=0
52  logical, save:: flop_count=.false.
53  logical, save:: do_2nd_div=.true.
54  logical, save:: do_smooth=.false.
55  logical, save:: do_maxwell=.true.
56 CONTAINS
57 
58 !===============================================================================
59 !> The construction with a double test on first_time avoids entering a critical
60 !> region once the variable has been set.
61 !===============================================================================
62 SUBROUTINE init (self)
63  class(mhd_t):: self
64  integer:: i, iv
65  real, save:: nu(6)=[0.1,1.0,0.0,0.5,0.5,0.5], csound=1.0, cdtd=0.5, courant=0.2
66  real(8), save:: gamma=1.4_8
67  logical, save:: first_time=.true.
68  character(len=16), save:: eos
69  namelist /stagger_params/ nu, courant, gamma, csound, cdtd, hardwire, eos, &
70  do_smooth, do_2nd_div, flop_count, verbose, do_maxwell
71  !----------------------------------------------------------------------------
72  call trace%begin('mhd_t%init')
73  !$omp critical (input_cr)
74  if (first_time) then
75  eos = self%eos
76  rewind(io%input)
77  read (io%input, stagger_params)
78  if (io%master) write (*, stagger_params)
79  call non_ideal%init
80  call timestep%init
81  ! fool-proofing
82  if (timestep%time_order <= 0) call mpi%abort("The Stagger solvers require `time_order` > 0. Abort!")
83  end if
84  !$omp end critical (input_cr)
85  self%eos = eos
86  self%nu = nu
87  self%gamma = gamma
88  self%csound = csound
89  self%courant = courant
90  self%cdtd = cdtd
91  !-----------------------------------------------------------------------------
92  ! Read IC parameters, which determine if this is HD or MHD. The IC init is
93  ! called again, from extras_t%init, when gamma is known and can be passed on.
94  !-----------------------------------------------------------------------------
95  call self%idx%init (5, self%mhd)
96  self%initial%mhd = self%mhd
97  call self%initial%init (self%kind)
98  self%mhd = self%initial%mhd
99  if (self%mhd) then
100  self%kind = 'stagger2_mhd_patch'
101  else
102  self%kind = 'stagger2_hd_patch'
103  self%idx%bx = -1
104  self%idx%by = -1
105  self%idx%bz = -1
106  end if
107  if (self%nv==0) then
108  if (self%mhd) then
109  self%nv = 8
110  else
111  self%nv = 5
112  end if
113  end if
114  !-----------------------------------------------------------------------------
115  ! Allocate memory and mesh, etc
116  !-----------------------------------------------------------------------------
117  call self%patch_t%init
118  call self%idx%init (5, self%mhd)
119  self%type = 'mhd_t'
120  do i=1,3
121  self%mesh(i)%h(self%idx%px+i-1) = -0.5
122  if (self%mhd) self%mesh(i)%h(self%idx%bx+i-1) = -0.5
123  end do
124  if (io%verbose>1 .and. self%track) then
125  do i=1,3
126  print '("h:",8f5.2)', self%mesh(i)%h
127  end do
128  end if
129  call self%gpatch_t%init
130  self%unsigned(self%idx%d) = .true.
131  self%pervolume(self%idx%s) = .true.
132  call self%extras_t%init
133  !-----------------------------------------------------------------------------
134  ! Self tests, only on one thread
135  !-----------------------------------------------------------------------------
136  !$omp critical (test_cr)
137  if (first_time) then
138  first_time = .false.
139  call stagger_test
140  call self%test
141  end if
142  !$omp end critical (test_cr)
143  call trace%end()
144 END SUBROUTINE init
145 
146 !===============================================================================
147 SUBROUTINE test (self)
148  class(mhd_t):: self
149  integer :: ix, iy, iz, m(3), n(3), n0(3), i1
150  real :: fx, fy, fz, eps
151  logical :: ok, allok
152  class(mhd_t), pointer:: tmp
153  !-----------------------------------------------------------------------------
154  ! Allocate mem to a temporary patch
155  !-----------------------------------------------------------------------------
156  call trace%begin('mhd_t%test')
157  allocate (tmp)
158  tmp%nv = self%nv
159  tmp%mhd = self%mhd
160  call tmp%idx%init (5, tmp%mhd)
161  call tmp%patch_t%init
162  n0 = 2*(tmp%ncell/2)
163  if (any(tmp%n /= n0)) then
164  print *,io%hs
165  print *, 'mhd_t%test not possible with n, no_mans_land =', tmp%n, tmp%no_mans_land
166  print *,io%hs
167  call tmp%dealloc
168  call trace%end()
169  return
170  end if
171  !-----------------------------------------------------------------------------
172  ! Give it initial values (linear compression), also in the guard zones
173  !-----------------------------------------------------------------------------
174  do iz=tmp%mesh(3)%lb,tmp%mesh(3)%ub
175  fz = sin(4.*cgs%pi*tmp%mesh(3)%r(iz)/tmp%size(3))
176  do iy=tmp%mesh(2)%lb,tmp%mesh(2)%ub
177  fy = sin(4.*cgs%pi*tmp%mesh(2)%r(iy)/tmp%size(2))
178  do ix=tmp%mesh(1)%lb,tmp%mesh(1)%ub
179  fx = sin(4.*cgs%pi*tmp%mesh(1)%r(ix)/tmp%size(1))
180  tmp%mem(ix,iy,iz,tmp%idx%d,tmp%it,1) = 1.0+0.1*(fx+fy+fz)
181  tmp%mem(ix,iy,iz,tmp%idx%s,tmp%it,1) = 0.1*(fx+fy+fz)
182  tmp%mem(ix,iy,iz,tmp%idx%px,tmp%it,1) = fx
183  tmp%mem(ix,iy,iz,tmp%idx%py,tmp%it,1) = fy
184  tmp%mem(ix,iy,iz,tmp%idx%pz,tmp%it,1) = fz
185  if (tmp%mhd) then
186  tmp%mem(ix,iy,iz,tmp%idx%bx,tmp%it,1) = fy
187  tmp%mem(ix,iy,iz,tmp%idx%by,tmp%it,1) = fz
188  tmp%mem(ix,iy,iz,tmp%idx%bz,tmp%it,1) = fx
189  end if
190  end do
191  end do
192  end do
193  !-----------------------------------------------------------------------------
194  ! Evaluate the PDE and print the time derivatives (which are placed in the 2nd
195  ! time slot and the 2nd variable slot)
196  !-----------------------------------------------------------------------------
197  tmp%do_force = .false.
198  call tmp%pde
199  tmp%do_force = .true.
200  m = (1+tmp%gn)/2
201  allok = .true.
202  eps = 3e-6/minval(tmp%ds)
203  print *,io%hl
204  print *,'minval(ds), eps =', minval(tmp%ds), eps
205  n = nint(tmp%size*0.5_8/tmp%ds)
206  do ix=tmp%mesh(1)%li,tmp%mesh(1)%ui
207  i1 = mod(ix-tmp%mesh(1)%li+n(1),2*n(1)) + tmp%mesh(1)%li
208  ok = all(abs(tmp%mem(ix,m(2),m(3),1:5,tmp%it,2)- &
209  tmp%mem(i1,m(2),m(3),1:5,tmp%it,2)) < eps)
210  allok = allok .and. ok
211  if (io%verbose>0 .or. io_unit%do_validate .or. .not.ok) &
212  print 1,ix,tmp%mem(ix,m(2),m(3),1:5,tmp%it,2), &
213  tmp%mem(i1,m(2),m(3),1:5,tmp%it,2)
214  1 format("mhd_t%test:",i4,1p,2(2x,5e14.5))
215  end do
216  do iy=tmp%mesh(2)%li,tmp%mesh(2)%ui
217  i1 = mod(iy-tmp%mesh(2)%li+n(2),2*n(2)) + tmp%mesh(2)%li
218  ok = all(abs(tmp%mem(m(1),iy,m(3),1:5,tmp%it,2)- &
219  tmp%mem(m(1),i1,m(3),1:5,tmp%it,2)) < eps)
220  allok = allok .and. ok
221  if (io%verbose>0 .or. io_unit%do_validate .or. .not.ok) &
222  print 1,iy,tmp%mem(m(1),iy,m(3),1:5,tmp%it,2), &
223  tmp%mem(m(1),i1,m(3),1:5,tmp%it,2)
224  end do
225  do iz=tmp%mesh(3)%li,tmp%mesh(3)%ui
226  i1 = mod(iz-tmp%mesh(3)%li+n(3),2*n(3)) + tmp%mesh(3)%li
227  ok = all(abs(tmp%mem(m(1),m(2),iz,1:5,tmp%it,2)- &
228  tmp%mem(m(1),m(2),i1,1:5,tmp%it,2)) < eps)
229  allok = allok .and. ok
230  if (io%verbose>0 .or. io_unit%do_validate .or. .not.ok) &
231  print 1,iz,tmp%mem(m(1),m(2),iz,1:5,tmp%it,2), &
232  tmp%mem(m(1),m(2),i1,1:5,tmp%it,2)
233  end do
234  if (allok) then
235  print *,'mhd_t%test passed'
236  print *,io%hl
237  else
238  print *,'mhd_t%test failed'
239  print *,io%hl
240  call io%abort('mhd_t%test')
241  end if
242  call tmp%dealloc
243  deallocate (tmp)
244  call trace%end()
245 END SUBROUTINE test
246 
247 !===============================================================================
248 SUBROUTINE pde (self)
249  class(mhd_t):: self
250  !.............................................................................
251  real, dimension(:,:,:), pointer:: d, s, dddt, dsdt
252  real, dimension(:,:,:), pointer:: dpdtx, dpdty, dpdtz, dbdtx, dbdty, dbdtz
253  real, dimension(:,:,:), pointer:: bx, by, bz, px, py, pz, phi
254  real, dimension(:,:,:), pointer:: ux, uy, uz, pg
255  real, dimension(:,:,:,:), pointer:: u, p, b, dpdt, dbdt
256  real, dimension(:,:,:), pointer:: ex, ey, ez, jx, jy, jz
257  real, dimension(:,:,:), pointer:: fdx, fdy, fdz, ldx, ldy, ldz, ddx, ddy, ddz
258  real, dimension(:,:,:), pointer:: fxx, fxy
259  real, dimension(:,:,:), pointer:: fyy, fyz
260  real, dimension(:,:,:), pointer:: fzx, fzz
261  real, dimension(:,:,:), allocatable:: du, lnd, ss, cs, pa, pp, pb, q, uu
262  real, dimension(:,:,:), allocatable:: txx, txy
263  real, dimension(:,:,:), allocatable:: tyy, tyz
264  real, dimension(:,:,:), allocatable:: tzx, tzz
265  real, dimension(:,:,:,:), allocatable, target:: emf, j
266  real, dimension(:,:,:,:), allocatable, target:: fx, fy, fz, fd, ld, dd, fm
267  integer :: omp_get_num_threads, omp_get_thread_num
268  logical, save:: first_time=.true.
269  !.............................................................................
270  real(8):: ds(3), dsmax
271  integer, save:: itimer=0
272  real:: u_max
273  !-----------------------------------------------------------------------------
274  call trace%begin('mhd_t%pde',itimer=itimer)
275  associate(i => self%idx)
276  self%d => self%mem(:,:,:,i%d ,self%it,1); self%dddt => self%mem(:,:,:,i%d ,self%it,2)
277  self%s => self%mem(:,:,:,i%s ,self%it,1); self%dsdt => self%mem(:,:,:,i%s ,self%it,2)
278  self%p => self%mem(:,:,:,i%px:i%pz,self%it,1); self%dpdt => self%mem(:,:,:,i%px:i%pz,self%it,2)
279  d=>self%d; dddt=>self%dddt; s=>self%s; dsdt=>self%dsdt;
280  p=>self%p; dpdt=>self%dpdt
281  px=>p(:,:,:,1); py=>p(:,:,:,2); pz=>p(:,:,:,3);
282  dpdtx=>dpdt(:,:,:,1); dpdty=>dpdt(:,:,:,2); dpdtz=>dpdt(:,:,:,3);
283  if (self%mhd) then
284  self%B => self%mem(:,:,:,i%bx:i%bz,self%it,1); self%dBdt => self%mem(:,:,:,i%bx:i%bz,self%it,2)
285  b=>self%B; dbdt=>self%dBdt
286  bx=>b(:,:,:,1); by=>b(:,:,:,2); bz=>b(:,:,:,3);
287  dbdtx=>dbdt(:,:,:,1); dbdty=>dbdt(:,:,:,2); dbdtz=>dbdt(:,:,:,3);
288  end if
289  end associate
290  call allocate_vectors_a (self%gn, ld, dd, emf, j, fd, fx, fy, fz, fm)
291  call allocate_scalars_a (self%gn, lnd, pa, du, ss, cs, pp, pb, q)
292  if (self%do_pebbles) then
293  u => self%vgas(:,:,:,:,self%it)
294  pg => self%pressure(:,:,:,self%it)
295  self%density => self%mem(:,:,:,1,:,1)
296  else
297  call allocate_scalars (self%gn, pg)
298  call allocate_vectors (self%gn, u)
299  end if
300  ux=>u(:,:,:,1); uy=>u(:,:,:,2); uz=>u(:,:,:,3)
301  ex=>emf(:,:,:,1); ey=>emf(:,:,:,2); ez=>emf(:,:,:,3)
302  jx=>j(:,:,:,1); jy=>j(:,:,:,2); jz=>j(:,:,:,3)
303  fxx=>fx(:,:,:,1); fxy=>fx(:,:,:,2)
304  fyy=>fy(:,:,:,2); fyz=>fy(:,:,:,3)
305  fzx=>fz(:,:,:,1); fzz=>fz(:,:,:,3)
306  fdx=>fd(:,:,:,1); fdy=>fd(:,:,:,2); fdz=>fd(:,:,:,3)
307  ldx=>ld(:,:,:,1); ldy=>ld(:,:,:,2); ldz=>ld(:,:,:,3)
308  ddx=>dd(:,:,:,1); ddy=>dd(:,:,:,2); ddz=>dd(:,:,:,3)
309  !-----------------------------------------------------------------------------
310  ! Compute velocity U=momentum/density, leaving down shifted density and log
311  ! density for reuse below. U is valid in (2:gn)
312  !-----------------------------------------------------------------------------
313  lnd = log(d) ! ops: 0 0 0 0 3
314  ld = down(lnd)
315  dd = exp(ld) ! ops: 0 0 0 0 3
316  u = p/dd ! ops: 0 0 3
317  !-----------------------------------------------------------------------------
318  ! Gas pressure.
319  !-----------------------------------------------------------------------------
320  ss = s/d ! ops: 0 0 1
321  call self%gas_pressure_sub (d, s, ss, pg)
322  self%pg => pg
323  !-----------------------------------------------------------------------------
324  ! Conservation of momentum -- first the diagonal part of the stress tensor
325  ! The diagonal fxx is valid in (2:gn-1)
326  !-----------------------------------------------------------------------------
327  call allocate_scalars_a (self%gn, txx, tyy, tzz)
328  ds = self%ds
329  txx = ddxup(ds,ux)
330  tyy = ddyup(ds,uy)
331  tzz = ddzup(ds,uz)
332  !-----------------------------------------------------------------------------
333  ! Divergence of velocity, and artificial pressure. The divergence is valid
334  ! in (2,gn-1)
335  !-----------------------------------------------------------------------------
336  du = txx+tyy+tzz ! convergence rate ! ops: 0 1 0
337  dsmax = maxval(ds,self%n > 2)
338  du = -dsmax*min(du,0.0) ! positive part
339  self%u_max = 0.0
340  ! Diagonal Reynolds stress
341  fxx = xup(ux)**2
342  fyy = yup(uy)**2
343  fzz = zup(uz)**2
344  call allocate_scalars_a (self%gn, uu)
345  uu = sqrt(fxx+fyy+fzz)
346  ! Artificial pressure, possibly with 3-point smoothing
347  if (do_smooth) then
348  pa = self%nu(2)*d*stagger%xyzsm(du)**2
349  else
350  pa = self%nu(2)*d*du**2 ! ops: 0 3 0
351  end if
352  !-----------------------------------------------------------------------------
353  ! MHHD stress diagonal
354  !-----------------------------------------------------------------------------
355  if (self%mhd) then
356  ! Diagonal (negative) Maxwell stress (borrowing emf for scratch)
357  ex = xup(bx)**2 ! mhd: 0 1 0
358  ey = yup(by)**2 ! mhd: 0 1 0
359  ez = zup(bz)**2 ! mhd: 0 1 0
360  ! Magnetic pressure and fast mode speed
361  pb = 0.5*(ex+ey+ez) ! mhd: 2 1 0
362  cs = sqrt((self%gamma*pg+2.0*pb + pa)/d) ! mhd: 1 1 0
363  if (verbose>2) then
364  u_max = self%fmaxval(cs,outer=.true.)
365  self%u_max = max(self%u_max,u_max)
366  print 1, self%id, 'u_max(ca) ', u_max, self%u_max
367  1 format(i6,2x,a,1p,2g12.3)
368  end if
369  ! Isotropic pressure part
370  pp = pg + pa + pb ! mhd: 1 0 0
371  ! Total diagonal stress
372  !du = (dsmax*self%nu(1))*cs + (dsmax*self%nu(7))*uu ! ops: 0 1 0
373  du = (dsmax*self%nu(1))*(cs + 2.*uu) ! ops: 0 1 0
374  !either use Maxwell stress tensor or Lorentz force
375  if (do_maxwell) then
376  fxx = pp + d*(fxx - du*txx) - ex ! mhd: 1 0 0
377  fyy = pp + d*(fyy - du*tyy) - ey ! mhd: 1 0 0
378  fzz = pp + d*(fzz - du*tzz) - ez ! mhd: 1 0 0
379  else
380  fxx = pp + d*(fxx - du*txx) ! mhd: 2 2 0
381  fyy = pp + d*(fyy - du*tyy) ! mhd: 2 2 0
382  fzz = pp + d*(fzz - du*tzz) ! mhd: 2 2 0
383  end if
384  !-----------------------------------------------------------------------------
385  ! HD stress diagonal
386  !-----------------------------------------------------------------------------
387  else
388  ! Sound speed + velocity
389  cs = sqrt((self%gamma*pg+pa)/d) ! ops: 1 1 1 1 0
390  if (verbose>2) then
391  u_max = self%fmaxval(cs,outer=.true.)
392  self%u_max = max(self%u_max,u_max)
393  print 1, self%id, 'u_max(cs) ', u_max, self%u_max
394  end if
395  ! Isotropic pressure part
396  pp = pg + pa ! ops: 1 0 0
397  ! Total diagonal stress, with density factor
398  !---------------------------------------------------------------------------
399  ! Incoming fxx is valid in (2:gn-1), as is du and ddxup(Ux)
400  !---------------------------------------------------------------------------
401  !du = (dsmax*self%nu(1))*cs + (dsmax*self%nu(7))*uu ! ops: 0 1 0
402  du = (dsmax*self%nu(1))*(cs + 2.*uu) ! ops: 0 1 0
403  fxx = pp + d*(fxx - du*txx) ! ops: 2 2 0
404  fyy = pp + d*(fyy - du*tyy) ! ops: 2 2 0
405  fzz = pp + d*(fzz - du*tzz) ! ops: 2 2 0
406  end if
407  if (first_time) &
408  print *, stagger%count, ' stagger calls', stagger%flops/product(real(self%n)),' flops per interior point'
409  u_max = self%fmaxval(cs + uu,outer=.true.)
410  self%u_max = max(self%u_max,u_max)
411  if (verbose>2) print 1, self%id, 'u_max(cs+u) ', u_max, self%u_max
412  !-----------------------------------------------------------------------------
413  ! du has units area per unit time, which gives velocity squared times mass
414  ! density per unit time = energy per unit volume and time
415  !-----------------------------------------------------------------------------
416  q = d*du*(txx**2 + tyy**2 + tzz**2) ! ops: 2 5 0
417  call deallocate_scalars_a (txx, tyy, tzz, uu)
418  if (verbose==1) print *, 'average(Q):', self%faver(q)
419  !-----------------------------------------------------------------------------
420  ! Conservation of mass
421  !-----------------------------------------------------------------------------
422  fd = -(self%nu(4)*dd*down(du))*ddown(ds,lnd) ! diffusive mass flux ! ops: 0 9 0
423  fm = p+fd ! mass flux ! ops: 3 0 0
424  dddt = -div(ds,fm) ! density derivative
425  u_max = self%cdtd*dsmax*self%fmaxval(abs(dddt/d),outer=.true.) ! ops: 0 0 1
426  self%u_max = max(self%u_max,u_max)
427  !if (verbose>2) print 1, self%id, 'u_max(dddt) ', u_max, self%u_max
428  !-----------------------------------------------------------------------------
429  ! Off-diagonal part of Reynolds stress, for now without the density factor.
430  ! Note that, because of the symmmetry, we define only three components
431  !-----------------------------------------------------------------------------
432  fxy = ydn(ux)*xdn(uy) ! ops: 0 1 0
433  fyz = zdn(uy)*ydn(uz) ! ops: 0 1 0
434  fzx = xdn(uz)*zdn(ux) ! ops: 0 1 0
435  !-----------------------------------------------------------------------------
436  ! Add viscous part -- edge-centered. Tij should really be divided by two, in
437  ! which case twice the product nu*Tij**2 should be added, so here a factor 0.5
438  ! is needed in the contribution to Q
439  !-----------------------------------------------------------------------------
440  call allocate_scalars_a (self%gn, txy, tyz, tzx)
441  txy = ddxdn(ds,uy)+ddydn(ds,ux) ! ops: 0 1 0
442  tyz = ddydn(ds,uz)+ddzdn(ds,uy) ! ops: 0 1 0
443  tzx = ddzdn(ds,ux)+ddxdn(ds,uz) ! ops: 0 1 0
444  fxy = fxy - self%nu(3)*0.5*xdn1(ydn1(du))*txy ! ops: 1 3 0
445  fyz = fyz - self%nu(3)*0.5*ydn1(zdn1(du))*tyz ! ops: 1 3 0
446  fzx = fzx - self%nu(3)*0.5*zdn1(xdn1(du))*tzx ! ops: 1 3 0
447  q = q + self%nu(3)*.25*d*du*(xup(yup(txy))**2 + & ! ops: 2 5 0
448  yup(zup(tyz))**2 + & ! ops: 1 1 0
449  zup(xup(tzx))**2) ! ops: 1 0 0
450  call deallocate_scalars_a (txy, tyz, tzx)
451  !-----------------------------------------------------------------------------
452  ! Multiply with edge-centered density
453  !-----------------------------------------------------------------------------
454  fxy = exp(ydn(ldx))*fxy ! ops: 0 1 0 0 1
455  fyz = exp(zdn(ldy))*fyz ! ops: 0 1 0 0 1
456  fzx = exp(xdn(ldz))*fzx ! ops: 0 1 0 0 1
457  !-----------------------------------------------------------------------------
458  ! Off-diagonal Maxwell stress
459  !-----------------------------------------------------------------------------
460  if (self%mhd.and.do_maxwell) then
461  fxy = fxy - ydn(bx)*xdn(by) ! mhd: 1 1 0
462  fyz = fyz - zdn(by)*ydn(bz) ! mhd: 1 1 0
463  fzx = fzx - xdn(bz)*zdn(bx) ! mhd: 1 1 0
464  end if
465  !-----------------------------------------------------------------------------
466  ! Consistent addition of mass flux related momentum flux, symmetrized. The
467  ! diagonal part is analogous, a centered version of \nu \rho U_i \grad\ln\rho,
468  ! while the off-diagonal part needs to be properly symmetrized.
469  !-----------------------------------------------------------------------------
470  if (self%nu(4) > 0.0) then
471  pa = d*du*self%nu(4) ! ops: 0 2 0
472  !---------------------------------------------------------------------------
473  ! fxx, pa, ddxup(ldx), and xup(Ux) are all valid in (2:gn-1)
474  !---------------------------------------------------------------------------
475  fxx = fxx - pa*ddxup(ds,ldx)*xup(ux) ! ops: 1 2 0
476  fyy = fyy - pa*ddyup(ds,ldy)*yup(uy) ! ops: 1 2 0
477  fzz = fzz - pa*ddzup(ds,ldz)*zup(uz) ! ops: 1 2 0
478  fxy = fxy + (xdn1(fdy)*ydn(ux)+ydn1(fdx)*xdn(uy)) ! ops: 2 2 0
479  fyz = fyz + (ydn1(fdz)*zdn(uy)+zdn1(fdy)*ydn(uz)) ! ops: 2 2 0
480  fzx = fzx + (zdn1(fdx)*xdn(uz)+xdn1(fdz)*zdn(ux)) ! ops: 2 2 0
481  end if
482  !-----------------------------------------------------------------------------
483  ! Momentum equations. The ddxdn(fxx) part is valid in (3:gn-1), and hence px
484  ! only needs guard zones in (1:2) and (gn:gn)
485  !-----------------------------------------------------------------------------
486  if (do_2nd_div) then
487  dpdtx = - ddxdn1(ds,fxx) - ddyup1(ds,fxy) - ddzup1(ds,fzx) ! ops: 2 0 0
488  dpdty = - ddxup1(ds,fxy) - ddydn1(ds,fyy) - ddzup1(ds,fyz) ! ops: 2 0 0
489  dpdtz = - ddxup1(ds,fzx) - ddyup1(ds,fyz) - ddzdn1(ds,fzz) ! ops: 2 0 0
490  else
491  dpdtx = - ddxdn(ds,fxx) - ddyup(ds,fxy) - ddzup(ds,fzx) ! ops: 2 0 0
492  dpdty = - ddxup(ds,fxy) - ddydn(ds,fyy) - ddzup(ds,fyz) ! ops: 2 0 0
493  dpdtz = - ddxup(ds,fzx) - ddyup(ds,fyz) - ddzdn(ds,fzz) ! ops: 2 0 0
494  end if
495  !-----------------------------------------------------------------------------
496  ! External force; e.g. point source gravity. Note that the density factor is
497  ! applied here, and should not be included in the function definition
498  !-----------------------------------------------------------------------------
499  if (self%do_force) then
500  if (allocated(self%force_per_unit_mass)) then
501  dpdt = dpdt + dd*self%force_per_unit_mass ! ops: 2 2 0
502  end if
503  if (allocated(self%force_per_unit_volume)) then
504  dpdt = dpdt + self%force_per_unit_volume
505  end if
506  end if
507  !-----------------------------------------------------------------------------
508  ! Conservation of entropy. The formulation that works best so far is to make
509  ! the diffusive entropy flux directly proportional to the diffusive mass flux,
510  ! plus a contribution proportional to the gradient of specific entropy ss.
511  !-----------------------------------------------------------------------------
512  if (self%gamma==1d0) then
513  dsdt = 0.0
514  else
515  fm = fm*down(ss) - (self%nu(5)*dd*down(du))*ddown(ds,ss) ! ops: 1 4 0
516  !---------------------------------------------------------------------------
517  ! Q is heating per unit volume (energy per unit volume and time), which when
518  ! divided by presssure (energy per unit volume) becomes entropy (unitless)
519  ! per unit mass. Multiplied by mass density it becomes mass-entropy per
520  ! unit volume, which is our units for s.
521  !---------------------------------------------------------------------------
522  dsdt = -div(ds,fm) + q*d/pg ! ops: 1 1 1
523  end if
524  !-----------------------------------------------------------------------------
525  ! Induction equation
526  !-----------------------------------------------------------------------------
527  if (self%mhd) then
528  !---------------------------------------------------------------------------
529  ! Electric current, resistive electric field and magnetic dissipation
530  !---------------------------------------------------------------------------
531  j = curl_dn(ds,b)
532  emf = self%nu(6)*edge(du)*j ! mhd: 0 6 0
533  q = yup(zup(jx*ex)) + zup(xup(jy*ey)) + xup(yup(jz*ez)) ! mhd: 2 3 0
534  !---------------------------------------------------------------------------
535  ! Lorentz force, used instead of Maxwell stress.
536  !---------------------------------------------------------------------------
537  if (.not.do_maxwell) then
538  !borrowing fxyz for scratch to store cell-centered B values
539  ! Bx -> fxx, By -> fyy, Bz -> fzz
540  fxx = xup(bx)
541  fyy = yup(by)
542  fzz = zup(bz)
543  dpdtx = dpdtx + zup(jy)*xdn(fzz) - yup(jz)*xdn(fyy)
544  dpdty = dpdty + xup(jz)*ydn(fxx) - zup(jx)*ydn(fzz)
545  dpdtz = dpdtz + yup(jx)*zdn(fyy) - xup(jy)*zdn(fxx)
546  end if
547  !---------------------------------------------------------------------------
548  ! Advective electric field on cell edges
549  !---------------------------------------------------------------------------
550  ex = ex + ydn(uz)*zdn(by) - zdn(uy)*ydn(bz) ! mhd: 2 2 0
551  ey = ey + zdn(ux)*xdn(bz) - xdn(uz)*zdn(bx) ! mhd: 2 2 0
552  ez = ez + xdn(uy)*ydn(bx) - ydn(ux)*xdn(by) ! mhd: 2 2 0
553  call non_ideal%update (self%gn, minval(ds), d, jx, jy, jz, bx, by, bz, q, emf, self%u_max)
554  if (self%gamma==1d0) then
555  dsdt = 0.0
556  else
557  dsdt = dsdt + q*d/pg ! mhd: 1 1 1
558  end if
559  dbdt = -curl_up(ds,emf)
560  end if
561  !-----------------------------------------------------------------------------
562  ! Time step limitation for rate of change dS/dt = (ds/dt-(s/d)*dddt)/d
563  !-----------------------------------------------------------------------------
564  u_max = self%cdtd*dsmax*self%fmaxval(abs((dsdt-s*dddt/d)/d),outer=.true.) ! ops: 1 1 2
565  self%u_max = max(self%u_max,u_max)
566  if (verbose>2) print 1, self%id, 'u_max(dsdt) ', u_max, self%u_max
567  !-----------------------------------------------------------------------------
568  call deallocate_vectors_a (ld, dd, emf, j, fd, fx, fy, fz, fm)
569  call deallocate_scalars_a (lnd, pa, du, ss, cs, pp, pb, q)
570  if (.not.self%do_pebbles) then
571  call deallocate_vectors (u)
572  call deallocate_scalars (pg)
573  end if
574  call trace%end (itimer)
575  if (mpi%master .and. (first_time .or. flop_count)) then
576  print *, stagger%count, ' stagger calls', stagger%flops/product(real(self%n)),' flops per interior point'
577  stagger%count=0
578  stagger%flops=0
579  first_time = .false.
580  end if
581 END SUBROUTINE pde
582 
583 !===============================================================================
584 !> Take a complete update step. This can be overloaded in higher levels, e.g.
585 !> in experiment_mod, where boundary conditions may be inserted
586 !===============================================================================
587 SUBROUTINE update (self)
588  class(mhd_t):: self
589  integer, save:: itimer=0
590  !-----------------------------------------------------------------------------
591  call trace%begin('mhd_t%update',itimer=itimer)
592  !-----------------------------------------------------------------------------
593  ! Output the state after guard zone downloads, selfgravity, and boundary
594  ! conditions. If another call is placed at an earlier point in the task handling
595  ! the call here will have no effect, since out_next has already been updated.
596  !-----------------------------------------------------------------------------
597  call self%output
598  call self%pde
599  call self%courant_condition
600  if (io%verbose>2) &
601  call self%check_density('before timestep')
602  call timestep%update (self%id, self%iit, self%dt, self%mem, self%mesh, self%lock)
603  call self%check_density(' after timestep')
604  !-----------------------------------------------------------------------------
605  call self%counter_update
606  call trace%end (itimer)
607 END SUBROUTINE update
608 
609 !===============================================================================
610 !> Add extra, MHD-related information to the patch .txt output file
611 !===============================================================================
612 SUBROUTINE output (self)
613  class(mhd_t):: self
614  character(len=64):: filename
615  logical:: exists
616  real(8):: out_next
617  !-----------------------------------------------------------------------------
618  out_next = self%out_next
619  call trace%begin('mhd_t%output')
620  call self%extras_t%output
621  if (self%time >= self%out_next .and. io%do_output .and. io%do_legacy) then
622  !$omp critical (out_cr)
623  write (filename,'(a,i4.4,"_",i4.4,".txt")') trim(io%outputname), &
624  self%id, self%iout
625  inquire(file=filename, exist=exists)
626  if (exists) then
627  open (io%data_unit, file=trim(filename), form='formatted', status='old', &
628  position='append')
629  write (io%data_unit,'(8f10.4)') self%nu
630  close (io%data_unit)
631  end if
632  !$omp end critical (out_cr)
633  end if
634  call trace%end()
635 END SUBROUTINE output
636 
637 !===============================================================================
638 FUNCTION gas_pressure (self) RESULT (pg)
639  class(mhd_t):: self
640  real, dimension(:,:,:), pointer:: d, s
641  real, dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
642  integer, save:: nprint=3, n(3)
643  !-----------------------------------------------------------------------------
644  if (trim(self%eos)=='ideal') then
645  d => self%mem(:,:,:,self%idx%d,self%it,1)
646  s => self%mem(:,:,:,self%idx%s,self%it,1)
647  if (self%gamma==1d0) then
648  s = 0.0
649  pg = d*self%csound**2
650  else
651  s => self%mem(:,:,:,self%idx%s,self%it,1)
652  pg = d**self%gamma*exp(s/d*(self%gamma-1))
653  !pg = exp(self%gamma*log(d))*exp(ss*(self%gamma-1.0))
654  end if
655  else
656  n=shape(d)
657 #if defined __PGI
658  call eos%lookup (n, x=s, y=d)
659 #else
660  call eos%lookup (n, x=s, y=d, pg=pg)
661 #endif
662  end if
663  if (io%master .and. nprint>0) then
664  nprint = nprint - 1
665  print *, 'mhd_t%gas_pressure: min, max =', self%fminval(pg), self%fmaxval(pg)
666  end if
667 END FUNCTION gas_pressure
668 
669 !===============================================================================
670 SUBROUTINE gas_pressure_sub (self, d, s, ss, pg)
671  class(mhd_t):: self
672  real, dimension(:,:,:), pointer:: d, s, pg
673  real, dimension(:,:,:):: ss
674  integer:: i1, i2, i3
675  integer, save:: nprint=3
676  !-----------------------------------------------------------------------------
677  if (trim(self%eos)=='ideal') then
678  if (self%gamma==1d0) then
679  ss = 0.0
680  pg = d*self%csound**2
681  else
682  pg = exp(self%gamma*log(d))*exp(ss*(self%gamma-1.0))
683  end if
684  else
685  print *,'yyy'
686  call eos%lookup (shape(d), x=s, y=d, pg=pg)
687  end if
688  if (io%master .and. nprint>0) then
689  nprint = nprint - 1
690  print *, 'mhd_t%gas_pressure: min, max =', self%fminval(pg), self%fmaxval(pg)
691  end if
692 END SUBROUTINE gas_pressure_sub
693 
694 END MODULE mhd_mod
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...
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
non-ideal MHD module
RAMSES Godunov solvers, use of guard zones; specifically in HLLD.
Definition: mhd_mod.f90:16
Equation of state module for any sort of tables, provided by a reader.
Definition: eos_mod.f90:4
Fundamental constants in CGS and SI units.
Definition: units_mod.f90:4
Simple initical condition module for tests.
Definition: initial_mod.f90:7
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Definition: index_mod.f90:7
Definition: io_mod.f90:4
6th order stagger operators, with self-test procedure
Definition: stagger_2nd.f90:4