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 aux_mod
14  USE bits_mod
15  USE omp_mod
16  USE math_mod
17  USE omp_timer_mod
18  USE mpi_mod
19  USE trace_mod
20  USE kinds_mod
21  USE units_mod
22  USE scaling_mod
23  USE scalar_mod
24  USE vector_mod
25  USE vector_ops
26  USE stagger_mod
27  USE extras_mod
28  USE index_mod
29  USE timestep_mod
30  USE eos_mod
31  USE non_ideal_mod
32  USE initial_mod
33  implicit none
34  private
35  type, extends(extras_t), public:: mhd_t
36  real, dimension(:,:,:), pointer :: d, dddt, e, dedt
37  real, dimension(:,:,:,:), pointer :: p, dpdt, b, dbdt
38  real, dimension(:,:,:), pointer :: lnd, ee, pg=>null(), tt=>null()
39  real, dimension(:,:,:), allocatable:: qrad
40  procedure(divb_clean1), pointer :: divb_clean => null()
41  real:: nu(6)
42  real:: cdtd
43  real:: maxj=0.0
44  logical:: mhd=.false.
45  logical:: first_time=.true.
46  contains
47  procedure:: pre_init
48  procedure:: init
49  procedure:: pde
50  procedure:: pre_update
51  procedure:: update
52  procedure:: output
53  procedure:: gas_pressure
54  procedure:: gas_pressure_sub
55  end type
56  logical, save:: unsigned=.true.
57  logical, save:: do_test=.false.
58  logical, save:: do_smooth=.false.
59  logical, save:: do_maxwell = .true.
60  logical, save:: do_temperature = .false.
61  integer, save:: verbose=0
62  integer, save:: divb_cleaner=2
63 CONTAINS
64 
65 !===============================================================================
66 !> Read the namelist parameters, store them, determine the type of the solver
67 !> (HD or MHD) and set the indices.
68 !===============================================================================
69 SUBROUTINE pre_init (self)
70  class(mhd_t), target:: self
71  type(index_t):: idx
72  integer:: i, iv
73  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
74  real(8), save:: gamma=1.4_8
75  logical, save:: first_time=.true.
76  logical, save:: mhd=.false.
77  character(len=16), save:: eos
78  namelist /stagger_params/ nu, courant, gamma, csound, cdtd, hardwire, eos, &
79  do_smooth, unsigned, do_maxwell, do_temperature, verbose, divb_cleaner, &
80  do_test
81  !----------------------------------------------------------------------------
82  call trace%begin('mhd_t%pre_init')
83  !$omp critical (input_cr)
84  if (first_time) then
85  first_time = .false.
86  eos = self%eos
87  mhd = self%mhd
88  rewind(io%input)
89  read (io%input, stagger_params)
90  if (io%master) write (*, stagger_params)
91  call non_ideal%init
92  call timestep%init
93  ! fool-proofing
94  if (timestep%time_order <= 0) call mpi%abort("The Stagger solvers require `time_order` > 0. Abort!")
95  end if
96  !$omp end critical (input_cr)
97  self%eos = eos
98  self%mhd = mhd
99  self%nu = nu
100  self%gamma = gamma
101  self%csound = csound
102  self%courant = courant
103  self%cdtd = cdtd
104  !-----------------------------------------------------------------------------
105  ! Read IC parameters, which determine if this is HD or MHD. The IC init is
106  ! called again, from extras_t%init, when gamma is known and can be passed on.
107  !-----------------------------------------------------------------------------
108  self%initial%mhd = self%mhd
109  call self%initial%init (self%kind, real(self%gamma))
110  self%mhd = self%initial%mhd
111  if (self%mhd) then
112  self%kind = 'stagger2e_mhd_patch'
113  else
114  self%kind = 'stagger2e_hd_patch'
115  end if
116  if (self%nv == 0) then
117  if (self%mhd) then
118  self%nv = 8
119  else
120  self%nv = 5
121  end if
122  else
123  self%nv = self%nv + 5
124  if (self%mhd) self%nv = self%nv + 3
125  end if
126  select case(divb_cleaner)
127  case(1)
128  self%divb_clean => divb_clean1
129  case(2)
130  self%divb_clean => divb_clean2
131  case default
132  end select
133  call self%idx%init (5, self%mhd)
134  ! ----------------------------------------------------------------------------
135  ! Call extras pre-init
136  ! ----------------------------------------------------------------------------
137  call self%extras_t%pre_init
138  call trace%end()
139 END SUBROUTINE pre_init
140 
141 !===============================================================================
142 !> The construction with a double test on first_time avoids entering a critical
143 !> region once the variable has been set.
144 !===============================================================================
145 SUBROUTINE init (self)
146  class(mhd_t):: self
147  integer :: i
148  ! ----------------------------------------------------------------------------
149  call trace%begin('mhd_t%init')
150  !-----------------------------------------------------------------------------
151  ! Allocate memory and mesh, etc
152  !-----------------------------------------------------------------------------
153  self%type = 'mhd_t'
154  call self%patch_t%init
155  do i=1,3
156  self%mesh(i)%h(self%idx%px+i-1) = -0.5
157  if (self%mhd) self%mesh(i)%h(self%idx%bx+i-1) = -0.5
158  end do
159  if (verbose>1 .and. self%track) then
160  do i=1,3
161  print '("h:",8f5.2)', self%mesh(i)%h
162  end do
163  end if
164  call self%patch_t%init
165  select case(divb_cleaner)
166  case(1)
167  self%divb_clean => divb_clean1
168  case(2)
169  self%divb_clean => divb_clean2
170  case default
171  call io%abort('mhd_mod:: divb cleaner method not understood')
172  end select
173  call self%gpatch_t%init
174  self%unsigned(self%idx%d) = unsigned
175  self%pervolume(self%idx%e) = unsigned
176  ! ----------------------------------------------------------------------------
177  ! init extras
178  ! ----------------------------------------------------------------------------
179  call self%extras_t%init
180  call trace%end()
181 END SUBROUTINE init
182 
183 !===============================================================================
184 SUBROUTINE pde (self)
185  USE scalar_mod
186  class(mhd_t):: self
187  !.............................................................................
188  real, dimension(:,:,:), pointer:: d, e, dddt, dedt
189  real, dimension(:,:,:), pointer:: dpdtx, dpdty, dpdtz, dbdtx, dbdty, dbdtz
190  real, dimension(:,:,:), pointer:: bx, by, bz, px, py, pz, phi, qr
191  real, dimension(:,:,:), pointer:: ux, uy, uz, pg
192  real, dimension(:,:,:,:), pointer:: u, p, b, dpdt, dbdt
193  real, dimension(:,:,:), pointer:: ex, ey, ez, jx, jy, jz
194  real, dimension(:,:,:), pointer:: fdx, fdy, fdz, ldx, ldy, ldz, ddx, ddy, ddz
195  real, dimension(:,:,:), pointer:: fxx, fxy
196  real, dimension(:,:,:), pointer:: fyy, fyz
197  real, dimension(:,:,:), pointer:: fzx, fzz
198  real, dimension(:,:,:), allocatable:: du, ss, cs, pa, pp, pb, q, uu, divu
199  real, dimension(:,:,:), allocatable, target:: ee, lnd
200  real, dimension(:,:,:), allocatable:: txx, txy
201  real, dimension(:,:,:), allocatable:: tyy, tyz
202  real, dimension(:,:,:), allocatable:: tzx, tzz
203  real, dimension(:,:,:,:), allocatable, target:: emf, j
204  real, dimension(:,:,:,:), allocatable, target:: fx, fy, fz, fd, ld, dd, fm
205  integer :: omp_get_num_threads, omp_get_thread_num
206  logical, save:: first_time=.true.
207  !.............................................................................
208  real(8):: ds(3), dsmax
209  real:: u_max
210  integer, save:: itimer=0
211  !-----------------------------------------------------------------------------
212  call trace%begin('mhd_t%pde',itimer=itimer)
213  associate(i => self%idx)
214  self%d => self%mem(:,:,:,i%d ,self%it,1); self%dddt => self%mem(:,:,:,i%d ,self%it,2)
215  self%e => self%mem(:,:,:,i%e ,self%it,1); self%dedt => self%mem(:,:,:,i%e ,self%it,2)
216  self%p => self%mem(:,:,:,i%px:i%pz,self%it,1); self%dpdt => self%mem(:,:,:,i%px:i%pz,self%it,2)
217  d=>self%d; dddt=>self%dddt; e=>self%e; dedt=>self%dedt;
218  p=>self%p; dpdt=>self%dpdt
219  px=>p(:,:,:,1); py=>p(:,:,:,2); pz=>p(:,:,:,3);
220  dpdtx=>dpdt(:,:,:,1); dpdty=>dpdt(:,:,:,2); dpdtz=>dpdt(:,:,:,3);
221  if (self%mhd) then
222  self%B => self%mem(:,:,:,i%bx:i%bz,self%it,1); self%dBdt => self%mem(:,:,:,i%bx:i%bz,self%it,2)
223  b=>self%B; dbdt=>self%dBdt
224  bx=>b(:,:,:,1); by=>b(:,:,:,2); bz=>b(:,:,:,3);
225  dbdtx=>dbdt(:,:,:,1); dbdty=>dbdt(:,:,:,2); dbdtz=>dbdt(:,:,:,3);
226  end if
227  end associate
228  call allocate_vectors_a (self%gn, ld, dd, emf, j, fd, fx, fy, fz, fm)
229  call allocate_scalars_a (self%gn, lnd, pa, du, ee, cs, pp, pb, q, uu, divu)
230  if (self%do_pebbles) then
231  u => self%vgas(:,:,:,:,self%it)
232  pg => self%pressure(:,:,:,self%it)
233  self%density => self%mem(:,:,:,1,:,1)
234  else
235  call allocate_scalars (self%gn, pg)
236  call allocate_vectors (self%gn, u)
237  end if
238  ux=>u(:,:,:,1); uy=>u(:,:,:,2); uz=>u(:,:,:,3)
239  ex=>emf(:,:,:,1); ey=>emf(:,:,:,2); ez=>emf(:,:,:,3)
240  jx=>j(:,:,:,1); jy=>j(:,:,:,2); jz=>j(:,:,:,3)
241  fxx=>fx(:,:,:,1); fxy=>fx(:,:,:,2)
242  fyy=>fy(:,:,:,2); fyz=>fy(:,:,:,3)
243  fzx=>fz(:,:,:,1); fzz=>fz(:,:,:,3)
244  fdx=>fd(:,:,:,1); fdy=>fd(:,:,:,2); fdz=>fd(:,:,:,3)
245  ldx=>ld(:,:,:,1); ldy=>ld(:,:,:,2); ldz=>ld(:,:,:,3)
246  ddx=>dd(:,:,:,1); ddy=>dd(:,:,:,2); ddz=>dd(:,:,:,3)
247  !-----------------------------------------------------------------------------
248  ! Compute velocity U=momentum/density, leaving down shifted density and log
249  ! density for reuse below. U is valid in (2:gn)
250  !-----------------------------------------------------------------------------
251  lnd = log(d) ! ops: 0 0 0 0 3
252  ld = down(lnd)
253  dd = exp(ld) ! ops: 0 0 0 0 3
254  u = p/dd ! ops: 0 0 3
255  !-----------------------------------------------------------------------------
256  ! Gas pressure. The 2nd table variable is entropy per unit mass, effectively
257  ! a log variable, from the table point of view.
258  !-----------------------------------------------------------------------------
259  ee = e/d ! ops: 0 0 1
260  self%lnd => lnd
261  self%ee => ee
262  call self%gas_pressure_sub (d, self%lnd, self%ee, pg)
263  self%pg => pg
264  !-----------------------------------------------------------------------------
265  ! Conservation of momentum -- first the diagonal part of the stress tensor
266  ! The diagonal fxx is valid in (2:gn-1)
267  !-----------------------------------------------------------------------------
268  call allocate_scalars_a (self%gn, txx, tyy, tzz)
269  ds = self%ds
270  txx = ddxup(ds,ux)
271  tyy = ddyup(ds,uy)
272  tzz = ddzup(ds,uz)
273  divu = txx+tyy+tzz ! convergence rate ! ops: 0 2 0
274  !-----------------------------------------------------------------------------
275  ! Divergence of velocity, and artificial pressure. The divergence is valid
276  ! in (2,gn-1)
277  !-----------------------------------------------------------------------------
278  dsmax = maxval(ds,self%n > 2)
279  du = -dsmax*min(divu,0.0) ! positive part
280  if (do_smooth) then
281  pa = self%nu(2)*d*xyzsm(du**2)
282  else
283  pa = self%nu(2)*d*du**2 ! ops: 0 3 0
284  end if
285  !-----------------------------------------------------------------------------
286  ! Diagonal Reynolds stress
287  !-----------------------------------------------------------------------------
288  fxx = xup(ux)**2
289  fyy = yup(uy)**2
290  fzz = zup(uz)**2
291  uu = sqrt(fxx+fyy+fzz)
292  self%u_max = 0.0
293  if (self%mhd) then
294  ! Diagonal Reynolds stress
295  ! Diagonal (negative) Maxwell stress (borrowing emf for scratch)
296  ex = xup(bx)**2 ! mhd: 0 1 0
297  ey = yup(by)**2 ! mhd: 0 1 0
298  ez = zup(bz)**2 ! mhd: 0 1 0
299  ! Magnetic pressure
300  pb = 0.5*(ex+ey+ez) ! mhd: 2 1 0
301  cs = sqrt((self%gamma*pg+2.0*pb + pa)/d) ! mhd: 1 1 0
302  if (verbose>3) then
303  u_max = self%fmaxval(cs,outer=.true.)
304  self%u_max = max(self%u_max,u_max)
305  print 1, 'u_max(ca) ', u_max, self%u_max
306  1 format(1x,a,1p,2g12.3)
307  end if
308  ! Isotropic pressure part
309  pp = pg + pa + pb ! mhd: 1 0 0
310  ! Total diagonal stress
311  du = (dsmax*self%nu(1))*(cs + 2.*uu) ! mhd: 0 2 0
312  !either use Maxwell stress tensor or Lorentz force
313  if (do_maxwell) then
314  fxx = pp + d*(fxx - du*txx) - ex ! mhd: 1 0 0
315  fyy = pp + d*(fyy - du*tyy) - ey ! mhd: 1 0 0
316  fzz = pp + d*(fzz - du*tzz) - ez ! mhd: 1 0 0
317  else
318  fxx = pp + d*(fxx - du*txx) ! mhd: 2 2 0
319  fyy = pp + d*(fyy - du*tyy) ! mhd: 2 2 0
320  fzz = pp + d*(fzz - du*tzz) ! mhd: 2 2 0
321  end if
322  else
323  ! Sound speed + velocity
324  cs = sqrt((self%gamma*pg+pa)/d) ! ops: 1 1 1 1 0
325  if (verbose>3) then
326  u_max = self%fmaxval(cs,outer=.true.)
327  self%u_max = max(self%u_max,u_max)
328  print 1, 'u_max(cs) ', u_max, self%u_max
329  end if
330  pp = pg + pa ! ops: 1 0 0
331  du = (dsmax*self%nu(1))*(cs + 2.*uu) ! ops: 0 2 0
332  !---------------------------------------------------------------------------
333  ! Total diagonal stress, with density factor
334  ! Incoming fxx is valid in (2:gn-1), as is du and ddxup(Ux)
335  !---------------------------------------------------------------------------
336  fxx = pp + d*(fxx - du*txx) ! ops: 2 2 0
337  fyy = pp + d*(fyy - du*tyy) ! ops: 2 2 0
338  fzz = pp + d*(fzz - du*tzz) ! ops: 2 2 0
339  end if
340  u_max = self%fmaxval(cs+uu,outer=.true.)
341  self%u_max = max(self%u_max,u_max)
342  if (verbose>2) print 1, 'u_max(cs+u)', u_max, self%u_max
343  !-----------------------------------------------------------------------------
344  ! du has units area per unit time, which gives velocity squared times mass
345  ! density per unit time = energy per unit volume and time
346  !-----------------------------------------------------------------------------
347  q = d*du*(txx**2 + tyy**2 + tzz**2) ! ops: 2 5 0
348  !if (verbose==1) print *, 'average(Q):', self%faver(Q), &
349  !self%faver(e-d**self%gamma/(self%gamma-1.0)), self%faver(fm(:,:,:,1))
350  !-----------------------------------------------------------------------------
351  ! Conservation of mass
352  !-----------------------------------------------------------------------------
353  fd = -(self%nu(4)*dd*down(du))*ddown(ds,lnd) ! diffusive mass flux ! ops: 0 9 0
354  fm = p+fd ! mass flux ! ops: 3 0 0
355  dddt = -div(ds,fm) ! density derivative
356  u_max = self%cdtd*dsmax*self%fmaxval(abs(dddt/d),outer=.true.) ! ops: 0 0 1
357  self%u_max = max(self%u_max,u_max)
358  if (verbose>2) print 1, 'u_max(dddt)', u_max, self%u_max
359  !-----------------------------------------------------------------------------
360  ! Off-diagonal part of Reynolds stress, for now without the density factor.
361  ! Note that, because of the symmmetry, we define only three components
362  !-----------------------------------------------------------------------------
363  fxy = ydn(ux)*xdn(uy) ! ops: 0 1 0
364  fyz = zdn(uy)*ydn(uz) ! ops: 0 1 0
365  fzx = xdn(uz)*zdn(ux) ! ops: 0 1 0
366  !-----------------------------------------------------------------------------
367  ! Add viscous part -- edge-centered
368  !-----------------------------------------------------------------------------
369  call allocate_scalars_a (self%gn, txy, tyz, tzx)
370  txy = ddxdn(ds,uy)+ddydn(ds,ux) ! ops: 0 1 0
371  tyz = ddydn(ds,uz)+ddzdn(ds,uy) ! ops: 0 1 0
372  tzx = ddzdn(ds,ux)+ddxdn(ds,uz) ! ops: 0 1 0
373  fxy = fxy - self%nu(3)*0.5*xdn1(ydn1(du))*txy ! ops: 1 3 0
374  fyz = fyz - self%nu(3)*0.5*ydn1(zdn1(du))*tyz ! ops: 1 3 0
375  fzx = fzx - self%nu(3)*0.5*zdn1(xdn1(du))*tzx ! ops: 1 3 0
376  q = q + self%nu(3)*0.5*d*du*(xup(yup(txy))**2 + & ! ops: 2 5 0
377  yup(zup(tyz))**2 + & ! ops: 1 1 0
378  zup(xup(tzx))**2) ! ops: 1 0 0
379  call deallocate_scalars_a (txy, tyz, tzx)
380  !-----------------------------------------------------------------------------
381  ! Multiply with edge-centered density
382  !-----------------------------------------------------------------------------
383  fxy = exp(ydn(ldx))*fxy ! ops: 0 1 0 0 1
384  fyz = exp(zdn(ldy))*fyz ! ops: 0 1 0 0 1
385  fzx = exp(xdn(ldz))*fzx ! ops: 0 1 0 0 1
386  !-----------------------------------------------------------------------------
387  ! Off-diagonal Maxwell stress
388  !-----------------------------------------------------------------------------
389  if (self%mhd.and.do_maxwell) then
390  fxy = fxy - ydn(bx)*xdn(by) ! mhd: 1 1 0
391  fyz = fyz - zdn(by)*ydn(bz) ! mhd: 1 1 0
392  fzx = fzx - xdn(bz)*zdn(bx) ! mhd: 1 1 0
393  end if
394  !-----------------------------------------------------------------------------
395  ! Consistent addition of mass flux related momentum flux, symmetrized. The
396  ! diagonal part is analogous, a centered version of \nu \rho U_i \grad\ln\rho,
397  ! while the off-diagonal part needs to be properly symmetrized.
398  !-----------------------------------------------------------------------------
399  if (self%nu(4) > 0.0) then
400  pa = d*du*self%nu(4) ! ops: 0 2 0
401  !---------------------------------------------------------------------------
402  ! fxx, pa, ddxup(ldx), and xup(Ux) are all valid in (2:gn-1)
403  !---------------------------------------------------------------------------
404  fxx = fxx - pa*ddxup(ds,ldx)*xup(ux) ! ops: 1 2 0
405  fyy = fyy - pa*ddyup(ds,ldy)*yup(uy) ! ops: 1 2 0
406  fzz = fzz - pa*ddzup(ds,ldz)*zup(uz) ! ops: 1 2 0
407  fxy = fxy + (xdn1(fdy)*ydn(ux)+ydn1(fdx)*xdn(uy)) ! ops: 2 2 0
408  fyz = fyz + (ydn1(fdz)*zdn(uy)+zdn1(fdy)*ydn(uz)) ! ops: 2 2 0
409  fzx = fzx + (zdn1(fdx)*xdn(uz)+xdn1(fdz)*zdn(ux)) ! ops: 2 2 0
410  end if
411  !-----------------------------------------------------------------------------
412  ! Momentum equations. The ddxdn(fxx) part is valid in (3:gn-1), and hence px
413  ! only needs guard zones in (1:2) and (gn:gn)
414  !-----------------------------------------------------------------------------
415  dpdtx = - ddxdn1(ds,fxx) - ddyup1(ds,fxy) - ddzup1(ds,fzx) ! ops: 2 0 0
416  dpdty = - ddxup1(ds,fxy) - ddydn1(ds,fyy) - ddzup1(ds,fyz) ! ops: 2 0 0
417  dpdtz = - ddxup1(ds,fzx) - ddyup1(ds,fyz) - ddzdn1(ds,fzz) ! ops: 2 0 0
418  !-----------------------------------------------------------------------------
419  ! External force; e.g. point source gravity. Note that the density factor is
420  ! applied here, and should not be included in the function definition
421  !-----------------------------------------------------------------------------
422  if (allocated(self%force_per_unit_mass)) then
423  dpdt = dpdt + dd*self%force_per_unit_mass ! ops: 2 2 0
424  end if
425  if (allocated(self%force_per_unit_volume)) then
426  dpdt = dpdt + self%force_per_unit_volume
427  end if
428  !-----------------------------------------------------------------------------
429  ! Conservation of energy; reuse the mass flux for the energy flux
430  !-----------------------------------------------------------------------------
431  if (self%gamma == 1d0) then ! unless isothermal
432  dedt = dddt*ee ! ops: 0 1 0
433  else
434  block
435  real:: afm
436  afm = self%faver(fm(:,:,:,1))
437  fm = fm*down(ee) - (self%nu(5)*dd*down(du))*ddown(ds,ee) ! ops: 1 4 0
438  dedt = -div(ds,fm) - pg*divu + q
439  if (verbose==1) print '(a,10f12.6)', 'average(Q):', &
440  self%time, self%time*self%faver(q), &
441  self%faver(e-d**self%gamma/(self%gamma-1.0)), &
442  self%faver(- pg*divu), &
443  afm, self%faver(fm(:,:,:,1))
444  end block
445  end if
446  call deallocate_scalars_a (txx, tyy, tzz)
447  !-----------------------------------------------------------------------------
448  ! Induction equation
449  !-----------------------------------------------------------------------------
450  if (self%mhd) then
451  !---------------------------------------------------------------------------
452  ! Electric current, resistive electric field and magnetic dissipation
453  !---------------------------------------------------------------------------
454  j = curl_dn(ds,b)
455  self%maxJ=maxval(abs(j))
456  emf = self%nu(6)*edge(du)*j ! mhd: 0 6 0
457  q = yup(zup(jx*ex)) + zup(xup(jy*ey)) + xup(yup(jz*ez)) ! mhd: 2 3 0
458  if (self%gamma/=1.0) &
459  dedt = dedt + q ! mhd: 1 0 0
460  !---------------------------------------------------------------------------
461  ! Lorentz force, used instead of Maxwell stress.
462  !---------------------------------------------------------------------------
463  if (.not.do_maxwell) then
464  ! borrowing fx for scratch cell centering Bx -> fxx, By -> fyy, Bz -> fzz
465  fxx = xup(bx)
466  fyy = yup(by)
467  fzz = zup(bz)
468  dpdtx = dpdtx + zup(jy)*xdn(fzz) - yup(jz)*xdn(fyy)
469  dpdty = dpdty + xup(jz)*ydn(fxx) - zup(jx)*ydn(fzz)
470  dpdtz = dpdtz + yup(jx)*zdn(fyy) - xup(jy)*zdn(fxx)
471  end if
472  !---------------------------------------------------------------------------
473  ! Advective electric field on cell edges
474  !---------------------------------------------------------------------------
475  ex = ex + ydn(uz)*zdn(by) - zdn(uy)*ydn(bz) ! mhd: 2 2 0
476  ey = ey + zdn(ux)*xdn(bz) - xdn(uz)*zdn(bx) ! mhd: 2 2 0
477  ez = ez + xdn(uy)*ydn(bx) - ydn(ux)*xdn(by) ! mhd: 2 2 0
478  call non_ideal%update (self%gn, minval(ds), d, jx, jy, jz, bx, by, bz, q, emf, self%u_max)
479  dbdt = -curl_up(ds,emf)
480  end if
481  !-----------------------------------------------------------------------------
482  ! External heating / cooling, if any; e.g. radiative or conductive
483  !-----------------------------------------------------------------------------
484  if (allocated(self%heating_per_unit_mass)) then
485  dedt = dedt + d*self%heating_per_unit_mass ! ops: 2 2 0
486  end if
487  if (allocated(self%heating_per_unit_volume)) then
488  dedt = dedt + self%heating_per_unit_volume
489  end if
490  !-----------------------------------------------------------------------------
491  ! Time step limitation for rate of change dE/dt = (de/dt-(e/d)*dddt)/d
492  !-----------------------------------------------------------------------------
493  if (self%gamma /= 1d0) then
494  u_max = self%cdtd*dsmax*self%fmaxval(abs((dedt-ee*dddt)/d),outer=.true.) ! ops: 1 1 1
495  self%u_max = max(self%u_max,u_max)
496  if (verbose>2) print 1, 'u_max(dedt)', u_max, self%u_max
497  end if
498  !-----------------------------------------------------------------------------
499  call deallocate_vectors_a (ld, dd, emf, j, fd, fx, fy, fz, fm)
500  call deallocate_scalars_a (lnd, pa, du, ee, cs, pp, pb, q, uu, divu)
501  if (.not.self%do_pebbles) then
502  call deallocate_vectors (u)
503  call deallocate_scalars (pg)
504  end if
505  call trace%end (itimer)
506  if (mpi%master .and. first_time) then
507  print *, stagger%count, ' stagger calls', stagger%flops/product(real(self%n)),' flops per interior point'
508  stagger%count=0
509  stagger%flops=0
510  first_time = .false.
511  end if
512 END SUBROUTINE pde
513 
514 !===============================================================================
515 !> This procedure is called by solver_t%update, before it calls mhd_t%update
516 !===============================================================================
517 SUBROUTINE pre_update (self)
518  class(mhd_t), target:: self
519  real, dimension(:,:,:), pointer:: d, e, lnd, ee
520  integer:: iz
521  !-----------------------------------------------------------------------------
522  ! Compute initial temperature
523  !-----------------------------------------------------------------------------
524  call trace%begin('mhd_t%pre_update')
525  if (do_temperature) then
526  if (.not.associated(self%tt)) &
527  allocate (self%tt(self%gn(1),self%gn(2),self%gn(3)))
528  allocate ( lnd(self%gn(1),self%gn(2),self%gn(3)))
529  allocate ( ee(self%gn(1),self%gn(2),self%gn(3)))
530  d => self%mem(:,:,:,self%idx%d,self%it,1)
531  e => self%mem(:,:,:,self%idx%e,self%it,1)
532  lnd = alog(d)
533  ee = e/d
534  call eos%lookup_table (shape(self%tt), lnd=lnd, ee=ee, tt=self%tt)
535  deallocate (lnd, ee)
536  if (verbose > 0) &
537  write(io%output,*) self%id, 'mhd_t%pre_update: tt =', &
538  minval(self%tt), maxval(self%tt)
539  call self%aux%register ('tt', self%tt)
540  end if
541  !-----------------------------------------------------------------------------
542  ! Now we can call the extras_t%pre_update, which may need the temperature
543  !-----------------------------------------------------------------------------
544  call self%extras_t%pre_update
545  call trace%end()
546 END SUBROUTINE pre_update
547 
548 !===============================================================================
549 !> Take a complete update step. This can be overloaded in higher levels, e.g.
550 !> in experiment_mod, where boundary conditions may be inserted
551 !===============================================================================
552 SUBROUTINE update (self)
553  class(mhd_t):: self
554  integer, save:: itimer=0
555  !-----------------------------------------------------------------------------
556  call trace%begin('mhd_t%update',itimer=itimer)
557  !-----------------------------------------------------------------------------
558  ! Output the state after guard zone downloads, selfgravity, and boundary
559  ! conditions. If another call is placed at an earlier point in the task handling
560  ! the call here will have no effect, since out_next has already been updated.
561  !-----------------------------------------------------------------------------
562  call self%output
563  if (self%gamma==1d0) then
564  self%mem(:,:,:,self%idx%e,self%it,1) = &
565  self%mem(:,:,:,self%idx%d,self%it,1)*self%csound**2
566  end if
567  if (self%mhd) call self%divB_clean
568  call self%pde
569  !-----------------------------------------------------------------------------
570  ! The call to courant_condition may optionally be interscepted in extras_t
571  !-----------------------------------------------------------------------------
572  call self%courant_condition
573  if (verbose>2) &
574  call self%check_density('before timestep')
575  call timestep%update (self%id, self%iit, self%dt, self%mem, self%mesh, self%lock)
576  call self%check_density(' after timestep')
577  call self%counter_update
578  !-----------------------------------------------------------------------------
579  call trace%end (itimer)
580 END SUBROUTINE update
581 
582 !===============================================================================
583 !> Perform a single iteration of a div(B) clean. The limit of stability is
584 !> eps=1./6. This is intended to prevent build-up of layers of large div(B)
585 !> between patches. An exact solution would be to solve an internal Laplace
586 !> equation, with boundary conditiosn = the boundary value of div(B), and then
587 !> subtract the gradient of the resulting potential.
588 !>
589 !> The problem can also be greatly reduced by using a higher order time
590 !> interpoation scheme, thus making the electric field computed from the
591 !> guard zone values of B and U a better match to the neighbors E-fields.
592 !===============================================================================
593 SUBROUTINE divb_clean1(self)
594  class(mhd_t):: self
595  !.............................................................................
596  real, dimension(:,:,:,:), pointer:: b
597  real, dimension(:,:,:) , pointer:: bx, by, bz, divb
598  real, dimension(:,:,:,:) , allocatable:: dndivb
599  real, parameter:: eps=0.11
600  !-----------------------------------------------------------------------------
601  b => self%mem(:,:,:,self%idx%bx:self%idx%bz,self%it,1)
602  bx => b(:,:,:,1)
603  by => b(:,:,:,2)
604  bz => b(:,:,:,3)
605  call allocate_scalars (self%gn, divb)
606  call allocate_vectors_a (self%gn, dndivb)
607  divb = ddup(bx,1) + ddup(by,2) + ddup(bz,3)
608  bx = bx + eps*dddn(divb,1)
609  by = by + eps*dddn(divb,2)
610  bz = bz + eps*dddn(divb,3)
611  call deallocate_vectors_a (dndivb)
612  call deallocate_scalars(divb)
613  return
614 contains
615  !-----------------------------------------------------------------------------
616  function ddup(f,i) result (g)
617  real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
618  integer:: i
619  !.............................................................................
620  g = cshift(f,1,i) - f
621  end function
622  !-----------------------------------------------------------------------------
623  function dddn(f,i) result (g)
624  real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
625  integer:: i
626  !.............................................................................
627  g = f - cshift(f,-1,i)
628  end function
629 END SUBROUTINE divb_clean1
630 
631 !===============================================================================
632 !> Enforce div(B)=0 in the guard zones, by succesively imposing the condition,
633 !> layer-by-layer, on the component perpendicular to each of the patch faces.
634 !===============================================================================
635 SUBROUTINE divb_clean2 (self)
636  class(mhd_t):: self
637  !.............................................................................
638  real, dimension(:,:,:), pointer:: bx, by, bz
639  real, dimension(:,:,:), allocatable:: divb
640  real, parameter:: eps=0.11
641  integer:: l(3), u(3), ix, iy, iz, i1, j1, k1
642  integer:: itimer=0
643  !-----------------------------------------------------------------------------
644  call trace%begin('mhd_t%divb_clean2', itimer=itimer)
645  bx => self%mem(:,:,:,self%idx%bx,self%it,1)
646  by => self%mem(:,:,:,self%idx%by,self%it,1)
647  bz => self%mem(:,:,:,self%idx%bz,self%it,1)
648  l = self%mesh%lb
649  u = self%mesh%ub-1
650  i1 = merge(1, 0, self%n(1) > 1)
651  j1 = merge(1, 0, self%n(2) > 1)
652  k1 = merge(1, 0, self%n(3) > 1)
653  !--- x-direction -------------------------------------------------------------
654  do ix=self%mesh(1)%lo,self%mesh(1)%lb,-1
655  bx(ix ,l(2):u(2),l(3):u(3)) = &
656  bx(ix+i1,l(2):u(2),l(3):u(3)) + ( &
657  by(ix,l(2)+j1:u(2)+j1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
658  bz(ix,l(2) :u(2) ,l(3)+k1:u(3)+k1) - bz(ix,l(2):u(2),l(3):u(3)))
659  end do
660  do ix=self%mesh(1)%ui,self%mesh(1)%ub-1
661  bx(ix+i1,l(2):u(2),l(3):u(3)) = &
662  bx(ix ,l(2):u(2),l(3):u(3)) - ( &
663  by(ix,l(2)+j1:u(2)+j1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
664  bz(ix,l(2) :u(2) ,l(3)+k1:u(3)+k1) - bz(ix,l(2):u(2),l(3):u(3)))
665  end do
666  !--- y-direction -------------------------------------------------------------
667  do iy=self%mesh(2)%lo,self%mesh(2)%lb,-1
668  by(l(1):u(1),iy ,l(3):u(3)) = &
669  by(l(1):u(1),iy+j1,l(3):u(3)) + ( &
670  bx(l(1)+i1:u(1)+i1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
671  bz(l(1) :u(1) ,iy,l(3)+k1:u(3)+k1) - bz(l(1):u(1),iy,l(3):u(3)))
672  end do
673  do iy=self%mesh(2)%ui,self%mesh(2)%ub-1
674  by(l(1):u(1),iy+j1,l(3):u(3)) = &
675  by(l(1):u(1),iy ,l(3):u(3)) - ( &
676  bx(l(1)+i1:u(1)+i1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
677  bz(l(1) :u(1) ,iy,l(3)+k1:u(3)+k1) - bz(l(1):u(1),iy,l(3):u(3)))
678  end do
679  !--- z-direction -------------------------------------------------------------
680  do iz=self%mesh(3)%lo,self%mesh(3)%lb,-1
681  bz(l(1):u(1),l(2):u(2),iz ) = &
682  bz(l(1):u(1),l(2):u(2),iz+k1) + ( &
683  bx(l(1)+i1:u(1)+i1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
684  by(l(1) :u(1) ,l(2)+j1:u(2)+j1,iz) - by(l(1):u(1),l(2):u(2),iz))
685  end do
686  do iz=self%mesh(3)%ui,self%mesh(3)%ub-1
687  bz(l(1):u(1),l(2):u(2),iz+k1) = &
688  bz(l(1):u(1),l(2):u(2),iz ) - ( &
689  bx(l(1)+i1:u(1)+i1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
690  by(l(1) :u(1) ,l(2)+j1:u(2)+j1,iz) - by(l(1):u(1),l(2):u(2),iz))
691  end do
692  call trace%end (itimer)
693 END SUBROUTINE divb_clean2
694 
695 !===============================================================================
696 !> Add extra, MHD-related information to the patch .txt output file
697 !===============================================================================
698 SUBROUTINE output (self)
699  class(mhd_t):: self
700  character(len=64):: filename
701  real(8):: out_next
702  logical:: exists
703  !-----------------------------------------------------------------------------
704  call trace%begin ('mhd_t%output')
705  out_next = self%out_next
706  call self%extras_t%output
707  if (self%time >= out_next .and. io%do_output .and. io%do_legacy) then
708  !$omp critical (out_cr)
709  write (filename,'(a,i4.4,"_",i4.4,".txt")') trim(io%outputname), &
710  self%id, self%iout
711  inquire(file=filename, exist=exists)
712  if (exists) then
713  open (io%data_unit, file=trim(filename), form='formatted', status='old', &
714  position='append')
715  write (io%data_unit,'(8f10.4)') self%nu
716  close (io%data_unit)
717  end if
718  !$omp end critical (out_cr)
719  end if
720  call trace%end()
721 END SUBROUTINE output
722 
723 !===============================================================================
724 FUNCTION gas_pressure (self) result(pg)
725  class(mhd_t):: self
726  real, dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
727  real, dimension(:,:,:), pointer:: d, lnd, e, ee
728  integer, save:: nprint=1
729  !-----------------------------------------------------------------------------
730  call trace%begin ('mhd_t%gas_pressure_sub')
731  d => self%mem(:,:,:,self%idx%d ,self%it,1)
732  e => self%mem(:,:,:,self%idx%e ,self%it,1)
733  call allocate_scalars(self%gn,ee)
734  ee = e/d
735  if (trim(self%eos)=='ideal') then
736  if (self%gamma==1d0) then
737  pg = d*self%csound**2
738  else
739  pg = d*ee*(self%gamma-1)
740  end if
741  else
742  if (do_temperature) then
743  call eos%lookup_table (shape(d), d=d, ee=ee, pg=pg, tt=self%tt)
744  else
745  call eos%lookup_table (shape(d), d=d, ee=ee, pg=pg)
746  end if
747  if (verbose > 1 .or. nprint > 0) then
748  nprint = nprint-1
749  write(io%output,*) self%id, 'lookup_table: lnd =', minval(lnd), maxval(lnd)
750  write(io%output,*) self%id, 'lookup_table: ee =', minval(ee) , maxval(ee)
751  write(io%output,*) self%id, 'lookup_table: pg =', minval(pg) , maxval(pg)
752  if (do_temperature) &
753  write(io%output,*) self%id, 'lookup_table: tt =', minval(self%tt), &
754  maxval(self%tt)
755  end if
756  end if
757  call deallocate_scalars(ee)
758  call trace%end()
759 END FUNCTION gas_pressure
760 
761 !===============================================================================
762 SUBROUTINE gas_pressure_sub (self, d, lnd, ee, pg)
763  class(mhd_t):: self
764  real, dimension(:,:,:):: pg
765  real, dimension(:,:,:), pointer:: d, lnd, ee
766  integer, save:: nprint=1
767  !-----------------------------------------------------------------------------
768  call trace%begin ('mhd_t%gas_pressure_sub')
769  if (trim(self%eos)=='ideal') then
770  if (self%gamma==1d0) then
771  pg = d*self%csound**2
772  else
773  pg = d*ee*(self%gamma-1)
774  end if
775  else
776  if (do_temperature) then
777  call eos%lookup_table (shape(lnd), lnd=lnd, ee=ee, pg=pg, tt=self%tt)
778  else
779  call eos%lookup_table (shape(lnd), lnd=lnd, ee=ee, pg=pg)
780  end if
781  if (verbose > 1 .or. nprint > 0) then
782  nprint = nprint-1
783  write(io%output,*) self%id, 'lookup_table: lnd =', minval(lnd), maxval(lnd)
784  write(io%output,*) self%id, 'lookup_table: ee =', minval(ee) , maxval(ee)
785  write(io%output,*) self%id, 'lookup_table: pg =', minval(pg) , maxval(pg)
786  if (do_temperature) &
787  write(io%output,*) self%id, 'lookup_table: tt =', minval(self%tt), &
788  maxval(self%tt)
789  end if
790  end if
791  call trace%end()
792 END SUBROUTINE gas_pressure_sub
793 
794 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
Fundamental constants in CGS and SI units.
Definition: math_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
Module with which one can register any number of pointers to real or integer arrays, and then output the contents of the links later.
Definition: aux_mod.f90:14
6th order stagger operators, with self-test procedure
Definition: stagger_2nd.f90:4