MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_mhd_phys.t
Go to the documentation of this file.
1 !> Magneto-hydrodynamics module
3  use mod_global_parameters, only: std_len
4  implicit none
5  private
6 
7  !> Whether an energy equation is used
8  logical, public, protected :: mhd_energy = .true.
9 
10  !> Whether thermal conduction is used
11  logical, public, protected :: mhd_thermal_conduction = .false.
12 
13  !> Whether radiative cooling is added
14  logical, public, protected :: mhd_radiative_cooling = .false.
15 
16  !> Whether viscosity is added
17  logical, public, protected :: mhd_viscosity = .false.
18 
19  !> Whether gravity is added
20  logical, public, protected :: mhd_gravity = .false.
21 
22  !> Whether Hall-MHD is used
23  logical, public, protected :: mhd_hall = .false.
24 
25  !> Whether particles module is added
26  logical, public, protected :: mhd_particles = .false.
27 
28  !> Whether magnetofriction is added
29  logical, public, protected :: mhd_magnetofriction = .false.
30 
31  !> Whether GLM-MHD is used
32  logical, public, protected :: mhd_glm = .false.
33 
34  !> Whether divB cleaning sources are added splitting from fluid solver
35  logical, public, protected :: source_split_divb = .false.
36 
37  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
38  !> taking values within [0, 1]
39  double precision, public :: mhd_glm_alpha = 0.5d0
40 
41  !> Use Boris approximation
42  character(len=20) :: mhd_boris_method = "none"
43 
44  integer, parameter :: boris_none = 0
45  integer, parameter :: boris_reduced_force = 1
46  integer, parameter :: boris_simplification = 2
47  integer :: mhd_boris_type = boris_none
48 
49  !> Speed of light for Boris' approximation. If negative, test changes to the
50  !> momentum equation with gamma_A = 1
51  double precision :: mhd_boris_c = 0.0d0
52 
53  !> MHD fourth order
54  logical, public, protected :: mhd_4th_order = .false.
55 
56  !> Number of tracer species
57  integer, public, protected :: mhd_n_tracer = 0
58 
59  !> Index of the density (in the w array)
60  integer, public, protected :: rho_
61 
62  !> Indices of the momentum density
63  integer, allocatable, public, protected :: mom(:)
64 
65  !> Index of the energy density (-1 if not present)
66  integer, public, protected :: e_
67 
68  !> Index of the gas pressure (-1 if not present) should equal e_
69  integer, public, protected :: p_
70 
71  !> Indices of the magnetic field
72  integer, allocatable, public, protected :: mag(:)
73 
74  !> Indices of the GLM psi
75  integer, public, protected :: psi_
76 
77  !> Indices of the tracers
78  integer, allocatable, public, protected :: tracer(:)
79 
80  !> The adiabatic index
81  double precision, public :: mhd_gamma = 5.d0/3.0d0
82 
83  !> The adiabatic constant
84  double precision, public :: mhd_adiab = 1.0d0
85 
86  !> The MHD resistivity
87  double precision, public :: mhd_eta = 0.0d0
88 
89  !> The MHD hyper-resistivity
90  double precision, public :: mhd_eta_hyper = 0.0d0
91 
92  !> TODO: what is this?
93  double precision, public :: mhd_etah = 0.0d0
94 
95  !> The small_est allowed energy
96  double precision, protected :: small_e
97 
98  !> The number of waves
99  integer :: nwwave=8
100 
101  !> Method type to clean divergence of B
102  character(len=std_len), public, protected :: typedivbfix = 'linde'
103 
104  !> Method type of constrained transport
105  character(len=std_len), public, protected :: type_ct = 'uct_contact'
106 
107  !> Whether divB is computed with a fourth order approximation
108  logical, public, protected :: mhd_divb_4thorder = .false.
109 
110  !> Method type in a integer for good performance
111  integer :: type_divb
112 
113  !> Coefficient of diffusive divB cleaning
114  double precision :: divbdiff = 0.8d0
115 
116  !> Update all equations due to divB cleaning
117  character(len=std_len) :: typedivbdiff = 'all'
118 
119  !> Use a compact way to add resistivity
120  logical :: compactres = .false.
121 
122  !> Add divB wave in Roe solver
123  logical, public :: divbwave = .true.
124 
125  !> Helium abundance over Hydrogen
126  double precision, public, protected :: he_abundance=0.1d0
127 
128  !> To control divB=0 fix for boundary
129  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
130 
131  !> To skip * layer of ghost cells during divB=0 fix for boundary
132  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
133 
134  !> B0 field is force-free
135  logical, public, protected :: b0field_forcefree=.true.
136 
137  !> gamma minus one and its inverse
138  double precision :: gamma_1, inv_gamma_1
139 
140  ! DivB cleaning methods
141  integer, parameter :: divb_none = 0
142  integer, parameter :: divb_multigrid = -1
143  integer, parameter :: divb_glm = 1
144  integer, parameter :: divb_powel = 2
145  integer, parameter :: divb_janhunen = 3
146  integer, parameter :: divb_linde = 4
147  integer, parameter :: divb_lindejanhunen = 5
148  integer, parameter :: divb_lindepowel = 6
149  integer, parameter :: divb_lindeglm = 7
150  integer, parameter :: divb_ct = 8
151 
152 
153  ! Public methods
154  public :: mhd_phys_init
155  public :: mhd_kin_en
156  public :: mhd_get_pthermal
157  public :: mhd_get_v
158  public :: mhd_get_v_idim
159  public :: mhd_to_conserved
160  public :: mhd_to_primitive
161  public :: mhd_get_csound2
162  public :: mhd_face_to_center
163  public :: get_divb
164  public :: get_current
165  public :: get_normalized_divb
166  public :: b_from_vector_potential
167  public :: mhd_mag_en_all
168  {^nooned
169  public :: mhd_clean_divb_multigrid
170  }
171 
172 contains
173 
174  !> Read this module"s parameters from a file
175  subroutine mhd_read_params(files)
177  character(len=*), intent(in) :: files(:)
178  integer :: n
179 
180  namelist /mhd_list/ mhd_energy, mhd_n_tracer, mhd_gamma, mhd_adiab,&
184  typedivbdiff, type_ct, compactres, divbwave, he_abundance, si_unit, b0field,&
187  mhd_boris_method, mhd_boris_c
188 
189  do n = 1, size(files)
190  open(unitpar, file=trim(files(n)), status="old")
191  read(unitpar, mhd_list, end=111)
192 111 close(unitpar)
193  end do
194 
195  end subroutine mhd_read_params
196 
197  !> Write this module's parameters to a snapsoht
198  subroutine mhd_write_info(fh)
200  integer, intent(in) :: fh
201  integer, parameter :: n_par = 1
202  double precision :: values(n_par)
203  character(len=name_len) :: names(n_par)
204  integer, dimension(MPI_STATUS_SIZE) :: st
205  integer :: er
206 
207  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
208 
209  names(1) = "gamma"
210  values(1) = mhd_gamma
211  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
212  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
213  end subroutine mhd_write_info
214 
215  subroutine mhd_angmomfix(fC,x,wnew,ixI^L,ixO^L,idim)
217  double precision, intent(in) :: x(ixi^s,1:ndim)
218  double precision, intent(inout) :: fC(ixi^s,1:nwflux,1:ndim), wnew(ixi^s,1:nw)
219  integer, intent(in) :: ixI^L, ixO^L
220  integer, intent(in) :: idim
221  integer :: hxO^L, kxC^L, iw
222  double precision :: inv_volume(ixi^s)
223 
224  call mpistop("to do")
225 
226  end subroutine mhd_angmomfix
227 
228  subroutine mhd_phys_init()
232  use mod_viscosity, only: viscosity_init
233  use mod_gravity, only: gravity_init
234  use mod_particles, only: particles_init
236  use mod_physics
237  {^nooned
239  }
240 
241  integer :: itr, idir
242 
243  call mhd_read_params(par_files)
244 
245  physics_type = "mhd"
247  ! set default gamma for polytropic/isothermal process
248  if(.not.mhd_energy) mhd_gamma=1.d0
250  if(ndim==1) typedivbfix='none'
251  select case (typedivbfix)
252  case ('none')
253  type_divb = divb_none
254  {^nooned
255  case ('multigrid')
256  type_divb = divb_multigrid
257  use_multigrid = .true.
258  mg%operator_type = mg_laplacian
260  }
261  case ('glm')
262  mhd_glm = .true.
263  need_global_cmax = .true.
264  type_divb = divb_glm
265  case ('powel', 'powell')
266  type_divb = divb_powel
267  case ('janhunen')
268  type_divb = divb_janhunen
269  case ('linde')
270  type_divb = divb_linde
271  case ('lindejanhunen')
272  type_divb = divb_lindejanhunen
273  case ('lindepowel')
274  type_divb = divb_lindepowel
275  case ('lindeglm')
276  mhd_glm = .true.
277  need_global_cmax = .true.
278  type_divb = divb_lindeglm
279  case ('ct')
280  type_divb = divb_ct
281  stagger_grid = .true.
282  case default
283  call mpistop('Unknown divB fix')
284  end select
285 
286  select case (mhd_boris_method)
287  case ("none")
288  mhd_boris_type = boris_none
289  case ("reduced_force")
290  mhd_boris_type = boris_reduced_force
291  case ("simplification")
292  mhd_boris_type = boris_simplification
293  case default
294  call mpistop("Unknown mhd_boris_method (none, reduced_force, simplification)")
295  end select
296 
297  ! Determine flux variables
298  rho_ = var_set_rho()
299 
300  allocate(mom(ndir))
301  mom(:) = var_set_momentum(ndir)
302 
303  ! Set index of energy variable
304  if (mhd_energy) then
305  nwwave = 8
306  e_ = var_set_energy() ! energy density
307  p_ = e_ ! gas pressure
308  else
309  nwwave = 7
310  e_ = -1
311  p_ = -1
312  end if
313 
314  allocate(mag(ndir))
315  mag(:) = var_set_bfield(ndir)
316 
317  if (mhd_glm) then
318  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
319  else
320  psi_ = -1
321  end if
322 
323  allocate(tracer(mhd_n_tracer))
324 
325  ! Set starting index of tracers
326  do itr = 1, mhd_n_tracer
327  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
328  end do
329 
330  ! determine number of stagger variables
331  if(stagger_grid) nws=ndim
332 
333  nvector = 2 ! No. vector vars
334  allocate(iw_vector(nvector))
335  iw_vector(1) = mom(1) - 1 ! TODO: why like this?
336  iw_vector(2) = mag(1) - 1 ! TODO: why like this?
337 
338  ! Check whether custom flux types have been defined
339  if (.not. allocated(flux_type)) then
340  allocate(flux_type(ndir, nw))
342  else if (any(shape(flux_type) /= [ndir, nw])) then
343  call mpistop("phys_check error: flux_type has wrong shape")
344  end if
345  do idir=1,ndir
346  if(ndim>1) flux_type(idir,mag(idir))=flux_tvdlf
347  end do
348  if(mhd_glm .and. ndim>1) flux_type(:,psi_)=flux_tvdlf
349 
350  select case (mhd_boris_method)
351  case ("none")
352  mhd_boris_type = boris_none
353  case ("reduced_force")
354  mhd_boris_type = boris_reduced_force
355  case ("simplification")
356  mhd_boris_type = boris_simplification
357  do idir = 1, ndir
358  phys_iw_methods(mom(idir))%inv_capacity => mhd_gamma2_alfven
359  end do
360  case default
361  call mpistop("Unknown mhd_boris_method (none, reduced_force, simplification)")
362  end select
363 
379 
380  if(type_divb==divb_glm) then
382  end if
383 
384  ! if using ct stagger grid, boundary divb=0 is not done here
385  if(stagger_grid) then
389  else if(ndim>1) then
391  end if
392 
393  ! Whether diagonal ghost cells are required for the physics
394  if(type_divb < divb_linde) phys_req_diagonal = .false.
395 
396  ! derive units from basic units
397  call mhd_physical_units()
398 
399  if(.not. mhd_energy .and. mhd_thermal_conduction) then
400  call mpistop("thermal conduction needs mhd_energy=T")
401  end if
402  if(.not. mhd_energy .and. mhd_radiative_cooling) then
403  call mpistop("radiative cooling needs mhd_energy=T")
404  end if
405 
406  ! initialize thermal conduction module
407  if (mhd_thermal_conduction) then
408  phys_req_diagonal = .true.
410  end if
411 
412  ! Initialize radiative cooling module
413  if (mhd_radiative_cooling) then
415  end if
416 
417  ! Initialize viscosity module
419 
420  ! Initialize gravity module
421  if(mhd_gravity) then
422  call gravity_init()
423  end if
424 
425  ! Initialize particles module
426  if(mhd_particles) then
427  call particles_init()
428  phys_req_diagonal = .true.
429  end if
430 
431  ! initialize magnetofriction module
432  if(mhd_magnetofriction) then
433  phys_req_diagonal = .true.
434  call magnetofriction_init()
435  end if
436 
437  ! For Hall, we need one more reconstructed layer since currents are computed
438  ! in getflux: assuming one additional ghost layer (two for FOURTHORDER) was
439  ! added in nghostcells.
440  if (mhd_hall) then
441  phys_req_diagonal = .true.
442  if (mhd_4th_order) then
444  else
446  end if
447  end if
448 
449  end subroutine mhd_phys_init
450 
451  subroutine mhd_check_params
453 
454  ! after user parameter setting
455  gamma_1=mhd_gamma-1.d0
456 
457  if (.not. mhd_energy) then
458  if (mhd_gamma <= 0.0d0) call mpistop ("Error: mhd_gamma <= 0")
459  if (mhd_adiab < 0.0d0) call mpistop ("Error: mhd_adiab < 0")
461  else
462  if (mhd_gamma <= 0.0d0 .or. mhd_gamma == 1.0d0) &
463  call mpistop ("Error: mhd_gamma <= 0 or mhd_gamma == 1")
464  inv_gamma_1=1.d0/gamma_1
465  small_e = small_pressure * inv_gamma_1
466  end if
467 
468  if (mhd_boris_type > 0 .and. abs(mhd_boris_c) <= 0.0d0) then
469  call mpistop("You have not specified mhd_boris_c")
470  end if
471 
472  end subroutine mhd_check_params
473 
474  subroutine mhd_physical_units()
476  double precision :: mp,kB,miu0
477  ! Derive scaling units
478  if(si_unit) then
479  mp=mp_si
480  kb=kb_si
481  miu0=miu0_si
482  else
483  mp=mp_cgs
484  kb=kb_cgs
485  miu0=4.d0*dpi
486  end if
487  if(unit_velocity==0) then
493  else
499  end if
500 
501  end subroutine mhd_physical_units
502 
503  subroutine mhd_check_w(primitive,ixI^L,ixO^L,w,flag)
505 
506  logical, intent(in) :: primitive
507  integer, intent(in) :: ixI^L, ixO^L
508  double precision, intent(in) :: w(ixi^s,nw)
509  integer, intent(inout) :: flag(ixi^s)
510  double precision :: tmp(ixi^s)
511 
512  flag(ixo^s)=0
513  where(w(ixo^s, rho_) < small_density) flag(ixo^s) = rho_
514 
515  if (mhd_energy) then
516  if (block%e_is_internal) then
517  where(w(ixo^s, e_) < small_pressure*inv_gamma_1) flag(ixo^s) = e_
518  else
519  if (primitive)then
520  where(w(ixo^s, e_) < small_pressure) flag(ixo^s) = e_
521  else
522  ! Calculate pressure=(gamma-1)*(e-0.5*(2ek+2eb))
523  tmp(ixo^s)=w(ixo^s,e_) - &
524  mhd_kin_en(w,ixi^l,ixo^l)-mhd_mag_en(w,ixi^l,ixo^l)
525  tmp(ixo^s)=gamma_1*tmp(ixo^s)
526  where(tmp(ixo^s) < small_pressure) flag(ixo^s) = e_
527  end if
528  end if
529  end if
530  end subroutine mhd_check_w
531 
532  !> Transform primitive variables into conservative ones
533  subroutine mhd_to_conserved(ixI^L,ixO^L,w,x)
535  integer, intent(in) :: ixI^L, ixO^L
536  double precision, intent(inout) :: w(ixi^s, nw)
537  double precision, intent(in) :: x(ixi^s, 1:ndim)
538  integer :: idir, itr
539 
541  call mhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'mhd_to_conserved')
542  end if
543 
544  if (mhd_energy) then
545  ! Calculate total energy from pressure, kinetic and magnetic energy
546  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
547  if(.not.block%e_is_internal) w(ixo^s,e_)=w(ixo^s,e_) + &
548  0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * w(ixo^s, rho_) + &
549  mhd_mag_en(w, ixi^l, ixo^l)
550  end if
551 
552  ! Convert velocity to momentum
553  do idir = 1, ndir
554  w(ixo^s, mom(idir)) = w(ixo^s, rho_) * w(ixo^s, mom(idir))
555  end do
556 
557  if (check_small_values .and. .not. small_values_use_primitive) then
558  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_conserved')
559  end if
560  end subroutine mhd_to_conserved
561 
562  !> Transform conservative variables into primitive ones
563  subroutine mhd_to_primitive(ixI^L,ixO^L,w,x)
565  integer, intent(in) :: ixI^L, ixO^L
566  double precision, intent(inout) :: w(ixi^s, nw)
567  double precision, intent(in) :: x(ixi^s, 1:ndim)
568  double precision :: inv_rho(ixo^s)
569  integer :: itr, idir
570 
571  if (check_small_values .and. .not. small_values_use_primitive) then
572  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive')
573  end if
574 
575  inv_rho = 1.0d0 / w(ixo^s, rho_)
576 
577  if (mhd_energy) then
578  ! Calculate pressure = (gamma-1) * (e-ek-eb)
579  if(.not.block%e_is_internal) then
580  w(ixo^s, p_) = w(ixo^s, e_) &
581  - mhd_kin_en(w, ixi^l, ixo^l, inv_rho) &
582  - mhd_mag_en(w, ixi^l, ixo^l)
583  end if
584  w(ixo^s, p_) = gamma_1*w(ixo^s, p_)
585  end if
586 
587  ! Convert momentum to velocity
588  do idir = 1, ndir
589  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
590  end do
591 
593  call mhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'mhd_to_primitive')
594  end if
595  end subroutine mhd_to_primitive
596 
597  subroutine mhd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
599  use mod_small_values
600  logical, intent(in) :: primitive
601  integer, intent(in) :: ixI^L,ixO^L
602  double precision, intent(inout) :: w(ixi^s,1:nw)
603  double precision, intent(in) :: x(ixi^s,1:ndim)
604  character(len=*), intent(in) :: subname
605 
606  double precision :: smallone
607  integer :: idir, flag(ixi^s)
608 
609  if (small_values_method == "ignore") return
610 
611  call mhd_check_w(primitive, ixi^l, ixo^l, w, flag)
612 
613  if (any(flag(ixo^s) /= 0)) then
614  select case (small_values_method)
615  case ("replace")
616  if (small_values_fix_iw(rho_)) then
617  where(flag(ixo^s) /= 0) w(ixo^s,rho_) = small_density
618  end if
619 
620  do idir = 1, ndir
621  if (small_values_fix_iw(mom(idir))) then
622  where(flag(ixo^s) /= 0) w(ixo^s, mom(idir)) = 0.0d0
623  end if
624  end do
625 
626  if (mhd_energy) then
627  if (small_values_fix_iw(e_)) then
628  if(primitive) then
629  where(flag(ixo^s) /= 0) w(ixo^s,e_) = small_pressure
630  else
631  where(flag(ixo^s) /= 0)
632  w(ixo^s,e_) = small_e + 0.5d0 * &
633  sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_) + &
634  mhd_mag_en(w, ixi^l, ixo^l)
635  end where
636  end if
637  end if
638  end if
639  case ("average")
640  call small_values_average(ixi^l, ixo^l, w, x, flag)
641  case default
642  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
643  end select
644  end if
645  end subroutine mhd_handle_small_values
646 
647  !> Convert energy to entropy
648  subroutine e_to_rhos(ixI^L,ixO^L,w,x)
650  integer, intent(in) :: ixI^L, ixO^L
651  double precision,intent(inout) :: w(ixi^s,nw)
652  double precision, intent(in) :: x(ixi^s,1:ndim)
653 
654  if (mhd_energy) then
655  if(.not.block%e_is_internal) &
656  w(ixo^s, e_) = w(ixo^s, e_) - mhd_kin_en(w, ixi^l, ixo^l) &
657  - mhd_mag_en(w, ixi^l, ixo^l)
658  w(ixo^s, e_) = gamma_1* w(ixo^s, rho_)**(1.0d0 - mhd_gamma) * &
659  w(ixo^s, e_)
660  else
661  call mpistop("e_to_rhos can not be used without energy equation!")
662  end if
663  end subroutine e_to_rhos
664 
665  !> Convert entropy to energy
666  subroutine rhos_to_e(ixI^L,ixO^L,w,x)
668  integer, intent(in) :: ixI^L, ixO^L
669  double precision :: w(ixi^s,nw)
670  double precision, intent(in) :: x(ixi^s,1:ndim)
671 
672  if (mhd_energy) then
673  w(ixo^s, e_) = w(ixo^s, rho_)**gamma_1 * w(ixo^s, e_) &
674  * inv_gamma_1
675  if(.not.block%e_is_internal) &
676  w(ixo^s, e_) =w(ixo^s, e_) + mhd_kin_en(w, ixi^l, ixo^l) + &
677  mhd_mag_en(w, ixi^l, ixo^l)
678  else
679  call mpistop("rhos_to_e can not be used without energy equation!")
680  end if
681  end subroutine rhos_to_e
682 
683  !> Calculate v vector
684  subroutine mhd_get_v(w,x,ixI^L,ixO^L,v)
686 
687  integer, intent(in) :: ixI^L, ixO^L
688  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
689  double precision, intent(out) :: v(ixi^s,ndir)
690 
691  integer :: idir
692 
693  do idir=1,ndir
694  v(ixo^s,idir) = w(ixo^s, mom(idir)) / w(ixo^s,rho_)
695  end do
696 
697  end subroutine mhd_get_v
698 
699  !> Calculate v component
700  subroutine mhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
702 
703  integer, intent(in) :: ixI^L, ixO^L, idim
704  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
705  double precision, intent(out) :: v(ixi^s)
706 
707  v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s,rho_)
708 
709  end subroutine mhd_get_v_idim
710 
711  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
712  subroutine mhd_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
714 
715  integer, intent(in) :: ixI^L, ixO^L, idim
716  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
717  double precision, intent(inout) :: cmax(ixi^s)
718 
719  call mhd_get_csound(w,x,ixi^l,ixo^l,idim,cmax)
720 
721  cmax(ixo^s)=abs(w(ixo^s,mom(idim))/w(ixo^s,rho_))+cmax(ixo^s)
722 
723  end subroutine mhd_get_cmax
724 
725  !> Estimating bounds for the minimum and maximum signal velocities
726  subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,cmax,cmin)
729 
730  integer, intent(in) :: ixI^L, ixO^L, idim
731  double precision, intent(in) :: wLC(ixi^s, nw), wRC(ixi^s, nw)
732  double precision, intent(in) :: wLp(ixi^s, nw), wRp(ixi^s, nw)
733  double precision, intent(in) :: x(ixi^s,1:ndim)
734  double precision, intent(inout) :: cmax(ixi^s)
735  double precision, intent(inout), optional :: cmin(ixi^s)
736 
737  double precision :: wmean(ixi^s,nw)
738  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
739  integer :: idimE,idimN
740 
741  if (typeboundspeed=='cmaxmean') then
742  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
743  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
744  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
745  if(present(cmin)) then
746  cmax(ixo^s)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
747  cmin(ixo^s)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
748  else
749  cmax(ixo^s)=abs(tmp1(ixo^s))+csoundr(ixo^s)
750  end if
751  else
752  ! This implements formula (10.52) from "Riemann Solvers and Numerical
753  ! Methods for Fluid Dynamics" by Toro.
754  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
755  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
756  tmp3(ixo^s)=1.d0/(sqrt(wlp(ixo^s,rho_))+sqrt(wrp(ixo^s,rho_)))
757  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
758  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
759  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
760  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
761  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
762  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
763  dmean(ixo^s)=sqrt(dmean(ixo^s))
764  if(present(cmin)) then
765  cmin(ixo^s)=umean(ixo^s)-dmean(ixo^s)
766  cmax(ixo^s)=umean(ixo^s)+dmean(ixo^s)
767  else
768  cmax(ixo^s)=abs(umean(ixo^s))+dmean(ixo^s)
769  end if
770  end if
771 
772  if(stagger_grid) then
773  ! calculate velocities related to different UCT schemes
774  select case(type_ct)
775  case('average')
776  case('uct_contact')
777  if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
778  ! get average normal velocity at cell faces
779  vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom(idim))+wrp(ixo^s,mom(idim)))
780  case('uct_hll')
781  if(.not.allocated(vcts%vbarC)) then
782  allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
783  allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
784  end if
785  ! Store magnitude of characteristics
786  if(present(cmin)) then
787  vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
788  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
789  else
790  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
791  vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
792  end if
793 
794  idimn=mod(idim,ndir)+1 ! 'Next' direction
795  idime=mod(idim+1,ndir)+1 ! Electric field direction
796  ! Store velocities
797  vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom(idimn))
798  vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(idimn))
799  vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
800  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
801  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
802 
803  vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom(idime))
804  vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(idime))
805  vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
806  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
807  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
808  case default
809  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
810  end select
811  end if
812 
813  end subroutine mhd_get_cbounds
814 
815  !> Calculate fast magnetosonic wave speed
816  subroutine mhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
818 
819  integer, intent(in) :: ixI^L, ixO^L, idim
820  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
821  double precision, intent(out):: csound(ixi^s)
822  double precision :: cfast2(ixi^s), AvMinCs2(ixi^s), b2(ixi^s), kmax
823  double precision :: inv_rho(ixo^s), gamma2(ixo^s)
824 
825  inv_rho=1.d0/w(ixo^s,rho_)
826 
827  if (mhd_boris_type == boris_reduced_force) then
828  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma2)
829  else
830  gamma2 = 1.0d0
831  end if
832 
833  call mhd_get_csound2(w,x,ixi^l,ixo^l,csound)
834 
835  ! store |B|^2 in v
836  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l) * gamma2
837 
838  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
839  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
840  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
841  * inv_rho * gamma2
842 
843  where(avmincs2(ixo^s)<zero)
844  avmincs2(ixo^s)=zero
845  end where
846 
847  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
848 
849  if (.not. mhd_hall) then
850  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
851  if (mhd_boris_type == boris_simplification) then
852  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
853  end if
854  else
855  ! take the Hall velocity into account:
856  ! most simple estimate, high k limit:
857  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
858  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
859  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
860  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
861  end if
862 
863  end subroutine mhd_get_csound
864 
865  !> Calculate fast magnetosonic wave speed
866  subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
868 
869  integer, intent(in) :: ixI^L, ixO^L, idim
870  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
871  double precision, intent(out):: csound(ixi^s)
872  double precision :: cfast2(ixi^s), AvMinCs2(ixi^s), b2(ixi^s), kmax
873  double precision :: inv_rho(ixo^s), gamma_A2(ixo^s)
874 
875  inv_rho=1.d0/w(ixo^s,rho_)
876 
877  if (mhd_boris_type == boris_reduced_force) then
878  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a2)
879  else
880  gamma_a2 = 1.0d0
881  end if
882 
883  if(mhd_energy) then
884  csound(ixo^s)=mhd_gamma*w(ixo^s,p_)*inv_rho
885  else
886  csound(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
887  end if
888  ! store |B|^2 in v
889  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l) * gamma_a2
890  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
891  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
892  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
893  * inv_rho * gamma_a2
894 
895  where(avmincs2(ixo^s)<zero)
896  avmincs2(ixo^s)=zero
897  end where
898 
899  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
900 
901  if (.not. mhd_hall) then
902  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
903  if (mhd_boris_type == boris_simplification) then
904  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
905  end if
906  else
907  ! take the Hall velocity into account:
908  ! most simple estimate, high k limit:
909  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
910  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
911  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
912  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
913  end if
914 
915  end subroutine mhd_get_csound_prim
916 
917  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
918  subroutine mhd_get_pthermal(w,x,ixI^L,ixO^L,pth)
920 
921  integer, intent(in) :: ixI^L, ixO^L
922  double precision, intent(in) :: w(ixi^s,nw)
923  double precision, intent(in) :: x(ixi^s,1:ndim)
924  double precision, intent(out):: pth(ixi^s)
925 
926  if(mhd_energy) then
927  if(block%e_is_internal) then
928  pth(ixo^s)=gamma_1*w(ixo^s,e_)
929  else
930  pth(ixo^s)=gamma_1*(w(ixo^s,e_)&
931  - mhd_kin_en(w,ixi^l,ixo^l)&
932  - mhd_mag_en(w,ixi^l,ixo^l))
933  end if
934  else
935  pth(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
936  end if
937  end subroutine mhd_get_pthermal
938 
939  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
940  !> csound2=gamma*p/rho
941  subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
943  integer, intent(in) :: ixI^L, ixO^L
944  double precision, intent(in) :: w(ixi^s,nw)
945  double precision, intent(in) :: x(ixi^s,1:ndim)
946  double precision, intent(out) :: csound2(ixi^s)
947 
948  if(mhd_energy) then
949  call mhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
950  csound2(ixo^s)=mhd_gamma*csound2(ixo^s)/w(ixo^s,rho_)
951  else
952  csound2(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
953  end if
954  end subroutine mhd_get_csound2
955 
956  !> Calculate total pressure within ixO^L including magnetic pressure
957  subroutine mhd_get_p_total(w,x,ixI^L,ixO^L,p)
959 
960  integer, intent(in) :: ixI^L, ixO^L
961  double precision, intent(in) :: w(ixi^s,nw)
962  double precision, intent(in) :: x(ixi^s,1:ndim)
963  double precision, intent(out) :: p(ixi^s)
964 
965  call mhd_get_pthermal(w,x,ixi^l,ixo^l,p)
966 
967  p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
968 
969  end subroutine mhd_get_p_total
970 
971  !> Calculate fluxes within ixO^L.
972  subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
974  use mod_usr_methods
975  use mod_geometry
976 
977  integer, intent(in) :: ixI^L, ixO^L, idim
978  ! conservative w
979  double precision, intent(in) :: wC(ixi^s,nw)
980  ! primitive w
981  double precision, intent(in) :: w(ixi^s,nw)
982  double precision, intent(in) :: x(ixi^s,1:ndim)
983  double precision,intent(out) :: f(ixi^s,nwflux)
984 
985  double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
986  double precision, allocatable:: vHall(:^d&,:)
987  integer :: idirmin, iw, idir, jdir, kdir
988 
989  if (mhd_hall) then
990  allocate(vhall(ixi^s,1:ndir))
991  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
992  end if
993 
994  if(b0field) tmp(ixo^s)=sum(block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
995 
996  if(mhd_energy) then
997  pgas=w(ixo^s,p_)
998  else
999  pgas=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
1000  end if
1001 
1002  ptotal = pgas + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1003 
1004  ! Get flux of density
1005  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
1006 
1007  ! Get flux of tracer
1008  do iw=1,mhd_n_tracer
1009  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
1010  end do
1011 
1012  ! Get flux of momentum
1013  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
1014  if (mhd_boris_type == boris_reduced_force) then
1015  do idir=1,ndir
1016  if(idim==idir) then
1017  f(ixo^s,mom(idir)) = pgas(ixo^s)
1018  else
1019  f(ixo^s,mom(idir)) = 0.0d0
1020  end if
1021  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+w(ixo^s,mom(idim))*wc(ixo^s,mom(idir))
1022  end do
1023  else
1024  ! Normal case (no Boris approximation)
1025  do idir=1,ndir
1026  if(idim==idir) then
1027  f(ixo^s,mom(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
1028  if(b0field) f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+tmp(ixo^s)
1029  else
1030  f(ixo^s,mom(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
1031  end if
1032  if (b0field) then
1033  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))&
1034  -w(ixo^s,mag(idir))*block%B0(ixo^s,idim,idim)&
1035  -w(ixo^s,mag(idim))*block%B0(ixo^s,idir,idim)
1036  end if
1037  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+w(ixo^s,mom(idim))*wc(ixo^s,mom(idir))
1038  end do
1039  end if
1040 
1041  ! Get flux of energy
1042  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
1043  if (mhd_energy) then
1044  if (block%e_is_internal) then
1045  f(ixo^s,e_)=w(ixo^s,mom(idim))*w(ixo^s,p_)
1046  if (mhd_hall) then
1047  call mpistop("solve pthermal not designed for Hall MHD")
1048  endif
1049  else
1050  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_) + ptotal(ixo^s))- &
1051  w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
1052 
1053  if (b0field) then
1054  f(ixo^s,e_) = f(ixo^s,e_) &
1055  + w(ixo^s,mom(idim)) * tmp(ixo^s) &
1056  - sum(w(ixo^s,mom(:))*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
1057  end if
1058 
1059  if (mhd_hall) then
1060  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
1061  if (mhd_etah>zero) then
1062  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
1063  sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
1064  - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
1065  if (b0field) then
1066  f(ixo^s,e_) = f(ixo^s,e_) &
1067  + vhall(ixo^s,idim) * tmp(ixo^s) &
1068  - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
1069  end if
1070  end if
1071  end if
1072 
1073  end if
1074  end if
1075 
1076  ! compute flux of magnetic field
1077  ! f_i[b_k]=v_i*b_k-v_k*b_i
1078  do idir=1,ndir
1079  if (idim==idir) then
1080  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
1081  if (mhd_glm) then
1082  f(ixo^s,mag(idir))=w(ixo^s,psi_)
1083  else
1084  f(ixo^s,mag(idir))=zero
1085  end if
1086  else
1087  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
1088 
1089  if (b0field) then
1090  f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
1091  +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
1092  -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
1093  end if
1094 
1095  if (mhd_hall) then
1096  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
1097  if (mhd_etah>zero) then
1098  if (b0field) then
1099  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1100  - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+block%B0(ixo^s,idim,idim)) &
1101  + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+block%B0(ixo^s,idir,idim))
1102  else
1103  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1104  - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
1105  + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
1106  end if
1107  end if
1108  end if
1109 
1110  end if
1111  end do
1112 
1113  if (mhd_glm) then
1114  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
1115  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
1116  end if
1117 
1118  end subroutine mhd_get_flux
1119 
1120  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
1121  subroutine mhd_add_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
1125  use mod_gravity, only: gravity_add_source
1126 
1127  integer, intent(in) :: ixI^L, ixO^L
1128  double precision, intent(in) :: qdt
1129  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1130  double precision, intent(inout) :: w(ixi^s,1:nw)
1131  logical, intent(in) :: qsourcesplit
1132  logical, intent(inout) :: active
1133 
1134  if (.not. qsourcesplit) then
1135  ! Source for solving internal energy
1136  if (mhd_energy .and. block%e_is_internal) then
1137  active = .true.
1138  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x)
1139  endif
1140 
1141  ! Source for B0 splitting
1142  if (b0field) then
1143  active = .true.
1144  call add_source_b0split(qdt,ixi^l,ixo^l,wct,w,x)
1145  end if
1146 
1147  ! Sources for resistivity in eqs. for e, B1, B2 and B3
1148  if (abs(mhd_eta)>smalldouble)then
1149  active = .true.
1150  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
1151  end if
1152 
1153  if (mhd_eta_hyper>0.d0)then
1154  active = .true.
1155  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
1156  end if
1157  end if
1158 
1159  {^nooned
1160  if(.not.source_split_divb .and. .not.qsourcesplit .and. istep==nstep) then
1161  ! Sources related to div B
1162  select case (type_divb)
1163  case (divb_none)
1164  ! Do nothing
1165  case (divb_glm)
1166  active = .true.
1167  call add_source_glm(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1168  case (divb_powel)
1169  active = .true.
1170  call add_source_powel(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1171  case (divb_janhunen)
1172  active = .true.
1173  call add_source_janhunen(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1174  case (divb_linde)
1175  active = .true.
1176  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1177  case (divb_lindejanhunen)
1178  active = .true.
1179  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1180  call add_source_janhunen(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1181  case (divb_lindepowel)
1182  active = .true.
1183  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1184  call add_source_powel(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1185  case (divb_lindeglm)
1186  active = .true.
1187  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1188  call add_source_glm(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1189  case (divb_ct)
1190  continue ! Do nothing
1191  case (divb_multigrid)
1192  continue ! Do nothing
1193  case default
1194  call mpistop('Unknown divB fix')
1195  end select
1196  else if(source_split_divb .and. qsourcesplit) then
1197  ! Sources related to div B
1198  select case (type_divb)
1199  case (divb_none)
1200  ! Do nothing
1201  case (divb_glm)
1202  active = .true.
1203  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
1204  case (divb_powel)
1205  active = .true.
1206  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
1207  case (divb_janhunen)
1208  active = .true.
1209  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
1210  case (divb_linde)
1211  active = .true.
1212  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1213  case (divb_lindejanhunen)
1214  active = .true.
1215  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1216  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
1217  case (divb_lindepowel)
1218  active = .true.
1219  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1220  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
1221  case (divb_lindeglm)
1222  active = .true.
1223  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1224  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
1225  case (divb_ct)
1226  continue ! Do nothing
1227  case (divb_multigrid)
1228  continue ! Do nothing
1229  case default
1230  call mpistop('Unknown divB fix')
1231  end select
1232  end if
1233  }
1234 
1235  if(mhd_radiative_cooling) then
1236  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,&
1237  w,x,qsourcesplit,active)
1238  end if
1239 
1240  if(mhd_viscosity) then
1241  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
1242  w,x,mhd_energy,qsourcesplit,active)
1243  end if
1244 
1245  if(mhd_gravity) then
1246  call gravity_add_source(qdt,ixi^l,ixo^l,wct,&
1247  w,x,mhd_energy,qsourcesplit,active)
1248  end if
1249 
1250  if (mhd_boris_type == boris_reduced_force) then
1251  call boris_add_source(qdt,ixi^l,ixo^l,wct,&
1252  w,x,qsourcesplit,active)
1253  end if
1254 
1255  end subroutine mhd_add_source
1256 
1257  subroutine boris_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
1258  qsourcesplit,active)
1260  use mod_usr_methods
1261 
1262  integer, intent(in) :: ixI^L, ixO^L
1263  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
1264  double precision, intent(in) :: wCT(ixi^s,1:nw)
1265  double precision, intent(inout) :: w(ixi^s,1:nw)
1266  logical, intent(in) :: qsourcesplit
1267  logical, intent(inout) :: active
1268 
1269  double precision :: JxB(ixi^s,3)
1270  double precision :: gamma_A2(ixo^s)
1271  integer :: idir
1272 
1273  ! Boris source term is always unsplit
1274  if (qsourcesplit) return
1275 
1276  call get_lorentz(ixi^l,ixo^l,w,jxb)
1277  call mhd_gamma2_alfven(ixi^l, ixo^l, wct, gamma_a2)
1278 
1279  do idir = 1, ndir
1280  w(ixo^s,mom(idir)) = w(ixo^s,mom(idir)) + &
1281  qdt * gamma_a2 * jxb(ixo^s, idir)
1282  end do
1283 
1284  end subroutine boris_add_source
1285 
1286  !> Compute the Lorentz force (JxB)
1287  subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
1289  integer, intent(in) :: ixI^L, ixO^L
1290  double precision, intent(in) :: w(ixi^s,1:nw)
1291  double precision, intent(inout) :: JxB(ixi^s,3)
1292  double precision :: a(ixi^s,3), b(ixi^s,3)
1293  integer :: idir, idirmin
1294  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1295  double precision :: current(ixi^s,7-2*ndir:3)
1296 
1297  b=0.0d0
1298  do idir = 1, ndir
1299  b(ixo^s, idir) = mhd_mag_i_all(w, ixi^l, ixo^l,idir)
1300  end do
1301 
1302  ! store J current in a
1303  call get_current(w,ixi^l,ixo^l,idirmin,current)
1304 
1305  a=0.0d0
1306  do idir=7-2*ndir,3
1307  a(ixo^s,idir)=current(ixo^s,idir)
1308  end do
1309 
1310  call cross_product(ixi^l,ixo^l,a,b,jxb)
1311  end subroutine get_lorentz
1312 
1313  !> Compute 1/(1+v_A^2/c^2) for Boris' approximation, where v_A is the Alfven
1314  !> velocity
1315  subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
1317  integer, intent(in) :: ixI^L, ixO^L
1318  double precision, intent(in) :: w(ixi^s, nw)
1319  double precision, intent(out) :: gamma_A2(ixo^s)
1320 
1321  if (mhd_boris_c < 0.0d0) then
1322  ! Good for testing the non-conservative momentum treatment
1323  gamma_a2 = 1.0d0
1324  else
1325  ! Compute the inverse of 1 + B^2/(rho * c^2)
1326  gamma_a2 = 1.0d0 / (1.0d0 + mhd_mag_en_all(w, ixi^l, ixo^l) / &
1327  (w(ixo^s, rho_) * mhd_boris_c**2))
1328  end if
1329  end subroutine mhd_gamma2_alfven
1330 
1331  !> Compute 1/sqrt(1+v_A^2/c^2) for Boris simplification, where v_A is the
1332  !> Alfven velocity
1333  function mhd_gamma_alfven(w, ixI^L, ixO^L) result(gamma_A)
1335  integer, intent(in) :: ixI^L, ixO^L
1336  double precision, intent(in) :: w(ixi^s, nw)
1337  double precision :: gamma_A(ixo^s)
1338 
1339  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a)
1340  gamma_a = sqrt(gamma_a)
1341  end function mhd_gamma_alfven
1342 
1343  subroutine internal_energy_add_source(qdt,ixI^L,ixO^L,wCT,w,x)
1345  use mod_geometry
1346 
1347  integer, intent(in) :: ixI^L, ixO^L
1348  double precision, intent(in) :: qdt
1349  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1350  double precision, intent(inout) :: w(ixi^s,1:nw)
1351  double precision :: pth(ixi^s),v(ixi^s,1:ndir),divv(ixi^s)
1352 
1353  call mhd_get_v(wct,x,ixi^l,ixi^l,v)
1354  call divvector(v,ixi^l,ixo^l,divv)
1355  call mhd_get_pthermal(wct,x,ixi^l,ixo^l,pth)
1356  w(ixo^s,e_)=w(ixo^s,e_)-qdt*pth(ixo^s)*divv(ixo^s)
1357 
1358  end subroutine internal_energy_add_source
1359 
1360  !> Source terms after split off time-independent magnetic field
1361  subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
1363 
1364  integer, intent(in) :: ixI^L, ixO^L
1365  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1366  double precision, intent(inout) :: w(ixi^s,1:nw)
1367 
1368  double precision :: a(ixi^s,3), b(ixi^s,3), axb(ixi^s,3)
1369  integer :: idir
1370 
1371  a=0.d0
1372  b=0.d0
1373  ! for force-free field J0xB0 =0
1374  if(.not.b0field_forcefree) then
1375  ! store B0 magnetic field in b
1376  b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
1377 
1378  ! store J0 current in a
1379  do idir=7-2*ndir,3
1380  a(ixo^s,idir)=block%J0(ixo^s,idir)
1381  end do
1382  call cross_product(ixi^l,ixo^l,a,b,axb)
1383  axb(ixo^s,:)=axb(ixo^s,:)*qdt
1384  ! add J0xB0 source term in momentum equations
1385  w(ixo^s,mom(1:ndir))=w(ixo^s,mom(1:ndir))+axb(ixo^s,1:ndir)
1386  end if
1387 
1388  if(mhd_energy) then
1389  if(.not.block%e_is_internal) then
1390  a=0.d0
1391  ! for free-free field -(vxB0) dot J0 =0
1392  b(ixo^s,:)=wct(ixo^s,mag(:))
1393  ! store full magnetic field B0+B1 in b
1394  if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
1395  ! store velocity in a
1396  do idir=1,ndir
1397  a(ixo^s,idir)=wct(ixo^s,mom(idir))/wct(ixo^s,rho_)
1398  end do
1399  call cross_product(ixi^l,ixo^l,a,b,axb)
1400  axb(ixo^s,:)=axb(ixo^s,:)*qdt
1401  ! add -(vxB) dot J0 source term in energy equation
1402  do idir=7-2*ndir,3
1403  w(ixo^s,e_)=w(ixo^s,e_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
1404  end do
1405  end if
1406  end if
1407 
1408  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
1409 
1410  end subroutine add_source_b0split
1411 
1412  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
1413  !> each direction, non-conservative. If the fourthorder precompiler flag is
1414  !> set, uses fourth order central difference for the laplacian. Then the
1415  !> stencil is 5 (2 neighbours).
1416  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
1418  use mod_usr_methods
1419  use mod_geometry
1420 
1421  integer, intent(in) :: ixI^L, ixO^L
1422  double precision, intent(in) :: qdt
1423  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1424  double precision, intent(inout) :: w(ixi^s,1:nw)
1425  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
1426  integer :: lxO^L, kxO^L
1427 
1428  double precision :: tmp(ixi^s),tmp2(ixi^s)
1429 
1430  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1431  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1432  double precision :: gradeta(ixi^s,1:ndim), Bf(ixi^s,1:ndir)
1433 
1434  ! Calculating resistive sources involve one extra layer
1435  if (mhd_4th_order) then
1436  ixa^l=ixo^l^ladd2;
1437  else
1438  ixa^l=ixo^l^ladd1;
1439  end if
1440 
1441  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1442  call mpistop("Error in add_source_res1: Non-conforming input limits")
1443 
1444  ! Calculate current density and idirmin
1445  call get_current(wct,ixi^l,ixo^l,idirmin,current)
1446 
1447  if (mhd_eta>zero)then
1448  eta(ixa^s)=mhd_eta
1449  gradeta(ixo^s,1:ndim)=zero
1450  else
1451  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
1452  ! assumes that eta is not function of current?
1453  do idim=1,ndim
1454  call gradient(eta,ixi^l,ixo^l,idim,tmp)
1455  gradeta(ixo^s,idim)=tmp(ixo^s)
1456  end do
1457  end if
1458 
1459  if(b0field) then
1460  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
1461  else
1462  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
1463  end if
1464 
1465  do idir=1,ndir
1466  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
1467  if (mhd_4th_order) then
1468  tmp(ixo^s)=zero
1469  tmp2(ixi^s)=bf(ixi^s,idir)
1470  do idim=1,ndim
1471  lxo^l=ixo^l+2*kr(idim,^d);
1472  jxo^l=ixo^l+kr(idim,^d);
1473  hxo^l=ixo^l-kr(idim,^d);
1474  kxo^l=ixo^l-2*kr(idim,^d);
1475  tmp(ixo^s)=tmp(ixo^s)+&
1476  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
1477  /(12.0d0 * dxlevel(idim)**2)
1478  end do
1479  else
1480  tmp(ixo^s)=zero
1481  tmp2(ixi^s)=bf(ixi^s,idir)
1482  do idim=1,ndim
1483  jxo^l=ixo^l+kr(idim,^d);
1484  hxo^l=ixo^l-kr(idim,^d);
1485  tmp(ixo^s)=tmp(ixo^s)+&
1486  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
1487  end do
1488  end if
1489 
1490  ! Multiply by eta
1491  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
1492 
1493  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
1494  if (mhd_eta<zero)then
1495  do jdir=1,ndim; do kdir=idirmin,3
1496  if (lvc(idir,jdir,kdir)/=0)then
1497  if (lvc(idir,jdir,kdir)==1)then
1498  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
1499  else
1500  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
1501  end if
1502  end if
1503  end do; end do
1504  end if
1505 
1506  ! Add sources related to eta*laplB-grad(eta) x J to B and e
1507  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
1508  if (mhd_energy) then
1509  w(ixo^s,e_)=w(ixo^s,e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
1510  end if
1511 
1512  end do ! idir
1513 
1514  if (mhd_energy) then
1515  ! de/dt+=eta*J**2
1516  tmp(ixo^s)=zero
1517  do idir=idirmin,3
1518  tmp(ixo^s)=tmp(ixo^s)+current(ixo^s,idir)**2
1519  end do
1520  w(ixo^s,e_)=w(ixo^s,e_)+qdt*eta(ixo^s)*tmp(ixo^s)
1521  end if
1522 
1523  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
1524 
1525  end subroutine add_source_res1
1526 
1527  !> Add resistive source to w within ixO
1528  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
1529  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
1531  use mod_usr_methods
1532  use mod_geometry
1533 
1534  integer, intent(in) :: ixI^L, ixO^L
1535  double precision, intent(in) :: qdt
1536  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1537  double precision, intent(inout) :: w(ixi^s,1:nw)
1538 
1539  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1540  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
1541  double precision :: tmpvec(ixi^s,1:3)
1542  integer :: ixA^L,idir,idirmin,idirmin1
1543 
1544  ixa^l=ixo^l^ladd2;
1545 
1546  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1547  call mpistop("Error in add_source_res2: Non-conforming input limits")
1548 
1549  ixa^l=ixo^l^ladd1;
1550  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
1551  ! Determine exact value of idirmin while doing the loop.
1552  call get_current(wct,ixi^l,ixa^l,idirmin,current)
1553 
1554  if (mhd_eta>zero)then
1555  eta(ixa^s)=mhd_eta
1556  else
1557  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
1558  end if
1559 
1560  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
1561  tmpvec(ixa^s,1:ndir)=zero
1562  do idir=idirmin,3
1563  tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
1564  end do
1565  curlj=0.d0
1566  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
1567  if(stagger_grid.and.ndim==2.and.ndir==3) then
1568  ! if 2.5D
1569  w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
1570  else
1571  w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
1572  end if
1573 
1574  if(mhd_energy) then
1575  ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
1576  ! de1/dt= eta J^2 - B1 dot curl(eta J)
1577  w(ixo^s,e_)=w(ixo^s,e_)+qdt*(sum(current(ixo^s,:)**2,dim=ndim+1)*eta(ixo^s)-&
1578  sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
1579  end if
1580 
1581  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
1582 
1583  end subroutine add_source_res2
1584 
1585  !> Add Hyper-resistive source to w within ixO
1586  !> Uses 9 point stencil (4 neighbours) in each direction.
1587  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
1589  use mod_geometry
1590 
1591  integer, intent(in) :: ixI^L, ixO^L
1592  double precision, intent(in) :: qdt
1593  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1594  double precision, intent(inout) :: w(ixi^s,1:nw)
1595  !.. local ..
1596  double precision :: current(ixi^s,7-2*ndir:3)
1597  double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
1598  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
1599 
1600  ixa^l=ixo^l^ladd3;
1601  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1602  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
1603 
1604  call get_current(wct,ixi^l,ixa^l,idirmin,current)
1605  tmpvec(ixa^s,1:ndir)=zero
1606  do jdir=idirmin,3
1607  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
1608  end do
1609 
1610  ixa^l=ixo^l^ladd2;
1611  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
1612 
1613  ixa^l=ixo^l^ladd1;
1614  tmpvec(ixa^s,1:ndir)=zero
1615  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
1616  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mhd_eta_hyper
1617 
1618  ixa^l=ixo^l;
1619  tmpvec2(ixa^s,1:ndir)=zero
1620  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
1621 
1622  do idir=1,ndir
1623  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
1624  end do
1625 
1626  if (mhd_energy) then
1627  ! de/dt= +div(B x Ehyper)
1628  ixa^l=ixo^l^ladd1;
1629  tmpvec2(ixa^s,1:ndir)=zero
1630  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
1631  tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
1632  + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
1633  end do; end do; end do
1634  tmp(ixo^s)=zero
1635  call divvector(tmpvec2,ixi^l,ixo^l,tmp)
1636  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
1637  end if
1638 
1639  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
1640 
1641  end subroutine add_source_hyperres
1642 
1643  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
1644  ! Add divB related sources to w within ixO
1645  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
1646  ! giving the EGLM-MHD scheme
1648  use mod_geometry
1649 
1650  integer, intent(in) :: ixI^L, ixO^L
1651  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1652  double precision, intent(inout) :: w(ixi^s,1:nw)
1653  double precision:: divb(ixi^s)
1654  integer :: idim,idir
1655  double precision :: gradPsi(ixi^s)
1656 
1657  ! We calculate now div B
1658  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1659 
1660  ! dPsi/dt = - Ch^2/Cp^2 Psi
1661  if (mhd_glm_alpha < zero) then
1662  w(ixo^s,psi_) = abs(mhd_glm_alpha)*wct(ixo^s,psi_)
1663  else
1664  ! implicit update of Psi variable
1665  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
1666  if(slab_uniform) then
1667  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
1668  else
1669  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
1670  end if
1671  end if
1672 
1673  ! gradient of Psi
1674  do idim=1,ndim
1675  select case(typegrad)
1676  case("central")
1677  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1678  case("limited")
1679  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1680  end select
1681  if (mhd_energy .and. .not.block%e_is_internal) then
1682  ! e = e -qdt (b . grad(Psi))
1683  w(ixo^s,e_) = w(ixo^s,e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
1684  end if
1685  end do
1686 
1687  ! m = m - qdt b div b
1688  do idir=1,ndir
1689  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
1690  end do
1691 
1692  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
1693 
1694  end subroutine add_source_glm
1695 
1696  !> Add divB related sources to w within ixO corresponding to Powel
1697  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
1699 
1700  integer, intent(in) :: ixI^L, ixO^L
1701  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1702  double precision, intent(inout) :: w(ixi^s,1:nw)
1703  double precision :: divb(ixi^s),v(ixi^s,1:ndir)
1704  integer :: idir
1705 
1706  ! We calculate now div B
1707  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1708 
1709  ! calculate velocity
1710  call mhd_get_v(wct,x,ixi^l,ixo^l,v)
1711 
1712  if (mhd_energy .and. .not.block%e_is_internal) then
1713  ! e = e - qdt (v . b) * div b
1714  w(ixo^s,e_)=w(ixo^s,e_)-&
1715  qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
1716  end if
1717 
1718  ! b = b - qdt v * div b
1719  do idir=1,ndir
1720  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
1721  end do
1722 
1723  ! m = m - qdt b div b
1724  do idir=1,ndir
1725  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
1726  end do
1727 
1728  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
1729 
1730  end subroutine add_source_powel
1731 
1732  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
1733  ! Add divB related sources to w within ixO
1734  ! corresponding to Janhunen, just the term in the induction equation.
1736 
1737  integer, intent(in) :: ixI^L, ixO^L
1738  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1739  double precision, intent(inout) :: w(ixi^s,1:nw)
1740  double precision :: divb(ixi^s)
1741  integer :: idir
1742 
1743  ! We calculate now div B
1744  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1745 
1746  ! b = b - qdt v * div b
1747  do idir=1,ndir
1748  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))/wct(ixo^s,rho_)*divb(ixo^s)
1749  end do
1750 
1751  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
1752 
1753  end subroutine add_source_janhunen
1754 
1755  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
1756  ! Add Linde's divB related sources to wnew within ixO
1758  use mod_geometry
1759 
1760  integer, intent(in) :: ixI^L, ixO^L
1761  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1762  double precision, intent(inout) :: w(ixi^s,1:nw)
1763  integer :: idim, idir, ixp^L, i^D, iside
1764  double precision :: divb(ixi^s),graddivb(ixi^s)
1765  logical, dimension(-1:1^D&) :: leveljump
1766 
1767  ! Calculate div B
1768  ixp^l=ixo^l^ladd1;
1769  call get_divb(wct,ixi^l,ixp^l,divb, mhd_divb_4thorder)
1770 
1771  ! for AMR stability, retreat one cell layer from the boarders of level jump
1772  {do i^db=-1,1\}
1773  if(i^d==0|.and.) cycle
1774  if(neighbor_type(i^d,saveigrid)==2 .or. neighbor_type(i^d,saveigrid)==4) then
1775  leveljump(i^d)=.true.
1776  else
1777  leveljump(i^d)=.false.
1778  end if
1779  {end do\}
1780 
1781  ixp^l=ixo^l;
1782  do idim=1,ndim
1783  select case(idim)
1784  {case(^d)
1785  do iside=1,2
1786  i^dd=kr(^dd,^d)*(2*iside-3);
1787  if (leveljump(i^dd)) then
1788  if (iside==1) then
1789  ixpmin^d=ixomin^d-i^d
1790  else
1791  ixpmax^d=ixomax^d-i^d
1792  end if
1793  end if
1794  end do
1795  \}
1796  end select
1797  end do
1798 
1799  ! Add Linde's diffusive terms
1800  do idim=1,ndim
1801  ! Calculate grad_idim(divb)
1802  select case(typegrad)
1803  case("central")
1804  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
1805  case("limited")
1806  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
1807  end select
1808 
1809  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
1810  if (slab_uniform) then
1811  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
1812  else
1813  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
1814  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
1815  end if
1816 
1817  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
1818 
1819  if (mhd_energy .and. typedivbdiff=='all' .and. .not.block%e_is_internal) then
1820  ! e += B_idim*eta*grad_idim(divb)
1821  w(ixp^s,e_)=w(ixp^s,e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
1822  end if
1823  end do
1824 
1825  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
1826 
1827  end subroutine add_source_linde
1828 
1829  !> Calculate div B within ixO
1830  subroutine get_divb(w,ixI^L,ixO^L,divb, fourthorder)
1833  use mod_geometry
1834 
1835  integer, intent(in) :: ixI^L, ixO^L
1836  double precision, intent(in) :: w(ixi^s,1:nw)
1837  double precision, intent(inout) :: divb(ixi^s)
1838  logical, intent(in), optional :: fourthorder
1839 
1840  double precision :: bvec(ixi^s,1:ndir)
1841  double precision :: divb_corner(ixi^s), sign
1842  double precision :: aux_vol(ixi^s)
1843  integer :: ixC^L, idir, ic^D, ix^L
1844 
1845  if(stagger_grid) then
1846  divb=0.d0
1847  do idir=1,ndim
1848  ixc^l=ixo^l-kr(idir,^d);
1849  divb(ixo^s)=divb(ixo^s)+block%ws(ixo^s,idir)*block%surfaceC(ixo^s,idir)-&
1850  block%ws(ixc^s,idir)*block%surfaceC(ixc^s,idir)
1851  end do
1852  divb(ixo^s)=divb(ixo^s)/block%dvolume(ixo^s)
1853  else
1854  bvec(ixi^s,:)=w(ixi^s,mag(:))
1855  select case(typediv)
1856  case("central")
1857  call divvector(bvec,ixi^l,ixo^l,divb,fourthorder)
1858  case("limited")
1859  call divvectors(bvec,ixi^l,ixo^l,divb)
1860  end select
1861  end if
1862 
1863  end subroutine get_divb
1864 
1865  !> get dimensionless div B = |divB| * volume / area / |B|
1866  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
1869 
1870  integer, intent(in) :: ixI^L, ixO^L
1871  double precision, intent(in) :: w(ixi^s,1:nw)
1872  double precision :: divb(ixi^s), dsurface(ixi^s)
1873 
1874  integer :: ixA^L,idims
1875 
1876  call get_divb(w,ixi^l,ixo^l,divb)
1877  if(slab_uniform) then
1878  divb(ixo^s)=0.5d0*abs(divb(ixo^s))/sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))/sum(1.d0/dxlevel(:))
1879  else
1880  ixamin^d=ixomin^d-1;
1881  ixamax^d=ixomax^d-1;
1882  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
1883  do idims=1,ndim
1884  ixa^l=ixo^l-kr(idims,^d);
1885  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
1886  end do
1887  divb(ixo^s)=abs(divb(ixo^s))/sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))*&
1888  block%dvolume(ixo^s)/dsurface(ixo^s)
1889  end if
1890 
1891  end subroutine get_normalized_divb
1892 
1893  !> Calculate idirmin and the idirmin:3 components of the common current array
1894  !> make sure that dxlevel(^D) is set correctly.
1895  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
1897  use mod_geometry
1898 
1899  integer, intent(in) :: ixO^L, ixI^L
1900  double precision, intent(in) :: w(ixi^s,1:nw)
1901  integer, intent(out) :: idirmin
1902  integer :: idir, idirmin0
1903 
1904  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1905  double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)
1906 
1907  idirmin0 = 7-2*ndir
1908 
1909  bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
1910 
1911  call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1912 
1913  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
1914  block%J0(ixo^s,idirmin0:3)
1915 
1916  end subroutine get_current
1917 
1918  !> If resistivity is not zero, check diffusion time limit for dt
1919  subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1921  use mod_usr_methods
1923  use mod_viscosity, only: viscosity_get_dt
1924  use mod_gravity, only: gravity_get_dt
1925 
1926  integer, intent(in) :: ixI^L, ixO^L
1927  double precision, intent(inout) :: dtnew
1928  double precision, intent(in) :: dx^D
1929  double precision, intent(in) :: w(ixi^s,1:nw)
1930  double precision, intent(in) :: x(ixi^s,1:ndim)
1931 
1932  integer :: idirmin,idim
1933  double precision :: dxarr(ndim)
1934  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1935 
1936  dtnew = bigdouble
1937 
1938  ^d&dxarr(^d)=dx^d;
1939  if (mhd_eta>zero)then
1940  dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mhd_eta
1941  else if (mhd_eta<zero)then
1942  call get_current(w,ixi^l,ixo^l,idirmin,current)
1943  call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
1944  dtnew=bigdouble
1945  do idim=1,ndim
1946  if(slab_uniform) then
1947  dtnew=min(dtnew,&
1948  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
1949  else
1950  dtnew=min(dtnew,&
1951  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
1952  end if
1953  end do
1954  end if
1955 
1956  if(mhd_eta_hyper>zero) then
1957  if(slab_uniform) then
1958  dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mhd_eta_hyper,dtnew)
1959  else
1960  dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mhd_eta_hyper,dtnew)
1961  end if
1962  end if
1963 
1964  if(mhd_radiative_cooling) then
1965  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1966  end if
1967 
1968  if(mhd_viscosity) then
1969  call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1970  end if
1971 
1972  if(mhd_gravity) then
1973  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1974  end if
1975 
1976  end subroutine mhd_get_dt
1977 
1978  ! Add geometrical source terms to w
1979  subroutine mhd_add_source_geom(qdt,ixI^L,ixO^L,wCT,w,x)
1981  use mod_geometry
1982 
1983  integer, intent(in) :: ixI^L, ixO^L
1984  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
1985  double precision, intent(inout) :: wCT(ixi^s,1:nw), w(ixi^s,1:nw)
1986 
1987  integer :: iw,idir, h1x^L{^NOONED, h2x^L}
1988  double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
1989 
1990  integer :: mr_,mphi_ ! Polar var. names
1991  integer :: br_,bphi_
1992 
1993  mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1994  br_=mag(1); bphi_=mag(1)-1+phi_
1995 
1996  select case (coordinate)
1997  case (cylindrical)
1998  if (angmomfix) then
1999  call mpistop("angmomfix not implemented yet in MHD")
2000  endif
2001  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp)
2002  if(phi_>0) then
2003  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
2004  wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/wct(ixo^s,rho_))
2005  w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
2006  -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/wct(ixo^s,rho_) &
2007  +wct(ixo^s,bphi_)*wct(ixo^s,br_))
2008  if(.not.stagger_grid) then
2009  w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
2010  (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
2011  -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
2012  /wct(ixo^s,rho_)
2013  end if
2014  else
2015  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
2016  end if
2017  if(mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
2018  case (spherical)
2019  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
2020  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp1)
2021  tmp(ixo^s)=tmp1(ixo^s)
2022  if(b0field) then
2023  tmp2(ixo^s)=sum(block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
2024  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
2025  end if
2026  ! m1
2027  tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
2028  *(block%surfaceC(ixo^s,1)-block%surfaceC(h1x^s,1))/block%dvolume(ixo^s)
2029  if(ndir>1) then
2030  do idir=2,ndir
2031  tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(idir))**2/wct(ixo^s,rho_)-wct(ixo^s,mag(idir))**2
2032  if(b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
2033  end do
2034  end if
2035  w(ixo^s,mom(1))=w(ixo^s,mom(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
2036  ! b1
2037  if(mhd_glm) then
2038  w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,psi_)
2039  end if
2040 
2041  {^nooned
2042  ! m2
2043  tmp(ixo^s)=tmp1(ixo^s)
2044  if(b0field) then
2045  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
2046  end if
2047  ! This will make hydrostatic p=const an exact solution
2048  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s) &
2049  *(block%surfaceC(ixo^s,2)-block%surfaceC(h2x^s,2)) &
2050  /block%dvolume(ixo^s)
2051  tmp(ixo^s)=-(wct(ixo^s,mom(1))*wct(ixo^s,mom(2))/wct(ixo^s,rho_) &
2052  -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
2053  if (b0field) then
2054  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
2055  +wct(ixo^s,mag(1))*block%B0(ixo^s,2,0)
2056  end if
2057  if(ndir==3) then
2058  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(3))**2/wct(ixo^s,rho_) &
2059  -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
2060  if (b0field) then
2061  tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
2062  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
2063  end if
2064  end if
2065  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
2066  ! b2
2067  if(.not.stagger_grid) then
2068  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
2069  -wct(ixo^s,mom(2))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_)
2070  if(b0field) then
2071  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,2,0) &
2072  -wct(ixo^s,mom(2))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_)
2073  end if
2074  if(mhd_glm) then
2075  tmp(ixo^s)=tmp(ixo^s) &
2076  + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,psi_)
2077  end if
2078  w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
2079  end if
2080  }
2081 
2082  if(ndir==3) then
2083  ! m3
2084  if(.not.angmomfix) then
2085  tmp(ixo^s)=-(wct(ixo^s,mom(3))*wct(ixo^s,mom(1))/wct(ixo^s,rho_) &
2086  -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
2087  -(wct(ixo^s,mom(2))*wct(ixo^s,mom(3))/wct(ixo^s,rho_) &
2088  -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
2089  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
2090  if (b0field) then
2091  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
2092  +wct(ixo^s,mag(1))*block%B0(ixo^s,3,0) {^nooned &
2093  +(block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
2094  +wct(ixo^s,mag(2))*block%B0(ixo^s,3,0)) &
2095  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
2096  end if
2097  w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
2098  else
2099  call mpistop("angmomfix not implemented yet in MHD")
2100  end if
2101  ! b3
2102  if(.not.stagger_grid) then
2103  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
2104  -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_) {^nooned &
2105  -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
2106  -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
2107  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
2108  if (b0field) then
2109  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,3,0) &
2110  -wct(ixo^s,mom(3))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_){^nooned &
2111  -(wct(ixo^s,mom(3))*block%B0(ixo^s,2,0) &
2112  -wct(ixo^s,mom(2))*block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
2113  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
2114  end if
2115  w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
2116  end if
2117  end if
2118  end select
2119  end subroutine mhd_add_source_geom
2120 
2121  !> Compute 2 times total magnetic energy
2122  function mhd_mag_en_all(w, ixI^L, ixO^L) result(mge)
2124  integer, intent(in) :: ixI^L, ixO^L
2125  double precision, intent(in) :: w(ixi^s, nw)
2126  double precision :: mge(ixo^s)
2127 
2128  if (b0field) then
2129  mge = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,block%iw0))**2, dim=ndim+1)
2130  else
2131  mge = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
2132  end if
2133  end function mhd_mag_en_all
2134 
2135  !> Compute full magnetic field by direction
2136  function mhd_mag_i_all(w, ixI^L, ixO^L,idir) result(mgf)
2138  integer, intent(in) :: ixI^L, ixO^L, idir
2139  double precision, intent(in) :: w(ixi^s, nw)
2140  double precision :: mgf(ixo^s)
2141 
2142  if (b0field) then
2143  mgf = w(ixo^s, mag(idir))+block%B0(ixo^s,idir,block%iw0)
2144  else
2145  mgf = w(ixo^s, mag(idir))
2146  end if
2147  end function mhd_mag_i_all
2148 
2149  !> Compute evolving magnetic energy
2150  function mhd_mag_en(w, ixI^L, ixO^L) result(mge)
2152  integer, intent(in) :: ixI^L, ixO^L
2153  double precision, intent(in) :: w(ixi^s, nw)
2154  double precision :: mge(ixo^s)
2155 
2156  mge = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
2157  end function mhd_mag_en
2158 
2159  !> compute kinetic energy
2160  function mhd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
2162  integer, intent(in) :: ixI^L, ixO^L
2163  double precision, intent(in) :: w(ixi^s, nw)
2164  double precision :: ke(ixo^s)
2165  double precision, intent(in), optional :: inv_rho(ixo^s)
2166 
2167  if (present(inv_rho)) then
2168  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
2169  else
2170  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
2171  end if
2172  end function mhd_kin_en
2173 
2174  subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
2176 
2177  integer, intent(in) :: ixI^L, ixO^L
2178  double precision, intent(in) :: w(ixi^s,nw)
2179  double precision, intent(in) :: x(ixi^s,1:ndim)
2180  double precision, intent(inout) :: vHall(ixi^s,1:3)
2181 
2182  integer :: idir, idirmin
2183  double precision :: current(ixi^s,7-2*ndir:3)
2184 
2185  ! Calculate current density and idirmin
2186  call get_current(w,ixi^l,ixo^l,idirmin,current)
2187  vhall(ixo^s,1:3) = zero
2188  vhall(ixo^s,idirmin:3) = - mhd_etah*current(ixo^s,idirmin:3)
2189  do idir = idirmin, 3
2190  vhall(ixo^s,idir) = vhall(ixo^s,idir)/w(ixo^s,rho_)
2191  end do
2192 
2193  end subroutine mhd_getv_hall
2194 
2195  subroutine mhd_getdt_hall(w,x,ixI^L,ixO^L,dx^D,dthall)
2197 
2198  integer, intent(in) :: ixI^L, ixO^L
2199  double precision, intent(in) :: dx^D
2200  double precision, intent(in) :: w(ixi^s,1:nw)
2201  double precision, intent(in) :: x(ixi^s,1:ndim)
2202  double precision, intent(out) :: dthall
2203  !.. local ..
2204  double precision :: dxarr(ndim)
2205  double precision :: bmag(ixi^s)
2206 
2207  dthall=bigdouble
2208 
2209  ! because we have that in cmax now:
2210  return
2211 
2212  ^d&dxarr(^d)=dx^d;
2213 
2214  if (.not. b0field) then
2215  bmag(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2, dim=ndim+1))
2216  bmag(ixo^s)=sqrt(sum((w(ixo^s,mag(:)) + block%B0(ixo^s,1:ndir,block%iw0))**2))
2217  end if
2218 
2219  if(slab_uniform) then
2220  dthall=dtdiffpar*minval(dxarr(1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
2221  else
2222  dthall=dtdiffpar*minval(block%ds(ixo^s,1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
2223  end if
2224 
2225  end subroutine mhd_getdt_hall
2226 
2227  subroutine mhd_modify_wlr(ixI^L,ixO^L,wLC,wRC,wLp,wRp,s,idir)
2229  integer, intent(in) :: ixI^L, ixO^L, idir
2230  double precision, intent(inout) :: wLC(ixi^s,1:nw), wRC(ixi^s,1:nw)
2231  double precision, intent(inout) :: wLp(ixi^s,1:nw), wRp(ixi^s,1:nw)
2232  type(state) :: s
2233  double precision :: dB(ixi^s), dPsi(ixi^s)
2234 
2235  if(stagger_grid) then
2236  wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2237  wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2238  wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2239  wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2240  else
2241  ! Solve the Riemann problem for the linear 2x2 system for normal
2242  ! B-field and GLM_Psi according to Dedner 2002:
2243  ! This implements eq. (42) in Dedner et al. 2002 JcP 175
2244  ! Gives the Riemann solution on the interface
2245  ! for the normal B component and Psi in the GLM-MHD system.
2246  ! 23/04/2013 Oliver Porth
2247  db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
2248  dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
2249 
2250  wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
2251  - 0.5d0/cmax_global * dpsi(ixo^s)
2252  wlp(ixo^s,psi_) = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
2253  - 0.5d0*cmax_global * db(ixo^s)
2254 
2255  wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
2256  wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
2257  end if
2258 
2259  end subroutine mhd_modify_wlr
2260 
2261  subroutine mhd_boundary_adjust
2263  integer :: iB, idim, iside, iigrid, igrid
2264  integer :: ixG^L, ixO^L, i^D
2265 
2266  ixg^l=ixg^ll;
2267  do iigrid=1,igridstail; igrid=igrids(iigrid);
2268  if(.not.phyboundblock(igrid)) cycle
2269  saveigrid=igrid
2270  block=>ps(igrid)
2271  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2272  do idim=1,ndim
2273  ! to avoid using as yet unknown corner info in more than 1D, we
2274  ! fill only interior mesh ranges of the ghost cell ranges at first,
2275  ! and progressively enlarge the ranges to include corners later
2276  do iside=1,2
2277  i^d=kr(^d,idim)*(2*iside-3);
2278  if (neighbor_type(i^d,igrid)/=1) cycle
2279  ib=(idim-1)*2+iside
2280  if(.not.boundary_divbfix(ib)) cycle
2281  if(any(typeboundary(:,ib)=="special")) then
2282  ! MF nonlinear force-free B field extrapolation and data driven
2283  ! require normal B of the first ghost cell layer to be untouched by
2284  ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
2285  select case (idim)
2286  {case (^d)
2287  if (iside==2) then
2288  ! maximal boundary
2289  ixomin^dd=ixgmax^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixgmin^dd;
2290  ixomax^dd=ixgmax^dd;
2291  else
2292  ! minimal boundary
2293  ixomin^dd=ixgmin^dd;
2294  ixomax^dd=ixgmin^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixgmax^dd;
2295  end if \}
2296  end select
2297  call fixdivb_boundary(ixg^l,ixo^l,ps(igrid)%w,ps(igrid)%x,ib)
2298  end if
2299  end do
2300  end do
2301  end do
2302 
2303  end subroutine mhd_boundary_adjust
2304 
2305  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
2307 
2308  integer, intent(in) :: ixG^L,ixO^L,iB
2309  double precision, intent(inout) :: w(ixg^s,1:nw)
2310  double precision, intent(in) :: x(ixg^s,1:ndim)
2311 
2312  double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
2313  integer :: ix^D,ixF^L
2314 
2315  select case(ib)
2316  case(1)
2317  ! 2nd order CD for divB=0 to set normal B component better
2318  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2319  {^iftwod
2320  ixfmin1=ixomin1+1
2321  ixfmax1=ixomax1+1
2322  ixfmin2=ixomin2+1
2323  ixfmax2=ixomax2-1
2324  if(slab_uniform) then
2325  dx1x2=dxlevel(1)/dxlevel(2)
2326  do ix1=ixfmax1,ixfmin1,-1
2327  w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
2328  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
2329  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
2330  enddo
2331  else
2332  do ix1=ixfmax1,ixfmin1,-1
2333  w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
2334  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
2335  +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
2336  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
2337  -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
2338  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
2339  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
2340  end do
2341  end if
2342  }
2343  {^ifthreed
2344  ixfmin1=ixomin1+1
2345  ixfmax1=ixomax1+1
2346  ixfmin2=ixomin2+1
2347  ixfmax2=ixomax2-1
2348  ixfmin3=ixomin3+1
2349  ixfmax3=ixomax3-1
2350  if(slab_uniform) then
2351  dx1x2=dxlevel(1)/dxlevel(2)
2352  dx1x3=dxlevel(1)/dxlevel(3)
2353  do ix1=ixfmax1,ixfmin1,-1
2354  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2355  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
2356  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
2357  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
2358  +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
2359  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
2360  end do
2361  else
2362  do ix1=ixfmax1,ixfmin1,-1
2363  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2364  ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
2365  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
2366  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
2367  +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
2368  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
2369  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
2370  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
2371  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
2372  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
2373  +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
2374  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
2375  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
2376  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
2377  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2378  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
2379  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
2380  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
2381  end do
2382  end if
2383  }
2384  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2385  case(2)
2386  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2387  {^iftwod
2388  ixfmin1=ixomin1-1
2389  ixfmax1=ixomax1-1
2390  ixfmin2=ixomin2+1
2391  ixfmax2=ixomax2-1
2392  if(slab_uniform) then
2393  dx1x2=dxlevel(1)/dxlevel(2)
2394  do ix1=ixfmin1,ixfmax1
2395  w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
2396  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
2397  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
2398  enddo
2399  else
2400  do ix1=ixfmin1,ixfmax1
2401  w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
2402  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
2403  -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
2404  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
2405  +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
2406  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
2407  /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
2408  end do
2409  end if
2410  }
2411  {^ifthreed
2412  ixfmin1=ixomin1-1
2413  ixfmax1=ixomax1-1
2414  ixfmin2=ixomin2+1
2415  ixfmax2=ixomax2-1
2416  ixfmin3=ixomin3+1
2417  ixfmax3=ixomax3-1
2418  if(slab_uniform) then
2419  dx1x2=dxlevel(1)/dxlevel(2)
2420  dx1x3=dxlevel(1)/dxlevel(3)
2421  do ix1=ixfmin1,ixfmax1
2422  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2423  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
2424  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
2425  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
2426  -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
2427  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
2428  end do
2429  else
2430  do ix1=ixfmin1,ixfmax1
2431  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2432  ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
2433  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
2434  block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
2435  -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
2436  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
2437  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
2438  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
2439  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
2440  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
2441  -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
2442  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
2443  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
2444  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
2445  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2446  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
2447  /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
2448  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
2449  end do
2450  end if
2451  }
2452  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2453  case(3)
2454  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2455  {^iftwod
2456  ixfmin1=ixomin1+1
2457  ixfmax1=ixomax1-1
2458  ixfmin2=ixomin2+1
2459  ixfmax2=ixomax2+1
2460  if(slab_uniform) then
2461  dx2x1=dxlevel(2)/dxlevel(1)
2462  do ix2=ixfmax2,ixfmin2,-1
2463  w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
2464  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
2465  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
2466  enddo
2467  else
2468  do ix2=ixfmax2,ixfmin2,-1
2469  w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
2470  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
2471  +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
2472  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
2473  -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
2474  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
2475  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
2476  end do
2477  end if
2478  }
2479  {^ifthreed
2480  ixfmin1=ixomin1+1
2481  ixfmax1=ixomax1-1
2482  ixfmin3=ixomin3+1
2483  ixfmax3=ixomax3-1
2484  ixfmin2=ixomin2+1
2485  ixfmax2=ixomax2+1
2486  if(slab_uniform) then
2487  dx2x1=dxlevel(2)/dxlevel(1)
2488  dx2x3=dxlevel(2)/dxlevel(3)
2489  do ix2=ixfmax2,ixfmin2,-1
2490  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
2491  ix2+1,ixfmin3:ixfmax3,mag(2)) &
2492  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
2493  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
2494  +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
2495  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
2496  end do
2497  else
2498  do ix2=ixfmax2,ixfmin2,-1
2499  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
2500  ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
2501  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
2502  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
2503  +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
2504  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2505  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
2506  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
2507  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2508  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
2509  +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
2510  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
2511  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
2512  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
2513  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2514  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
2515  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
2516  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
2517  end do
2518  end if
2519  }
2520  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2521  case(4)
2522  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2523  {^iftwod
2524  ixfmin1=ixomin1+1
2525  ixfmax1=ixomax1-1
2526  ixfmin2=ixomin2-1
2527  ixfmax2=ixomax2-1
2528  if(slab_uniform) then
2529  dx2x1=dxlevel(2)/dxlevel(1)
2530  do ix2=ixfmin2,ixfmax2
2531  w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
2532  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
2533  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
2534  end do
2535  else
2536  do ix2=ixfmin2,ixfmax2
2537  w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
2538  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
2539  -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
2540  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
2541  +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
2542  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
2543  /block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
2544  end do
2545  end if
2546  }
2547  {^ifthreed
2548  ixfmin1=ixomin1+1
2549  ixfmax1=ixomax1-1
2550  ixfmin3=ixomin3+1
2551  ixfmax3=ixomax3-1
2552  ixfmin2=ixomin2-1
2553  ixfmax2=ixomax2-1
2554  if(slab_uniform) then
2555  dx2x1=dxlevel(2)/dxlevel(1)
2556  dx2x3=dxlevel(2)/dxlevel(3)
2557  do ix2=ixfmin2,ixfmax2
2558  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
2559  ix2-1,ixfmin3:ixfmax3,mag(2)) &
2560  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
2561  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
2562  -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
2563  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
2564  end do
2565  else
2566  do ix2=ixfmin2,ixfmax2
2567  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
2568  ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
2569  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
2570  block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
2571  -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
2572  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2573  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
2574  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
2575  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2576  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
2577  -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
2578  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
2579  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
2580  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
2581  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2582  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
2583  /block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
2584  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
2585  end do
2586  end if
2587  }
2588  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2589  {^ifthreed
2590  case(5)
2591  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2592  ixfmin1=ixomin1+1
2593  ixfmax1=ixomax1-1
2594  ixfmin2=ixomin2+1
2595  ixfmax2=ixomax2-1
2596  ixfmin3=ixomin3+1
2597  ixfmax3=ixomax3+1
2598  if(slab_uniform) then
2599  dx3x1=dxlevel(3)/dxlevel(1)
2600  dx3x2=dxlevel(3)/dxlevel(2)
2601  do ix3=ixfmax3,ixfmin3,-1
2602  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
2603  ixfmin2:ixfmax2,ix3+1,mag(3)) &
2604  +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
2605  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
2606  +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
2607  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
2608  end do
2609  else
2610  do ix3=ixfmax3,ixfmin3,-1
2611  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
2612  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
2613  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
2614  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
2615  +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
2616  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2617  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
2618  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
2619  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2620  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
2621  +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
2622  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
2623  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
2624  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
2625  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
2626  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
2627  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
2628  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
2629  end do
2630  end if
2631  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2632  case(6)
2633  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2634  ixfmin1=ixomin1+1
2635  ixfmax1=ixomax1-1
2636  ixfmin2=ixomin2+1
2637  ixfmax2=ixomax2-1
2638  ixfmin3=ixomin3-1
2639  ixfmax3=ixomax3-1
2640  if(slab_uniform) then
2641  dx3x1=dxlevel(3)/dxlevel(1)
2642  dx3x2=dxlevel(3)/dxlevel(2)
2643  do ix3=ixfmin3,ixfmax3
2644  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
2645  ixfmin2:ixfmax2,ix3-1,mag(3)) &
2646  -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
2647  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
2648  -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
2649  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
2650  end do
2651  else
2652  do ix3=ixfmin3,ixfmax3
2653  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
2654  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
2655  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
2656  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
2657  -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
2658  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2659  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
2660  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
2661  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2662  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
2663  -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
2664  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
2665  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
2666  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
2667  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
2668  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
2669  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
2670  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
2671  end do
2672  end if
2673  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2674  }
2675  case default
2676  call mpistop("Special boundary is not defined for this region")
2677  end select
2678 
2679  end subroutine fixdivb_boundary
2680 
2681  {^nooned
2682  subroutine mhd_clean_divb_multigrid(qdt, qt, active)
2686  use mod_geometry
2687 
2688  double precision, intent(in) :: qdt !< Current time step
2689  double precision, intent(in) :: qt !< Current time
2690  logical, intent(inout) :: active !< Output if the source is active
2691  integer :: iigrid, igrid, id
2692  integer :: n, nc, lvl, ix^L, ixC^L, idim
2693  type(tree_node), pointer :: pnode
2694  double precision :: tmp(ixg^t), grad(ixg^t, ndim)
2695  double precision :: res
2696  double precision, parameter :: max_residual = 1d-3
2697  double precision, parameter :: residual_reduction = 1d-10
2698  integer, parameter :: max_its = 50
2699  double precision :: residual_it(max_its), max_divb
2700 
2701  mg%operator_type = mg_laplacian
2702 
2703  ! Set boundary conditions
2704  do n = 1, 2*ndim
2705  idim = (n+1)/2
2706  select case (typeboundary(mag(idim), n))
2707  case ('symm')
2708  ! d/dx B = 0, take phi = 0
2709  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
2710  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2711  case ('asymm')
2712  ! B = 0, so grad(phi) = 0
2713  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
2714  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2715  case ('cont')
2716  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
2717  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2718  case ('special')
2719  ! Assume Dirichlet boundary conditions, derivative zero
2720  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
2721  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2722  case ('periodic')
2723  ! Nothing to do here
2724  case default
2725  print *, "divb_multigrid warning: unknown b.c.: ", &
2726  trim(typeboundary(mag(idim), n))
2727  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
2728  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2729  end select
2730  end do
2731 
2732  ix^l=ixm^ll^ladd1;
2733  max_divb = 0.0d0
2734 
2735  ! Store divergence of B as right-hand side
2736  do iigrid = 1, igridstail
2737  igrid = igrids(iigrid);
2738  pnode => igrid_to_node(igrid, mype)%node
2739  id = pnode%id
2740  lvl = mg%boxes(id)%lvl
2741  nc = mg%box_size_lvl(lvl)
2742 
2743  ! Geometry subroutines expect this to be set
2744  block => ps(igrid)
2745  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2746 
2747  call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^ll, ixm^ll, tmp, &
2749  mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(ixm^t)
2750  max_divb = max(max_divb, maxval(abs(tmp(ixm^t))))
2751  end do
2752 
2753  ! Solve laplacian(phi) = divB
2754  if(stagger_grid) then
2755  call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
2756  mpi_max, icomm, ierrmpi)
2757 
2758  if (mype == 0) print *, "Performing multigrid divB cleaning"
2759  if (mype == 0) print *, "iteration vs residual"
2760  ! Solve laplacian(phi) = divB
2761  do n = 1, max_its
2762  call mg_fas_fmg(mg, n>1, max_res=residual_it(n))
2763  if (mype == 0) write(*, "(I4,E11.3)") n, residual_it(n)
2764  if (residual_it(n) < residual_reduction * max_divb) exit
2765  end do
2766  if (mype == 0 .and. n > max_its) then
2767  print *, "divb_multigrid warning: not fully converged"
2768  print *, "current amplitude of divb: ", residual_it(max_its)
2769  print *, "multigrid smallest grid: ", &
2770  mg%domain_size_lvl(:, mg%lowest_lvl)
2771  print *, "note: smallest grid ideally has <= 8 cells"
2772  print *, "multigrid dx/dy/dz ratio: ", mg%dr(:, 1)/mg%dr(1, 1)
2773  print *, "note: dx/dy/dz should be similar"
2774  end if
2775  else
2776  do n = 1, max_its
2777  call mg_fas_vcycle(mg, max_res=res)
2778  if (res < max_residual) exit
2779  end do
2780  if (res > max_residual) call mpistop("divb_multigrid: no convergence")
2781  end if
2782 
2783 
2784  ! Correct the magnetic field
2785  do iigrid = 1, igridstail
2786  igrid = igrids(iigrid);
2787  pnode => igrid_to_node(igrid, mype)%node
2788  id = pnode%id
2789 
2790  ! Geometry subroutines expect this to be set
2791  block => ps(igrid)
2792  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2793 
2794  ! Compute the gradient of phi
2795  tmp(ix^s) = mg%boxes(id)%cc({:,}, mg_iphi)
2796 
2797  if(stagger_grid) then
2798  do idim =1, ndim
2799  ixcmin^d=ixmlo^d-kr(idim,^d);
2800  ixcmax^d=ixmhi^d;
2801  call gradientx(tmp,ps(igrid)%x,ixg^ll,ixc^l,idim,grad(ixg^t,idim),.false.)
2802  ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
2803  end do
2804  call mhd_face_to_center(ixm^ll,ps(igrid))
2805  else
2806  do idim = 1, ndim
2807  call gradient(tmp,ixg^ll,ixm^ll,idim,grad(ixg^t, idim))
2808  end do
2809  ! Apply the correction B* = B - gradient(phi)
2810  tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
2811  ps(igrid)%w(ixm^t, mag(1:ndim)) = &
2812  ps(igrid)%w(ixm^t, mag(1:ndim)) - grad(ixm^t, :)
2813  end if
2814 
2815  ! Determine magnetic energy difference
2816  tmp(ixm^t) = 0.5_dp * (sum(ps(igrid)%w(ixm^t, &
2817  mag(1:ndim))**2, dim=ndim+1) - tmp(ixm^t))
2818 
2819  if (mhd_energy) then
2820  ! Keep thermal pressure the same
2821  ps(igrid)%w(ixm^t, e_) = ps(igrid)%w(ixm^t, e_) + tmp(ixm^t)
2822 
2823  ! Possible alternative: keep total pressure the same
2824  ! ps(igrid)%w(ixM^T, e_) = ps(igrid)%w(ixM^T, e_) + &
2825  ! (mhd_gamma-2) * inv_gamma_1 * tmp(ixM^T)
2826  end if
2827  end do
2828 
2829  active = .true.
2830 
2831  end subroutine mhd_clean_divb_multigrid
2832  }
2833 
2834  subroutine mhd_update_faces(ixI^L,ixO^L,qdt,wprim,fC,fE,sCT,s)
2836 
2837  integer, intent(in) :: ixI^L, ixO^L
2838  double precision, intent(in) :: qdt
2839  ! cell-center primitive variables
2840  double precision, intent(in) :: wprim(ixi^s,1:nw)
2841  type(state) :: sCT, s
2842  double precision, intent(in) :: fC(ixi^s,1:nwflux,1:ndim)
2843  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
2844 
2845  select case(type_ct)
2846  case('average')
2847  call update_faces_average(ixi^l,ixo^l,qdt,fc,fe,sct,s)
2848  case('uct_contact')
2849  call update_faces_contact(ixi^l,ixo^l,qdt,wprim,fc,fe,sct,s)
2850  case('uct_hll')
2851  call update_faces_hll(ixi^l,ixo^l,qdt,fe,sct,s)
2852  case default
2853  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
2854  end select
2855 
2856  end subroutine mhd_update_faces
2857 
2858  !> get electric field though averaging neighors to update faces in CT
2859  subroutine update_faces_average(ixI^L,ixO^L,qdt,fC,fE,sCT,s)
2862  use mod_usr_methods
2863 
2864  integer, intent(in) :: ixI^L, ixO^L
2865  double precision, intent(in) :: qdt
2866  type(state) :: sCT, s
2867  double precision, intent(in) :: fC(ixi^s,1:nwflux,1:ndim)
2868  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
2869 
2870  integer :: hxC^L,ixC^L,jxC^L,ixCm^L
2871  integer :: idim1,idim2,idir,iwdim1,iwdim2
2872  double precision :: circ(ixi^s,1:ndim)
2873  ! current on cell edges
2874  double precision :: jce(ixi^s,7-2*ndim:3)
2875 
2876  associate(bfaces=>s%ws,x=>s%x)
2877 
2878  ! Calculate contribution to FEM of each edge,
2879  ! that is, estimate value of line integral of
2880  ! electric field in the positive idir direction.
2881  ixcmax^d=ixomax^d;
2882  ixcmin^d=ixomin^d-1;
2883 
2884  ! if there is resistivity, get eta J
2885  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
2886 
2887  fe=zero
2888 
2889  do idim1=1,ndim
2890  iwdim1 = mag(idim1)
2891  do idim2=1,ndim
2892  iwdim2 = mag(idim2)
2893  do idir=7-2*ndim,3! Direction of line integral
2894  ! Allow only even permutations
2895  if (lvc(idim1,idim2,idir)==1) then
2896  ! Assemble indices
2897  jxc^l=ixc^l+kr(idim1,^d);
2898  hxc^l=ixc^l+kr(idim2,^d);
2899  ! Interpolate to edges
2900  fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
2901  -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
2902 
2903  ! add current component of electric field at cell edges E=-vxB+eta J
2904  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2905 
2906  fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
2907 
2908  if (.not.slab) then
2909  where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
2910  fe(ixc^s,idir)=zero
2911  end where
2912  end if
2913  end if
2914  end do
2915  end do
2916  end do
2917 
2918  ! allow user to change inductive electric field, especially for boundary driven applications
2919  if(associated(usr_set_electric_field)) &
2920  call usr_set_electric_field(ixi^l,ixo^l,qdt,fe,sct)
2921 
2922  circ(ixi^s,1:ndim)=zero
2923 
2924  ! Calculate circulation on each face
2925 
2926  do idim1=1,ndim ! Coordinate perpendicular to face
2927  do idim2=1,ndim
2928  do idir=7-2*ndim,3 ! Direction of line integral
2929  ! Assemble indices
2930  hxc^l=ixc^l-kr(idim2,^d);
2931  ! Add line integrals in direction idir
2932  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2933  +lvc(idim1,idim2,idir)&
2934  *(fe(ixc^s,idir)&
2935  -fe(hxc^s,idir))
2936  end do
2937  end do
2938  end do
2939 
2940  ! Divide by the area of the face to get dB/dt
2941  do idim1=1,ndim
2942  ixcmax^d=ixomax^d;
2943  ixcmin^d=ixomin^d-kr(idim1,^d);
2944  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
2945  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2946  elsewhere
2947  circ(ixc^s,idim1)=zero
2948  end where
2949  ! Time update
2950  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2951  end do
2952 
2953  end associate
2954 
2955  end subroutine update_faces_average
2956 
2957  !> update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
2958  subroutine update_faces_contact(ixI^L,ixO^L,qdt,wp,fC,fE,sCT,s)
2961  use mod_usr_methods
2962 
2963  integer, intent(in) :: ixI^L, ixO^L
2964  double precision, intent(in) :: qdt
2965  ! cell-center primitive variables
2966  double precision, intent(in) :: wp(ixi^s,1:nw)
2967  type(state) :: sCT, s
2968  double precision, intent(in) :: fC(ixi^s,1:nwflux,1:ndim)
2969  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
2970 
2971  double precision :: circ(ixi^s,1:ndim)
2972  ! electric field at cell centers
2973  double precision :: ECC(ixi^s,7-2*ndim:3)
2974  ! gradient of E at left and right side of a cell face
2975  double precision :: EL(ixi^s),ER(ixi^s)
2976  ! gradient of E at left and right side of a cell corner
2977  double precision :: ELC(ixi^s),ERC(ixi^s)
2978  ! current on cell edges
2979  double precision :: jce(ixi^s,7-2*ndim:3)
2980  integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
2981  integer :: idim1,idim2,idir,iwdim1,iwdim2
2982 
2983  associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
2984 
2985  ecc=0.d0
2986  ! Calculate electric field at cell centers
2987  do idim1=1,ndim; do idim2=1,ndim; do idir=7-2*ndim,3
2988  if(lvc(idim1,idim2,idir)==1)then
2989  ecc(ixi^s,idir)=ecc(ixi^s,idir)+wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
2990  else if(lvc(idim1,idim2,idir)==-1) then
2991  ecc(ixi^s,idir)=ecc(ixi^s,idir)-wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
2992  endif
2993  enddo; enddo; enddo
2994 
2995  ! if there is resistivity, get eta J
2996  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
2997 
2998  ! Calculate contribution to FEM of each edge,
2999  ! that is, estimate value of line integral of
3000  ! electric field in the positive idir direction.
3001  fe=zero
3002  ! evaluate electric field along cell edges according to equation (41)
3003  do idim1=1,ndim
3004  iwdim1 = mag(idim1)
3005  do idim2=1,ndim
3006  iwdim2 = mag(idim2)
3007  do idir=7-2*ndim,3 ! Direction of line integral
3008  ! Allow only even permutations
3009  if (lvc(idim1,idim2,idir)==1) then
3010  ixcmax^d=ixomax^d;
3011  ixcmin^d=ixomin^d+kr(idir,^d)-1;
3012  ! Assemble indices
3013  jxc^l=ixc^l+kr(idim1,^d);
3014  hxc^l=ixc^l+kr(idim2,^d);
3015  ! average cell-face electric field to cell edges
3016  fe(ixc^s,idir)=quarter*&
3017  (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
3018  -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
3019 
3020  ! add slope in idim2 direction from equation (50)
3021  ixamin^d=ixcmin^d;
3022  ixamax^d=ixcmax^d+kr(idim1,^d);
3023  el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
3024  hxc^l=ixa^l+kr(idim2,^d);
3025  er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
3026  where(vnorm(ixc^s,idim1)>0.d0)
3027  elc(ixc^s)=el(ixc^s)
3028  else where(vnorm(ixc^s,idim1)<0.d0)
3029  elc(ixc^s)=el(jxc^s)
3030  else where
3031  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
3032  end where
3033  hxc^l=ixc^l+kr(idim2,^d);
3034  where(vnorm(hxc^s,idim1)>0.d0)
3035  erc(ixc^s)=er(ixc^s)
3036  else where(vnorm(hxc^s,idim1)<0.d0)
3037  erc(ixc^s)=er(jxc^s)
3038  else where
3039  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
3040  end where
3041  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
3042 
3043  ! add slope in idim1 direction from equation (50)
3044  jxc^l=ixc^l+kr(idim2,^d);
3045  ixamin^d=ixcmin^d;
3046  ixamax^d=ixcmax^d+kr(idim2,^d);
3047  el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
3048  hxc^l=ixa^l+kr(idim1,^d);
3049  er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
3050  where(vnorm(ixc^s,idim2)>0.d0)
3051  elc(ixc^s)=el(ixc^s)
3052  else where(vnorm(ixc^s,idim2)<0.d0)
3053  elc(ixc^s)=el(jxc^s)
3054  else where
3055  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
3056  end where
3057  hxc^l=ixc^l+kr(idim1,^d);
3058  where(vnorm(hxc^s,idim2)>0.d0)
3059  erc(ixc^s)=er(ixc^s)
3060  else where(vnorm(hxc^s,idim2)<0.d0)
3061  erc(ixc^s)=er(jxc^s)
3062  else where
3063  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
3064  end where
3065  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
3066 
3067  ! add current component of electric field at cell edges E=-vxB+eta J
3068  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
3069 
3070  ! times time step and edge length
3071  fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
3072  if (.not.slab) then
3073  where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
3074  fe(ixc^s,idir)=zero
3075  end where
3076  end if
3077  end if
3078  end do
3079  end do
3080  end do
3081 
3082  ! allow user to change inductive electric field, especially for boundary driven applications
3083  if(associated(usr_set_electric_field)) &
3084  call usr_set_electric_field(ixi^l,ixo^l,qdt,fe,sct)
3085 
3086  circ(ixi^s,1:ndim)=zero
3087 
3088  ! Calculate circulation on each face
3089  do idim1=1,ndim ! Coordinate perpendicular to face
3090  ixcmax^d=ixomax^d;
3091  ixcmin^d=ixomin^d-kr(idim1,^d);
3092  do idim2=1,ndim
3093  do idir=7-2*ndim,3 ! Direction of line integral
3094  ! Assemble indices
3095  hxc^l=ixc^l-kr(idim2,^d);
3096  ! Add line integrals in direction idir
3097  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
3098  +lvc(idim1,idim2,idir)&
3099  *(fe(ixc^s,idir)&
3100  -fe(hxc^s,idir))
3101  end do
3102  end do
3103  ! Divide by the area of the face to get dB/dt
3104  ixcmax^d=ixomax^d;
3105  ixcmin^d=ixomin^d-kr(idim1,^d);
3106  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
3107  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
3108  elsewhere
3109  circ(ixc^s,idim1)=zero
3110  end where
3111  ! Time update cell-face magnetic field component
3112  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
3113  end do
3114 
3115  end associate
3116 
3117  end subroutine update_faces_contact
3118 
3119  !> update faces
3120  subroutine update_faces_hll(ixI^L,ixO^L,qdt,fE,sCT,s)
3123  use mod_usr_methods
3124 
3125  integer, intent(in) :: ixI^L, ixO^L
3126  double precision, intent(in) :: qdt
3127  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
3128  type(state) :: sCT, s
3129 
3130  double precision :: vtilL(ixi^s,2)
3131  double precision :: vtilR(ixi^s,2)
3132  double precision :: btilL(s%ixgs^s,ndim)
3133  double precision :: btilR(s%ixgs^s,ndim)
3134  double precision :: cp(ixi^s,2)
3135  double precision :: cm(ixi^s,2)
3136  double precision :: circ(ixi^s,1:ndim)
3137  ! current on cell edges
3138  double precision :: jce(ixi^s,7-2*ndim:3)
3139  integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
3140  integer :: idim1,idim2,idir
3141 
3142  associate(bfaces=>s%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
3143  cbarmax=>vcts%cbarmax)
3144 
3145  ! Calculate contribution to FEM of each edge,
3146  ! that is, estimate value of line integral of
3147  ! electric field in the positive idir direction.
3148 
3149  ! Loop over components of electric field
3150 
3151  ! idir: electric field component we need to calculate
3152  ! idim1: directions in which we already performed the reconstruction
3153  ! idim2: directions in which we perform the reconstruction
3154 
3155  ! if there is resistivity, get eta J
3156  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
3157 
3158  fe=zero
3159 
3160  do idir=7-2*ndim,3
3161  ! Indices
3162  ! idir: electric field component
3163  ! idim1: one surface
3164  ! idim2: the other surface
3165  ! cyclic permutation: idim1,idim2,idir=1,2,3
3166  ! Velocity components on the surface
3167  ! follow cyclic premutations:
3168  ! Sx(1),Sx(2)=y,z ; Sy(1),Sy(2)=z,x ; Sz(1),Sz(2)=x,y
3169 
3170  ixcmax^d=ixomax^d;
3171  ixcmin^d=ixomin^d-1+kr(idir,^d);
3172 
3173  ! Set indices and directions
3174  idim1=mod(idir,3)+1
3175  idim2=mod(idir+1,3)+1
3176 
3177  jxc^l=ixc^l+kr(idim1,^d);
3178  ixcp^l=ixc^l+kr(idim2,^d);
3179 
3180  ! Reconstruct transverse transport velocities
3181  call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
3182  vtill(ixi^s,2),vtilr(ixi^s,2))
3183 
3184  call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
3185  vtill(ixi^s,1),vtilr(ixi^s,1))
3186 
3187  ! Reconstruct magnetic fields
3188  ! Eventhough the arrays are larger, reconstruct works with
3189  ! the limits ixG.
3190  call reconstruct(ixi^l,ixc^l,idim2,bfaces(ixi^s,idim1),&
3191  btill(ixi^s,idim1),btilr(ixi^s,idim1))
3192 
3193  call reconstruct(ixi^l,ixc^l,idim1,bfaces(ixi^s,idim2),&
3194  btill(ixi^s,idim2),btilr(ixi^s,idim2))
3195 
3196  ! Take the maximum characteristic
3197 
3198  cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
3199  cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
3200 
3201  cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
3202  cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
3203 
3204 
3205  ! Calculate eletric field
3206  fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
3207  + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
3208  - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
3209  /(cp(ixc^s,1)+cm(ixc^s,1)) &
3210  +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
3211  + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
3212  - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
3213  /(cp(ixc^s,2)+cm(ixc^s,2))
3214 
3215  ! add current component of electric field at cell edges E=-vxB+eta J
3216  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
3217 
3218  fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
3219 
3220  if (.not.slab) then
3221  where(abs(x(ixc^s,r_)+half*dxlevel(r_)).lt.1.0d-9)
3222  fe(ixc^s,idir)=zero
3223  end where
3224  end if
3225 
3226  end do
3227 
3228  ! allow user to change inductive electric field, especially for boundary driven applications
3229  if(associated(usr_set_electric_field)) &
3230  call usr_set_electric_field(ixi^l,ixo^l,qdt,fe,sct)
3231 
3232  circ(ixi^s,1:ndim)=zero
3233 
3234  ! Calculate circulation on each face: interal(fE dot dl)
3235 
3236  do idim1=1,ndim ! Coordinate perpendicular to face
3237  ixcmax^d=ixomax^d;
3238  ixcmin^d=ixomin^d-kr(idim1,^d);
3239  do idim2=1,ndim
3240  do idir=7-2*ndim,3 ! Direction of line integral
3241  ! Assemble indices
3242  hxc^l=ixc^l-kr(idim2,^d);
3243  ! Add line integrals in direction idir
3244  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
3245  +lvc(idim1,idim2,idir)&
3246  *(fe(ixc^s,idir)&
3247  -fe(hxc^s,idir))
3248  end do
3249  end do
3250  end do
3251 
3252  ! Divide by the area of the face to get dB/dt
3253  do idim1=1,ndim
3254  ixcmax^d=ixomax^d;
3255  ixcmin^d=ixomin^d-kr(idim1,^d);
3256  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
3257  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
3258  elsewhere
3259  circ(ixc^s,idim1)=zero
3260  end where
3261  ! Time update
3262  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
3263  end do
3264 
3265  end associate
3266  end subroutine update_faces_hll
3267 
3268  !> calculate eta J at cell edges
3269  subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
3271  use mod_usr_methods
3272  use mod_geometry
3273 
3274  integer, intent(in) :: ixI^L, ixO^L
3275  type(state), intent(in) :: sCT, s
3276  ! current on cell edges
3277  double precision :: jce(ixi^s,7-2*ndim:3)
3278 
3279  ! current on cell centers
3280  double precision :: jcc(ixi^s,7-2*ndim:3)
3281  ! location at cell faces
3282  double precision :: xs(ixgs^t,1:ndim)
3283  ! resistivity
3284  double precision :: eta(ixi^s)
3285  double precision :: gradi(ixgs^t)
3286  integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
3287 
3288  associate(x=>s%x,dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
3289  ! calculate current density at cell edges
3290  jce=0.d0
3291  do idim1=1,ndim
3292  do idim2=1,ndim
3293  do idir=7-2*ndim,3
3294  if (lvc(idim1,idim2,idir)==0) cycle
3295  ixcmax^d=ixomax^d;
3296  ixcmin^d=ixomin^d+kr(idir,^d)-1;
3297  ixbmax^d=ixcmax^d-kr(idir,^d)+1;
3298  ixbmin^d=ixcmin^d;
3299  ! current at transverse faces
3300  xs(ixb^s,:)=x(ixb^s,:)
3301  xs(ixb^s,idim2)=x(ixb^s,idim2)+half*dx(ixb^s,idim2)
3302  call gradientx(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi,.true.)
3303  if (lvc(idim1,idim2,idir)==1) then
3304  jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
3305  else
3306  jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
3307  end if
3308  end do
3309  end do
3310  end do
3311  ! get resistivity
3312  if(mhd_eta>zero)then
3313  jce(ixc^s,:)=jce(ixc^s,:)*mhd_eta
3314  else
3315  ixa^l=ixo^l^ladd1;
3316  call get_current(wct,ixi^l,ixo^l,idirmin,jcc)
3317  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,jcc,eta)
3318  ! calcuate eta on cell edges
3319  do idir=7-2*ndim,3
3320  ixcmax^d=ixomax^d;
3321  ixcmin^d=ixomin^d+kr(idir,^d)-1;
3322  jcc(ixc^s,idir)=0.d0
3323  {do ix^db=0,1\}
3324  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
3325  ixamin^d=ixcmin^d+ix^d;
3326  ixamax^d=ixcmax^d+ix^d;
3327  jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
3328  {end do\}
3329  jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
3330  jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
3331  enddo
3332  end if
3333 
3334  end associate
3335  end subroutine get_resistive_electric_field
3336 
3337  !> calculate cell-center values from face-center values
3338  subroutine mhd_face_to_center(ixO^L,s)
3340  ! Non-staggered interpolation range
3341  integer, intent(in) :: ixO^L
3342  type(state) :: s
3343 
3344  integer :: fxO^L, gxO^L, hxO^L, jxO^L, kxO^L, idim
3345 
3346  associate(w=>s%w, ws=>s%ws)
3347 
3348  ! calculate cell-center values from face-center values in 2nd order
3349  do idim=1,ndim
3350  ! Displace index to the left
3351  ! Even if ixI^L is the full size of the w arrays, this is ok
3352  ! because the staggered arrays have an additional place to the left.
3353  hxo^l=ixo^l-kr(idim,^d);
3354  ! Interpolate to cell barycentre using arithmetic average
3355  ! This might be done better later, to make the method less diffusive.
3356  w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
3357  (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
3358  +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
3359  end do
3360 
3361  ! calculate cell-center values from face-center values in 4th order
3362  !do idim=1,ndim
3363  ! gxO^L=ixO^L-2*kr(idim,^D);
3364  ! hxO^L=ixO^L-kr(idim,^D);
3365  ! jxO^L=ixO^L+kr(idim,^D);
3366 
3367  ! ! Interpolate to cell barycentre using fourth order central formula
3368  ! w(ixO^S,mag(idim))=(0.0625d0/s%surface(ixO^S,idim))*&
3369  ! ( -ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
3370  ! +9.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
3371  ! +9.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
3372  ! -ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) )
3373  !end do
3374 
3375  ! calculate cell-center values from face-center values in 6th order
3376  !do idim=1,ndim
3377  ! fxO^L=ixO^L-3*kr(idim,^D);
3378  ! gxO^L=ixO^L-2*kr(idim,^D);
3379  ! hxO^L=ixO^L-kr(idim,^D);
3380  ! jxO^L=ixO^L+kr(idim,^D);
3381  ! kxO^L=ixO^L+2*kr(idim,^D);
3382 
3383  ! ! Interpolate to cell barycentre using sixth order central formula
3384  ! w(ixO^S,mag(idim))=(0.00390625d0/s%surface(ixO^S,idim))* &
3385  ! ( +3.0d0*ws(fxO^S,idim)*s%surfaceC(fxO^S,idim) &
3386  ! -25.0d0*ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
3387  ! +150.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
3388  ! +150.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
3389  ! -25.0d0*ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) &
3390  ! +3.0d0*ws(kxO^S,idim)*s%surfaceC(kxO^S,idim) )
3391  !end do
3392 
3393  end associate
3394 
3395  end subroutine mhd_face_to_center
3396 
3397  !> calculate magnetic field from vector potential
3398  subroutine b_from_vector_potential(ixIs^L, ixI^L, ixO^L, ws, x)
3401 
3402  integer, intent(in) :: ixIs^L, ixI^L, ixO^L
3403  double precision, intent(inout) :: ws(ixis^s,1:nws)
3404  double precision, intent(in) :: x(ixi^s,1:ndim)
3405 
3406  double precision :: Adummy(ixi^s,1:3)
3407 
3408  call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, adummy)
3409 
3410  end subroutine b_from_vector_potential
3411 
3412 end module mod_mhd_phys
integer, parameter cylindrical
Definition: mod_geometry.t:9
type(iw_methods), dimension(max_nw) phys_iw_methods
Special methods defined per variable.
Definition: mod_physics.t:43
procedure(special_resistivity), pointer usr_special_resistivity
double precision function, dimension(ixo^s) mhd_mag_en(w, ixIL, ixOL)
Compute evolving magnetic energy.
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_gravity.t:74
procedure(sub_global_source), pointer phys_global_source
Definition: mod_physics.t:56
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
Definition: mod_mhd_phys.t:29
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition: mod_mhd_phys.t:66
subroutine update_faces_hll(ixIL, ixOL, qdt, fE, sCT, s)
update faces
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_mhd_phys.t:69
subroutine update_faces_average(ixIL, ixOL, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
logical small_values_use_primitive
Use primitive variables when correcting small values.
subroutine viscosity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
Definition: mod_viscosity.t:86
subroutine add_source_res2(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 5 point stencil (2 neighbours) in each direction, conservative.
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
subroutine mhd_angmomfix(fC, x, wnew, ixIL, ixOL, idim)
Definition: mod_mhd_phys.t:216
double precision, public mhd_eta
The MHD resistivity.
Definition: mod_mhd_phys.t:87
double precision unit_length
Physical scaling factor for length.
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition: mod_forest.t:32
double precision unit_density
Physical scaling factor for density.
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:14
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:289
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges...
Definition: mod_geometry.t:519
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
subroutine mhd_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
Definition: mod_geometry.t:333
logical need_global_cmax
need global maximal wave speed
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer ndir
Number of spatial dimensions (components) for vector variables.
logical, public, protected mhd_viscosity
Whether viscosity is added.
Definition: mod_mhd_phys.t:17
integer nstep
How many sub-steps the time integrator takes.
subroutine e_to_rhos(ixIL, ixOL, w, x)
Convert energy to entropy.
Definition: mod_mhd_phys.t:649
integer coordinate
Definition: mod_geometry.t:6
procedure(sub_check_w), pointer phys_check_w
Definition: mod_physics.t:58
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:65
logical, public, protected mhd_gravity
Whether gravity is added.
Definition: mod_mhd_phys.t:20
procedure(sub_small_values), pointer phys_handle_small_values
Definition: mod_physics.t:63
Module with basic grid data structures.
Definition: mod_forest.t:2
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
integer istep
Index of the sub-step in a multi-step time integrator.
character(len=std_len) typeboundspeed
Which type of TVDLF method to use.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
Definition: mod_mhd_phys.t:132
subroutine mhd_getv_hall(w, x, ixIL, ixOL, vHall)
procedure(sub_boundary_adjust), pointer phys_boundary_adjust
Definition: mod_physics.t:60
subroutine add_source_hyperres(qdt, ixIL, ixOL, wCT, w, x)
Add Hyper-resistive source to w within ixO Uses 9 point stencil (4 neighbours) in each direction...
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl...
Definition: mod_physics.t:22
type(mg_t) mg
Data structure containing the multigrid tree.
logical b0field
split magnetic field as background B0 field
subroutine mhd_get_csound(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
Definition: mod_mhd_phys.t:817
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor, but its use is kept to a minimum.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition: mod_mhd_phys.t:126
subroutine update_faces_contact(ixIL, ixOL, qdt, wp, fC, fE, sCT, s)
update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
double precision small_density
subroutine add_source_res1(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in each direction, non-conservative. If the fourthorder precompiler flag is set, uses fourth order central difference for the laplacian. Then the stencil is 5 (2 neighbours).
subroutine mhd_check_w(primitive, ixIL, ixOL, w, flag)
Definition: mod_mhd_phys.t:504
subroutine particles_init()
Initialize particle data and parameters.
Definition: mod_particles.t:17
logical, public, protected mhd_energy
Whether an energy equation is used.
Definition: mod_mhd_phys.t:8
logical, public, protected mhd_divb_4thorder
Whether divB is computed with a fourth order approximation.
Definition: mod_mhd_phys.t:108
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag)
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
Definition: mod_mhd_phys.t:39
procedure(sub_update_faces), pointer phys_update_faces
Definition: mod_physics.t:64
integer, public, protected psi_
Indices of the GLM psi.
Definition: mod_mhd_phys.t:75
subroutine mhd_boundary_adjust
integer nghostcells
Number of ghost cells surrounding a grid.
subroutine radiative_cooling_init(phys_gamma, He_abund)
Radiative cooling initialization.
double precision, public mhd_gamma
The adiabatic index.
Definition: mod_mhd_phys.t:81
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
Module with all the methods that users can customize in AMRVAC.
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
Definition: mod_mhd_phys.t:14
subroutine, public small_values_error(w, x, ixIL, ixOL, w_flag, subname)
subroutine add_source_b0split(qdt, ixIL, ixOL, wCT, w, x)
Source terms after split off time-independent magnetic field.
subroutine mhd_check_params
Definition: mod_mhd_phys.t:452
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
integer, dimension(:), allocatable, public, protected mag
Indices of the magnetic field.
Definition: mod_mhd_phys.t:72
logical, dimension(:), allocatable small_values_fix_iw
Whether to apply small value fixes to certain variables.
double precision small_pressure
subroutine mhd_get_p_total(w, x, ixIL, ixOL, p)
Calculate total pressure within ixO^L including magnetic pressure.
Definition: mod_mhd_phys.t:958
subroutine mhd_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
logical stagger_grid
True for using stagger grid.
procedure(sub_angmomfix), pointer phys_angmomfix
Definition: mod_physics.t:62
Thermal conduction for HD and MHD.
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:46
integer ierrmpi
A global MPI error return code.
integer ixm
the mesh range (within a block with ghost cells)
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
Definition: mod_geometry.t:574
subroutine, public mhd_get_v(w, x, ixIL, ixOL, v)
Calculate v vector.
Definition: mod_mhd_phys.t:685
character(len=std_len), dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_velocity
Physical scaling factor for velocity.
subroutine mhd_write_info(fh)
Write this module&#39;s parameters to a snapsoht.
Definition: mod_mhd_phys.t:199
subroutine mhd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities.
Definition: mod_mhd_phys.t:727
procedure(set_electric_field), pointer usr_set_electric_field
subroutine mhd_getdt_hall(w, x, ixIL, ixOL, dxD, dthall)
subroutine, public get_current(w, ixIL, ixOL, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
subroutine, public mhd_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
Definition: mod_mhd_phys.t:534
subroutine mhd_get_csound_prim(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
Definition: mod_mhd_phys.t:867
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Definition: mod_mhd_phys.t:35
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_mhd_phys.t:63
double precision, public mhd_etah
TODO: what is this?
Definition: mod_mhd_phys.t:93
integer, parameter unitpar
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision, dimension(:,:), allocatable dx
integer, parameter flux_default
Indicates a normal flux.
Definition: mod_physics.t:29
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
Definition: mod_mhd_phys.t:129
character(len=std_len) typegrad
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
Definition: mod_mhd_phys.t:102
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
subroutine mhd_update_faces(ixIL, ixOL, qdt, wprim, fC, fE, sCT, s)
integer, dimension(:), allocatable, parameter d
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:448
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:51
integer mype
The rank of the current MPI task.
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
subroutine mhd_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
Definition: mod_mhd_phys.t:713
subroutine gravity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
Definition: mod_gravity.t:37
subroutine reconstruct(ixIL, ixCL, idir, q, qL, qR)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
Module for including gravity in (magneto)hydrodynamics simulations.
Definition: mod_gravity.t:2
subroutine mhd_physical_units()
Definition: mod_mhd_phys.t:475
subroutine rhos_to_e(ixIL, ixOL, w, x)
Convert entropy to energy.
Definition: mod_mhd_phys.t:667
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:53
subroutine internal_energy_add_source(qdt, ixIL, ixOL, wCT, w, x)
subroutine get_lorentz(ixIL, ixOL, w, JxB)
Compute the Lorentz force (JxB)
character(len=std_len), public, protected type_ct
Method type of constrained transport.
Definition: mod_mhd_phys.t:105
Module for handling problematic values in simulations, such as negative pressures.
double precision unit_pressure
Physical scaling factor for pressure.
double precision unit_magneticfield
Physical scaling factor for magnetic field.
logical use_particles
Use particles module or not.
double precision unit_temperature
Physical scaling factor for temperature.
subroutine mhd_modify_wlr(ixIL, ixOL, wLC, wRC, wLp, wRp, s, idir)
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_mhd_phys.t:60
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:54
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
subroutine, public mhd_get_csound2(w, x, ixIL, ixOL, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
Definition: mod_mhd_phys.t:942
Module containing all the particle routines.
Definition: mod_particles.t:2
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
subroutine, public mhd_face_to_center(ixOL, s)
calculate cell-center values from face-center values
double precision function, dimension(ixo^s), public mhd_kin_en(w, ixIL, ixOL, inv_rho)
compute kinetic energy
procedure(sub_check_params), pointer phys_check_params
Definition: mod_physics.t:45
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
Definition: mod_mhd_phys.t:90
double precision, dimension(ndim) dxlevel
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user&#39;s field
subroutine mhd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
Definition: mod_mhd_phys.t:78
subroutine gravity_init()
Initialize the module.
Definition: mod_gravity.t:27
logical, public, protected mhd_particles
Whether particles module is added.
Definition: mod_mhd_phys.t:26
integer phys_wider_stencil
To use wider stencils in flux calculations. A value of 1 will extend it by one cell in both direction...
Definition: mod_physics.t:18
subroutine, public mhd_get_pthermal(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L.
Definition: mod_mhd_phys.t:919
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
subroutine gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
Definition: mod_geometry.t:390
logical use_multigrid
Use multigrid (only available in 2D and 3D)
subroutine mhd_handle_small_values(primitive, w, x, ixIL, ixOL, subname)
Definition: mod_mhd_phys.t:598
logical slab
Cartesian geometry or not.
subroutine, public mhd_phys_init()
Definition: mod_mhd_phys.t:229
integer icomm
The MPI communicator.
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixIL, ixOL)
Compute 2 times total magnetic energy.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical, public, protected mhd_4th_order
MHD fourth order.
Definition: mod_mhd_phys.t:54
procedure(sub_get_v_idim), pointer phys_get_v_idim
Definition: mod_physics.t:52
subroutine mhd_gamma2_alfven(ixIL, ixOL, w, gamma_A2)
Compute 1/(1+v_A^2/c^2) for Boris&#39; approximation, where v_A is the Alfven velocity.
character(len=std_len) typediv
logical, public divbwave
Add divB wave in Roe solver.
Definition: mod_mhd_phys.t:123
integer, public, protected mhd_n_tracer
Number of tracer species.
Definition: mod_mhd_phys.t:57
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:49
subroutine boris_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
The module add viscous source terms and check time step.
Definition: mod_viscosity.t:10
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_add_source), pointer phys_add_source
Definition: mod_physics.t:55
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:47
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
Definition: mod_viscosity.t:57
subroutine thermal_conduction_init(phys_gamma)
Initialize the module.
integer, parameter spherical
Definition: mod_geometry.t:10
module radiative cooling – add optically thin radiative cooling for HD and MHD
double precision unit_time
Physical scaling factor for time.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
Definition: mod_mhd_phys.t:11
logical, public, protected mhd_glm
Whether GLM-MHD is used.
Definition: mod_mhd_phys.t:32
logical phys_energy
Solve energy equation or not.
logical, public, protected b0field_forcefree
B0 field is force-free.
Definition: mod_mhd_phys.t:135
integer, dimension(:, :), allocatable flux_type
Array per direction per variable, which can be used to specify that certain fluxes have to be treated...
Definition: mod_physics.t:26
logical, dimension(:), allocatable phyboundblock
True if a block has any physical boundary.
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
Definition: mod_mhd_phys.t:2
character(len=20), public small_values_method
How to handle small values.
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision function, dimension(ixo^s) mhd_gamma_alfven(w, ixIL, ixOL)
Compute 1/sqrt(1+v_A^2/c^2) for Boris simplification, where v_A is the Alfven velocity.
double precision, public mhd_adiab
The adiabatic constant.
Definition: mod_mhd_phys.t:84
double precision unit_numberdensity
Physical scaling factor for number density.
logical, public, protected mhd_hall
Whether Hall-MHD is used.
Definition: mod_mhd_phys.t:23
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition: mod_physics.t:59
subroutine, public mhd_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
Definition: mod_mhd_phys.t:564
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:50
procedure(sub_modify_wlr), pointer phys_modify_wlr
Definition: mod_physics.t:48
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:61
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
Definition: mod_physics.t:31
subroutine mhd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L.
Definition: mod_mhd_phys.t:973
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
double precision function, dimension(ixo^s) mhd_mag_i_all(w, ixIL, ixOL, idir)
Compute full magnetic field by direction.
subroutine, public mhd_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate v component.
Definition: mod_mhd_phys.t:701