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  !> type of TC used: 1: adapted module (mhd implementation), 2: adapted module (hd implementation)
14  integer, parameter, private :: mhd_tc =1
15  integer, parameter, private :: hd_tc =2
16  integer, protected :: use_mhd_tc = mhd_tc
17 
18  !> Whether radiative cooling is added
19  logical, public, protected :: mhd_radiative_cooling = .false.
20 
21  !> Whether viscosity is added
22  logical, public, protected :: mhd_viscosity = .false.
23 
24  !> Whether gravity is added
25  logical, public, protected :: mhd_gravity = .false.
26 
27  !> Whether Hall-MHD is used
28  logical, public, protected :: mhd_hall = .false.
29 
30  !> Whether Ambipolar term is used
31  logical, public, protected :: mhd_ambipolar = .false.
32 
33  !> Whether Ambipolar term is implemented using supertimestepping
34  logical, public, protected :: mhd_ambipolar_sts = .false.
35 
36  !> Whether Ambipolar term is implemented explicitly
37  logical, public, protected :: mhd_ambipolar_exp = .false.
38 
39  !> Whether particles module is added
40  logical, public, protected :: mhd_particles = .false.
41 
42  !> Whether magnetofriction is added
43  logical, public, protected :: mhd_magnetofriction = .false.
44 
45  !> Whether GLM-MHD is used
46  logical, public, protected :: mhd_glm = .false.
47 
48  !> Whether TRAC method is used
49  logical, public, protected :: mhd_trac = .false.
50 
51  !> Which TRAC method is used
52  integer, public, protected :: mhd_trac_type=1
53 
54  !> Height of the mask used in the TRAC method
55  double precision, public, protected :: mhd_trac_mask = 0.d0
56 
57  !> Distance between two adjacent traced magnetic field lines (in finest cell size)
58  integer, public, protected :: mhd_trac_finegrid=4
59 
60  !> Whether auxiliary internal energy is solved
61  logical, public, protected :: mhd_solve_eaux = .false.
62 
63  !> Whether internal energy is solved instead of total energy
64  logical, public, protected :: mhd_internal_e = .false.
65 
66  !> Whether divB cleaning sources are added splitting from fluid solver
67  logical, public, protected :: source_split_divb = .false.
68 
69  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
70  !> taking values within [0, 1]
71  double precision, public :: mhd_glm_alpha = 0.5d0
72 
73  !> Use Boris approximation
74  character(len=20) :: mhd_boris_method = "none"
75 
76  integer, parameter :: boris_none = 0
77  integer, parameter :: boris_reduced_force = 1
78  integer, parameter :: boris_simplification = 2
79  integer :: mhd_boris_type = boris_none
80 
81  !> Speed of light for Boris' approximation. If negative, test changes to the
82  !> momentum equation with gamma_A = 1
83  double precision :: mhd_boris_c = 0.0d0
84 
85  !> MHD fourth order
86  logical, public, protected :: mhd_4th_order = .false.
87 
88  !> Number of tracer species
89  integer, public, protected :: mhd_n_tracer = 0
90 
91  !> Index of the density (in the w array)
92  integer, public, protected :: rho_
93 
94  !> Indices of the momentum density
95  integer, allocatable, public, protected :: mom(:)
96 
97  !> Index of the energy density (-1 if not present)
98  integer, public, protected :: e_
99 
100  !> Index of the gas pressure (-1 if not present) should equal e_
101  integer, public, protected :: p_
102 
103  !> Indices of the magnetic field
104  integer, allocatable, public, protected :: mag(:)
105 
106  !> Indices of the GLM psi
107  integer, public, protected :: psi_
108 
109  !> Indices of auxiliary internal energy
110  integer, public, protected :: eaux_
111  integer, public, protected :: paux_
112 
113  !> Index of the cutoff temperature for the TRAC method
114  integer, public, protected :: tcoff_
115  integer, public, protected :: tweight_
116 
117  !> Indices of the tracers
118  integer, allocatable, public, protected :: tracer(:)
119 
120  !> The adiabatic index
121  double precision, public :: mhd_gamma = 5.d0/3.0d0
122 
123  !> The adiabatic constant
124  double precision, public :: mhd_adiab = 1.0d0
125 
126  !> The MHD resistivity
127  double precision, public :: mhd_eta = 0.0d0
128 
129  !> The MHD hyper-resistivity
130  double precision, public :: mhd_eta_hyper = 0.0d0
131 
132  !> TODO: what is this?
133  double precision, public :: mhd_etah = 0.0d0
134 
135  !> The MHD ambipolar coefficient
136  double precision, public :: mhd_eta_ambi = 0.0d0
137 
138  !> The small_est allowed energy
139  double precision, protected :: small_e
140 
141  !> The number of waves
142  integer :: nwwave=8
143 
144  !> Method type to clean divergence of B
145  character(len=std_len), public, protected :: typedivbfix = 'linde'
146 
147  !> Method type of constrained transport
148  character(len=std_len), public, protected :: type_ct = 'uct_contact'
149 
150  !> Whether divB is computed with a fourth order approximation
151  logical, public, protected :: mhd_divb_4thorder = .false.
152 
153  !> Method type in a integer for good performance
154  integer :: type_divb
155 
156  !> Coefficient of diffusive divB cleaning
157  double precision :: divbdiff = 0.8d0
158 
159  !> Update all equations due to divB cleaning
160  character(len=std_len) :: typedivbdiff = 'all'
161 
162  !> Use a compact way to add resistivity
163  logical :: compactres = .false.
164 
165  !> Add divB wave in Roe solver
166  logical, public :: divbwave = .true.
167 
168  !> clean initial divB
169  logical, public :: clean_initial_divb = .false.
170 
171  !> Helium abundance over Hydrogen
172  double precision, public, protected :: he_abundance=0.1d0
173 
174  !> To control divB=0 fix for boundary
175  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
176 
177  !> To skip * layer of ghost cells during divB=0 fix for boundary
178  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
179 
180  !> B0 field is force-free
181  logical, public, protected :: b0field_forcefree=.true.
182 
183  !> Whether an total energy equation is used
184  logical :: total_energy = .true.
185 
186  !> gamma minus one and its inverse
187  double precision :: gamma_1, inv_gamma_1
188 
189  ! DivB cleaning methods
190  integer, parameter :: divb_none = 0
191  integer, parameter :: divb_multigrid = -1
192  integer, parameter :: divb_glm = 1
193  integer, parameter :: divb_powel = 2
194  integer, parameter :: divb_janhunen = 3
195  integer, parameter :: divb_linde = 4
196  integer, parameter :: divb_lindejanhunen = 5
197  integer, parameter :: divb_lindepowel = 6
198  integer, parameter :: divb_lindeglm = 7
199  integer, parameter :: divb_ct = 8
200 
201 
202  ! Public methods
203  public :: mhd_phys_init
204  public :: mhd_kin_en
205  public :: mhd_get_pthermal
206  public :: mhd_get_v
207  public :: mhd_get_v_idim
208  public :: mhd_to_conserved
209  public :: mhd_to_primitive
210  public :: mhd_get_csound2
211  public :: mhd_face_to_center
212  public :: get_divb
213  public :: get_current
214  !> needed public if we want to use the ambipolar coefficient in the user file
215  public :: multiplyambicoef
216  public :: get_normalized_divb
217  public :: b_from_vector_potential
218  public :: mhd_mag_en_all
219  {^nooned
220  public :: mhd_clean_divb_multigrid
221  }
222 
223 
224  !define the subroutine interface for the ambipolar mask
225  abstract interface
226 
227  subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
229  integer, intent(in) :: ixi^l, ixo^l
230  double precision, intent(in) :: x(ixi^s,1:ndim)
231  double precision, intent(in) :: w(ixi^s,1:nw)
232  double precision, intent(inout) :: res(ixi^s)
233  end subroutine mask_subroutine
234 
235  end interface
236 
237  procedure(mask_subroutine), pointer :: usr_mask_ambipolar => null()
238  public :: usr_mask_ambipolar
239 
240 contains
241 
242  !> Read this module"s parameters from a file
243  subroutine mhd_read_params(files)
245  use mod_particles, only: particles_eta, particles_etah
246  character(len=*), intent(in) :: files(:)
247  integer :: n
248 
249  namelist /mhd_list/ mhd_energy, mhd_n_tracer, mhd_gamma, mhd_adiab,&
253  typedivbdiff, type_ct, compactres, divbwave, he_abundance, si_unit, b0field,&
255  particles_eta, particles_etah,&
257  mhd_boris_method, mhd_boris_c, clean_initial_divb, mhd_solve_eaux, mhd_internal_e, &
259 
260  do n = 1, size(files)
261  open(unitpar, file=trim(files(n)), status="old")
262  read(unitpar, mhd_list, end=111)
263 111 close(unitpar)
264  end do
265 
266  end subroutine mhd_read_params
267 
268  !> Write this module's parameters to a snapsoht
269  subroutine mhd_write_info(fh)
271  integer, intent(in) :: fh
272  integer, parameter :: n_par = 1
273  double precision :: values(n_par)
274  character(len=name_len) :: names(n_par)
275  integer, dimension(MPI_STATUS_SIZE) :: st
276  integer :: er
277 
278  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
279 
280  names(1) = "gamma"
281  values(1) = mhd_gamma
282  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
283  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
284  end subroutine mhd_write_info
285 
286  subroutine mhd_angmomfix(fC,x,wnew,ixI^L,ixO^L,idim)
288  double precision, intent(in) :: x(ixI^S,1:ndim)
289  double precision, intent(inout) :: fC(ixI^S,1:nwflux,1:ndim), wnew(ixI^S,1:nw)
290  integer, intent(in) :: ixI^L, ixO^L
291  integer, intent(in) :: idim
292  integer :: hxO^L, kxC^L, iw
293  double precision :: inv_volume(ixI^S)
294 
295  call mpistop("to do")
296 
297  end subroutine mhd_angmomfix
298 
299  subroutine mhd_phys_init()
303  use mod_viscosity, only: viscosity_init
304  use mod_gravity, only: gravity_init
305  use mod_particles, only: particles_init, particles_eta, particles_etah
307  use mod_physics
309 
310  {^nooned
312  }
313 
314  integer :: itr, idir
315 
316  call mhd_read_params(par_files)
317 
318  if(mhd_internal_e.and.mhd_solve_eaux) then
319  mhd_solve_eaux=.false.
320  if(mype==0) write(*,*) 'WARNING: set mhd_solve_eaux=F when mhd_internal_e=T'
321  end if
322 
323  if(.not. mhd_energy) then
324  if(mhd_internal_e) then
325  mhd_internal_e=.false.
326  if(mype==0) write(*,*) 'WARNING: set mhd_internal_e=F when mhd_energy=F'
327  end if
328  if(mhd_solve_eaux) then
329  mhd_solve_eaux=.false.
330  if(mype==0) write(*,*) 'WARNING: set mhd_solve_eaux=F when mhd_energy=F'
331  end if
332  if(mhd_thermal_conduction) then
333  mhd_thermal_conduction=.false.
334  if(mype==0) write(*,*) 'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
335  end if
336  if(mhd_radiative_cooling) then
337  mhd_radiative_cooling=.false.
338  if(mype==0) write(*,*) 'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
339  end if
340  if(mhd_trac) then
341  mhd_trac=.false.
342  if(mype==0) write(*,*) 'WARNING: set mhd_trac=F when mhd_energy=F'
343  end if
344  end if
345 
346  physics_type = "mhd"
352 
353  if(mhd_energy.and..not.mhd_internal_e) then
354  total_energy=.true.
355  else
356  total_energy=.false.
357  end if
358  phys_total_energy=total_energy
359 
360  {^ifoned
361  if(mhd_trac .and. mhd_trac_type .gt. 2) then
362  mhd_trac_type=1
363  if(mype==0) write(*,*) 'WARNING: reset mhd_trac_type=1 for 1D simulation'
364  end if
365  }
366  if(mhd_trac .and. mhd_trac_type .le. 4) then
367  mhd_trac_mask=bigdouble
368  if(mype==0) write(*,*) 'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
369  end if
371 
373 
374  ! set default gamma for polytropic/isothermal process
376  if(ndim==1) typedivbfix='none'
377  select case (typedivbfix)
378  case ('none')
379  type_divb = divb_none
380  {^nooned
381  case ('multigrid')
382  type_divb = divb_multigrid
383  use_multigrid = .true.
384  mg%operator_type = mg_laplacian
386  }
387  case ('glm')
388  mhd_glm = .true.
389  need_global_cmax = .true.
390  type_divb = divb_glm
391  case ('powel', 'powell')
392  type_divb = divb_powel
393  case ('janhunen')
394  type_divb = divb_janhunen
395  case ('linde')
396  type_divb = divb_linde
397  case ('lindejanhunen')
398  type_divb = divb_lindejanhunen
399  case ('lindepowel')
400  type_divb = divb_lindepowel
401  case ('lindeglm')
402  mhd_glm = .true.
403  need_global_cmax = .true.
404  type_divb = divb_lindeglm
405  case ('ct')
406  type_divb = divb_ct
407  stagger_grid = .true.
408  case default
409  call mpistop('Unknown divB fix')
410  end select
411 
412  ! Determine flux variables
413  rho_ = var_set_rho()
414 
415  allocate(mom(ndir))
416  mom(:) = var_set_momentum(ndir)
417 
418  ! Set index of energy variable
419  if (mhd_energy) then
420  nwwave = 8
421  e_ = var_set_energy() ! energy density
422  p_ = e_ ! gas pressure
423  else
424  nwwave = 7
425  e_ = -1
426  p_ = -1
427  end if
428 
429  allocate(mag(ndir))
430  mag(:) = var_set_bfield(ndir)
431 
432  ! set auxiliary internal energy variable
433  if(mhd_energy .and. mhd_solve_eaux) then
434  eaux_ = var_set_internal_energy()
435  paux_ = eaux_
436  else
437  eaux_ = -1
438  paux_ = -1
439  end if
440 
441  if (mhd_glm) then
442  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
443  else
444  psi_ = -1
445  end if
446 
447  allocate(tracer(mhd_n_tracer))
448 
449  ! Set starting index of tracers
450  do itr = 1, mhd_n_tracer
451  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
452  end do
453 
454  ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
455  tweight_ = -1
456  if(mhd_trac) then
457  tcoff_ = var_set_wextra()
458  iw_tcoff=tcoff_
459  if(mhd_trac_type .ge. 3) then
460  tweight_ = var_set_wextra()
461  endif
462  else
463  tcoff_ = -1
464  end if
465 
466  ! set number of variables which need update ghostcells
467  nwgc=nwflux
468 
469  ! determine number of stagger variables
470  if(stagger_grid) nws=ndim
471 
472  nvector = 2 ! No. vector vars
473  allocate(iw_vector(nvector))
474  iw_vector(1) = mom(1) - 1 ! TODO: why like this?
475  iw_vector(2) = mag(1) - 1 ! TODO: why like this?
476 
477  ! Check whether custom flux types have been defined
478  if (.not. allocated(flux_type)) then
479  allocate(flux_type(ndir, nwflux))
481  else if (any(shape(flux_type) /= [ndir, nwflux])) then
482  call mpistop("phys_check error: flux_type has wrong shape")
483  end if
484 
485  if(nwflux>mag(ndir)) then
486  ! for flux of eaux and tracers, using hll flux
487  flux_type(:,mag(ndir)+1:nwflux)=flux_hll
488  end if
489 
490  if(ndim>1) then
491  if(mhd_glm) then
493  do idir=1,ndir
494  flux_type(idir,mag(idir))=flux_special
495  end do
496  else
497  do idir=1,ndir
498  flux_type(idir,mag(idir))=flux_tvdlf
499  end do
500  end if
501  end if
502 
503  select case (mhd_boris_method)
504  case ("none")
505  mhd_boris_type = boris_none
506  case ("reduced_force")
507  mhd_boris_type = boris_reduced_force
508  case ("simplification")
509  mhd_boris_type = boris_simplification
510  do idir = 1, ndir
511  phys_iw_methods(mom(idir))%inv_capacity => mhd_gamma2_alfven
512  end do
513  case default
514  call mpistop("Unknown mhd_boris_method (none, reduced_force, simplification)")
515  end select
516 
529  if(mhd_solve_eaux) then
532  else
535  end if
543 
544  if(type_divb==divb_glm) then
546  end if
547 
548  ! if using ct stagger grid, boundary divb=0 is not done here
549  if(stagger_grid) then
554  else if(ndim>1) then
556  end if
557 
558  {^nooned
559  ! clean initial divb
561  }
562 
563  ! Whether diagonal ghost cells are required for the physics
564  if(type_divb < divb_linde) phys_req_diagonal = .false.
565 
566  ! derive units from basic units
567  call mhd_physical_units()
568 
569  if(.not. mhd_energy .and. mhd_thermal_conduction) then
570  call mpistop("thermal conduction needs mhd_energy=T")
571  end if
572  if(.not. mhd_energy .and. mhd_radiative_cooling) then
573  call mpistop("radiative cooling needs mhd_energy=T")
574  end if
575 
576  ! resistive MHD needs diagonal ghost cells
577  if(mhd_eta/=0.d0) phys_req_diagonal = .true.
578 
579  ! initialize thermal conduction module
580  if (mhd_thermal_conduction) then
581  phys_req_diagonal = .true.
582  if(use_mhd_tc .eq. mhd_tc) then
583 
584  if(mhd_internal_e) then
586  else
587  if(mhd_solve_eaux) then
589  else
591  endif
592  endif
593 
594  else if(use_mhd_tc .eq. hd_tc) then
595  if(mhd_internal_e) then
597  else
598  if(mhd_solve_eaux) then
600  else
602  endif
603  endif
604  endif
605  end if
606 
607  ! Initialize radiative cooling module
608  if (mhd_radiative_cooling) then
610  end if
611 
612  ! Initialize viscosity module
614 
615  ! Initialize gravity module
616  if(mhd_gravity) then
617  call gravity_init()
618  end if
619 
620  ! Initialize particles module
621  if(mhd_particles) then
622  call particles_init()
623  phys_req_diagonal = .true.
624  if(mype==0) then
625  write(*,*) '*****Using particles: with mhd_eta, mhd_etah :', mhd_eta, mhd_etah
626  write(*,*) '*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
627  end if
628  end if
629 
630  ! initialize magnetofriction module
631  if(mhd_magnetofriction) then
632  phys_req_diagonal = .true.
633  call magnetofriction_init()
634  end if
635 
636  ! For Hall, we need one more reconstructed layer since currents are computed
637  ! in mhd_get_flux: assuming one additional ghost layer (two for FOURTHORDER) was
638  ! added in nghostcells.
639  if(mhd_hall) then
640  phys_req_diagonal = .true.
641  if(mhd_4th_order) then
643  else
645  end if
646  end if
647 
648  if(mhd_ambipolar) then
649  phys_req_diagonal = .true.
650  if(mhd_ambipolar_sts) then
651  call sts_init()
652  if(mhd_internal_e) then
654  ndir,mag(1),ndir,.true.)
655  else
657  mag(ndir)-mom(ndir),mag(1),ndir,.true.)
658  end if
659  else
660  mhd_ambipolar_exp=.true.
661  ! For flux ambipolar term, we need one more reconstructed layer since currents are computed
662  ! in mhd_get_flux: assuming one additional ghost layer (two for FOURTHORDER) was
663  ! added in nghostcells.
664  if(mhd_4th_order) then
666  else
668  end if
669  end if
670  end if
671 
672  end subroutine mhd_phys_init
673 
674  subroutine mhd_check_params
676 
677  ! after user parameter setting
678  gamma_1=mhd_gamma-1.d0
679  if (.not. mhd_energy) then
680  if (mhd_gamma <= 0.0d0) call mpistop ("Error: mhd_gamma <= 0")
681  if (mhd_adiab < 0.0d0) call mpistop ("Error: mhd_adiab < 0")
683  else
684  if (mhd_gamma <= 0.0d0 .or. mhd_gamma == 1.0d0) &
685  call mpistop ("Error: mhd_gamma <= 0 or mhd_gamma == 1")
686  inv_gamma_1=1.d0/gamma_1
687  small_e = small_pressure * inv_gamma_1
688  end if
689 
690  if (mhd_boris_type > 0 .and. abs(mhd_boris_c) <= 0.0d0) then
691  call mpistop("You have not specified mhd_boris_c")
692  end if
693 
694  end subroutine mhd_check_params
695 
696  subroutine mhd_physical_units()
698  double precision :: mp,kB,miu0,c_lightspeed
699  ! Derive scaling units
700  if(si_unit) then
701  mp=mp_si
702  kb=kb_si
703  miu0=miu0_si
704  c_lightspeed=c_si
705  else
706  mp=mp_cgs
707  kb=kb_cgs
708  miu0=4.d0*dpi ! G^2 cm^2 dyne^-1
709  c_lightspeed=const_c
710  end if
711  if(unit_numberdensity/=1.d0) then
713  else if(unit_density/=1.d0) then
715  end if
716  if(unit_temperature/=1.d0) then
720  else if(unit_velocity/=1.d0) then
724  else if(unit_pressure/=1.d0) then
728  else if(unit_magneticfield/=1.d0) then
732  end if
733  if(unit_length/=1.d0) then
735  else if(unit_time/=1.d0) then
737  end if
738  ! Additional units needed for the particles
739  c_norm=c_lightspeed/unit_velocity
741  if (.not. si_unit) unit_charge = unit_charge*const_c
743 
744  end subroutine mhd_physical_units
745 
746  subroutine mhd_check_w(primitive,ixI^L,ixO^L,w,flag)
748 
749  logical, intent(in) :: primitive
750  integer, intent(in) :: ixI^L, ixO^L
751  double precision, intent(in) :: w(ixI^S,nw)
752  double precision :: tmp(ixI^S)
753  logical, intent(inout) :: flag(ixI^S,1:nw)
754 
755  flag=.false.
756  where(w(ixo^s, rho_) < small_density) flag(ixo^s,rho_) = .true.
757 
758  if(mhd_energy) then
759  if(primitive) then
760  where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
761  else
762  if(mhd_internal_e) then
763  where(w(ixo^s, e_) < small_e) flag(ixo^s,e_) = .true.
764  else
765  tmp(ixo^s)=w(ixo^s,e_)-&
766  mhd_kin_en(w,ixi^l,ixo^l)-mhd_mag_en(w,ixi^l,ixo^l)
767  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
768  end if
769  end if
770  end if
771 
772  end subroutine mhd_check_w
773 
774  !> Transform primitive variables into conservative ones
775  subroutine mhd_to_conserved(ixI^L,ixO^L,w,x)
777  integer, intent(in) :: ixi^l, ixo^l
778  double precision, intent(inout) :: w(ixi^s, nw)
779  double precision, intent(in) :: x(ixi^s, 1:ndim)
780  integer :: idir, itr
781 
782  !if (fix_small_values) then
783  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
784  !end if
785 
786  ! Calculate total energy from pressure, kinetic and magnetic energy
787  if(mhd_energy) then
788  if(mhd_internal_e) then
789  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
790  else
791  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
792  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
793  +mhd_mag_en(w, ixi^l, ixo^l)
794  if(mhd_solve_eaux) w(ixo^s,eaux_)=w(ixo^s,paux_)*inv_gamma_1
795  end if
796  end if
797 
798  ! Convert velocity to momentum
799  do idir = 1, ndir
800  w(ixo^s, mom(idir)) = w(ixo^s, rho_) * w(ixo^s, mom(idir))
801  end do
802  end subroutine mhd_to_conserved
803 
804  !> Transform conservative variables into primitive ones
805  subroutine mhd_to_primitive(ixI^L,ixO^L,w,x)
807  integer, intent(in) :: ixi^l, ixo^l
808  double precision, intent(inout) :: w(ixi^s, nw)
809  double precision, intent(in) :: x(ixi^s, 1:ndim)
810  double precision :: inv_rho(ixo^s)
811  integer :: itr, idir
812 
813  if (fix_small_values) then
814  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive')
815  end if
816 
817  inv_rho=1.0d0/w(ixo^s,rho_)
818 
819  ! Calculate pressure = (gamma-1) * (e-ek-eb)
820  if(mhd_energy) then
821  if(mhd_internal_e) then
822  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
823  else
824  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
825  -mhd_kin_en(w,ixi^l,ixo^l,inv_rho)&
826  -mhd_mag_en(w,ixi^l,ixo^l))
827  if(mhd_solve_eaux) w(ixo^s,paux_)=w(ixo^s,eaux_)*gamma_1
828  end if
829  end if
830 
831  ! Convert momentum to velocity
832  do idir = 1, ndir
833  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
834  end do
835 
836  end subroutine mhd_to_primitive
837 
838  !> Transform internal energy to total energy
839  subroutine mhd_ei_to_e(ixI^L,ixO^L,w,x)
841  integer, intent(in) :: ixI^L, ixO^L
842  double precision, intent(inout) :: w(ixI^S, nw)
843  double precision, intent(in) :: x(ixI^S, 1:ndim)
844 
845  ! Calculate total energy from internal, kinetic and magnetic energy
846  w(ixo^s,e_)=w(ixo^s,e_)&
847  +mhd_kin_en(w,ixi^l,ixo^l)&
848  +mhd_mag_en(w,ixi^l,ixo^l)
849 
850  end subroutine mhd_ei_to_e
851 
852  !> Transform total energy to internal energy
853  subroutine mhd_e_to_ei(ixI^L,ixO^L,w,x)
855  integer, intent(in) :: ixI^L, ixO^L
856  double precision, intent(inout) :: w(ixI^S, nw)
857  double precision, intent(in) :: x(ixI^S, 1:ndim)
858 
859  ! Calculate ei = e - ek - eb
860  w(ixo^s,e_)=w(ixo^s,e_)&
861  -mhd_kin_en(w,ixi^l,ixo^l)&
862  -mhd_mag_en(w,ixi^l,ixo^l)
863 
864  if(fix_small_values) then
865  call mhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,'mhd_e_to_ei')
866  end if
867 
868  end subroutine mhd_e_to_ei
869 
870  !> Update eaux and transform internal energy to total energy
871  subroutine mhd_ei_to_e_aux(ixI^L,ixO^L,w,x)
873  integer, intent(in) :: ixI^L, ixO^L
874  double precision, intent(inout) :: w(ixI^S, nw)
875  double precision, intent(in) :: x(ixI^S, 1:ndim)
876 
877  w(ixo^s,eaux_)=w(ixo^s,e_)
878  ! Calculate total energy from internal, kinetic and magnetic energy
879  w(ixo^s,e_)=w(ixo^s,e_)&
880  +mhd_kin_en(w,ixi^l,ixo^l)&
881  +mhd_mag_en(w,ixi^l,ixo^l)
882 
883  end subroutine mhd_ei_to_e_aux
884 
885  !> Transform total energy to internal energy via eaux as internal energy
886  subroutine mhd_e_to_ei_aux(ixI^L,ixO^L,w,x)
888  integer, intent(in) :: ixI^L, ixO^L
889  double precision, intent(inout) :: w(ixI^S, nw)
890  double precision, intent(in) :: x(ixI^S, 1:ndim)
891 
892  w(ixo^s,e_)=w(ixo^s,eaux_)
893 
894  end subroutine mhd_e_to_ei_aux
895 
896  subroutine mhd_energy_synchro(ixI^L,ixO^L,w,x)
898  integer, intent(in) :: ixI^L,ixO^L
899  double precision, intent(in) :: x(ixI^S,1:ndim)
900  double precision, intent(inout) :: w(ixI^S,1:nw)
901 
902  double precision :: pth1(ixI^S),pth2(ixI^S),alfa(ixI^S),beta(ixI^S)
903  double precision, parameter :: beta_low=0.005d0,beta_high=0.05d0
904 
905 ! double precision :: vtot(ixI^S),cs2(ixI^S),mach(ixI^S)
906 ! double precision, parameter :: mach_low=20.d0,mach_high=200.d0
907 
908  ! get magnetic energy
909  alfa(ixo^s)=mhd_mag_en(w,ixi^l,ixo^l)
910  pth1(ixo^s)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l)-alfa(ixo^s))
911  pth2(ixo^s)=w(ixo^s,eaux_)*gamma_1
912  ! get plasma beta
913  beta(ixo^s)=min(pth1(ixo^s),pth2(ixo^s))/alfa(ixo^s)
914 
915  ! whether Mach number should be another criterion ?
916 ! vtot(ixO^S)=sum(w(ixO^S,mom(:))**2,dim=ndim+1)
917 ! call mhd_get_csound2(w,x,ixI^L,ixO^L,cs2)
918 ! mach(ixO^S)=sqrt(vtot(ixO^S)/cs2(ixO^S))/w(ixO^S,rho_)
919  where(beta(ixo^s) .ge. beta_high)
920 ! where(beta(ixO^S) .ge. beta_high .and. mach(ixO^S) .le. mach_low)
921  w(ixo^s,eaux_)=pth1(ixo^s)*inv_gamma_1
922  else where(beta(ixo^s) .le. beta_low)
923 ! else where(beta(ixO^S) .le. beta_low .or. mach(ixO^S) .ge. mach_high)
924  w(ixo^s,e_)=w(ixo^s,e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,eaux_)
925  else where
926  alfa(ixo^s)=dlog(beta(ixo^s)/beta_low)/dlog(beta_high/beta_low)
927 ! alfa(ixO^S)=min(dlog(beta(ixO^S)/beta_low)/dlog(beta_high/beta_low),
928 ! dlog(mach_high(ixO^S)/mach(ixO^S))/dlog(mach_high/mach_low))
929  w(ixo^s,eaux_)=(pth2(ixo^s)*(one-alfa(ixo^s))&
930  +pth1(ixo^s)*alfa(ixo^s))*inv_gamma_1
931  w(ixo^s,e_)=w(ixo^s,e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,eaux_)
932  end where
933  end subroutine mhd_energy_synchro
934 
935  subroutine mhd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
937  use mod_small_values
938  logical, intent(in) :: primitive
939  integer, intent(in) :: ixI^L,ixO^L
940  double precision, intent(inout) :: w(ixI^S,1:nw)
941  double precision, intent(in) :: x(ixI^S,1:ndim)
942  character(len=*), intent(in) :: subname
943 
944  integer :: idir
945  logical :: flag(ixI^S,1:nw)
946 
947  if(small_values_method == "ignore") return
948 
949  call mhd_check_w(primitive, ixi^l, ixo^l, w, flag)
950 
951  if(any(flag)) then
952  select case (small_values_method)
953  case ("replace")
954  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
955  do idir = 1, ndir
956  if(small_values_fix_iw(mom(idir))) then
957  where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
958  end if
959  end do
960 
961  if(mhd_energy) then
962  if(primitive) then
963  where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
964  else if(mhd_internal_e) then
965  where(flag(ixo^s,e_))
966  w(ixo^s,e_)=small_e
967  end where
968  else
969  where(flag(ixo^s,e_))
970  w(ixo^s,e_) = small_e+&
971  mhd_kin_en(w,ixi^l,ixo^l)+&
972  mhd_mag_en(w,ixi^l,ixo^l)
973  end where
974  if(mhd_solve_eaux) then
975  where(flag(ixo^s,e_))
976  w(ixo^s,eaux_)=small_e
977  end where
978  end if
979  end if
980  end if
981  case ("average")
982  call small_values_average(ixi^l, ixo^l, w, x, flag)
983  case default
984  if(.not.primitive) then
985  !convert w to primitive
986  ! Calculate pressure = (gamma-1) * (e-ek-eb)
987  if(mhd_energy) then
988  if(mhd_internal_e) then
989  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
990  else
991  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
992  -mhd_kin_en(w,ixi^l,ixo^l)&
993  -mhd_mag_en(w,ixi^l,ixo^l))
994  if(mhd_solve_eaux) w(ixo^s,paux_)=w(ixo^s,eaux_)*gamma_1
995  end if
996  end if
997  ! Convert momentum to velocity
998  do idir = 1, ndir
999  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/w(ixo^s,rho_)
1000  end do
1001  end if
1002  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
1003  end select
1004  end if
1005  end subroutine mhd_handle_small_values
1006 
1007  !> Convert energy to entropy
1008  subroutine e_to_rhos(ixI^L,ixO^L,w,x)
1010  integer, intent(in) :: ixI^L, ixO^L
1011  double precision,intent(inout) :: w(ixI^S,nw)
1012  double precision, intent(in) :: x(ixI^S,1:ndim)
1013 
1014  if(mhd_energy) then
1015  if(.not.mhd_internal_e) &
1016  w(ixo^s, e_)=w(ixo^s, e_)-mhd_kin_en(w, ixi^l, ixo^l) &
1017  -mhd_mag_en(w, ixi^l, ixo^l)
1018  w(ixo^s, e_)=gamma_1*w(ixo^s, rho_)**(1.0d0 - mhd_gamma)&
1019  *w(ixo^s, e_)
1020  else
1021  call mpistop("e_to_rhos can not be used without energy equation!")
1022  end if
1023  end subroutine e_to_rhos
1024 
1025  !> Convert entropy to energy
1026  subroutine rhos_to_e(ixI^L,ixO^L,w,x)
1028  integer, intent(in) :: ixI^L, ixO^L
1029  double precision :: w(ixI^S,nw)
1030  double precision, intent(in) :: x(ixI^S,1:ndim)
1031 
1032  if(mhd_energy) then
1033  w(ixo^s, e_)=w(ixo^s, rho_)**gamma_1 * w(ixo^s, e_) &
1034  * inv_gamma_1
1035  if(.not.mhd_internal_e) &
1036  w(ixo^s, e_)=w(ixo^s, e_)+mhd_kin_en(w, ixi^l, ixo^l) + &
1037  mhd_mag_en(w, ixi^l, ixo^l)
1038  else
1039  call mpistop("rhos_to_e can not be used without energy equation!")
1040  end if
1041  end subroutine rhos_to_e
1042 
1043  !> Calculate v vector
1044  subroutine mhd_get_v(w,x,ixI^L,ixO^L,v)
1046 
1047  integer, intent(in) :: ixi^l, ixo^l
1048  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
1049  double precision, intent(out) :: v(ixi^s,ndir)
1050 
1051  integer :: idir
1052 
1053  do idir=1,ndir
1054  v(ixo^s,idir) = w(ixo^s, mom(idir)) / w(ixo^s,rho_)
1055  end do
1056 
1057  end subroutine mhd_get_v
1058 
1059  !> Calculate v component
1060  subroutine mhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
1062 
1063  integer, intent(in) :: ixi^l, ixo^l, idim
1064  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
1065  double precision, intent(out) :: v(ixi^s)
1066 
1067  v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s,rho_)
1068 
1069  end subroutine mhd_get_v_idim
1070 
1071  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
1072  subroutine mhd_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1074 
1075  integer, intent(in) :: ixI^L, ixO^L, idim
1076  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1077  double precision, intent(inout) :: cmax(ixI^S)
1078 
1079  call mhd_get_csound(w,x,ixi^l,ixo^l,idim,cmax)
1080 
1081  cmax(ixo^s)=abs(w(ixo^s,mom(idim))/w(ixo^s,rho_))+cmax(ixo^s)
1082 
1083  end subroutine mhd_get_cmax
1084 
1085  subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
1087 
1088  integer, intent(in) :: ixI^L, ixO^L
1089  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1090  double precision, intent(inout) :: a2max(ndim)
1091  double precision :: a2(ixI^S,ndim,nw)
1092  integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
1093 
1094  a2=zero
1095  do i = 1,ndim
1096  !> 4th order
1097  hxo^l=ixo^l-kr(i,^d);
1098  gxo^l=hxo^l-kr(i,^d);
1099  jxo^l=ixo^l+kr(i,^d);
1100  kxo^l=jxo^l+kr(i,^d);
1101  a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1102  -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1103  a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
1104  end do
1105  end subroutine mhd_get_a2max
1106 
1107  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
1108  subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1110  use mod_geometry
1111  integer, intent(in) :: ixI^L,ixO^L
1112  double precision, intent(in) :: x(ixI^S,1:ndim)
1113  double precision, intent(inout) :: w(ixI^S,1:nw)
1114  double precision, intent(out) :: Tco_local,Tmax_local
1115 
1116  double precision, parameter :: trac_delta=0.25d0
1117  double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1118  double precision, dimension(ixI^S,1:ndir) :: bunitvec
1119  double precision, dimension(ixI^S,1:ndim) :: gradT
1120  double precision :: Bdir(ndim)
1121  double precision :: ltrc,ltrp,altr(ixI^S)
1122  integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
1123  integer :: jxP^L,hxP^L,ixP^L
1124  logical :: lrlt(ixI^S)
1125 
1126  if(mhd_internal_e) then
1127  tmp1(ixi^s)=w(ixi^s,e_)
1128  else
1129  tmp1(ixi^s)=w(ixi^s,e_)-0.5d0*(sum(w(ixi^s,iw_mom(:))**2,dim=ndim+1)/&
1130  w(ixi^s,rho_)+sum(w(ixi^s,iw_mag(:))**2,dim=ndim+1))
1131  end if
1132  te(ixi^s)=tmp1(ixi^s)/w(ixi^s,rho_)*(mhd_gamma-1.d0)
1133  tco_local=zero
1134  tmax_local=maxval(te(ixo^s))
1135 
1136  {^ifoned
1137  select case(mhd_trac_type)
1138  case(0)
1139  !> test case, fixed cutoff temperature
1140  block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
1141  case(1)
1142  hxo^l=ixo^l-1;
1143  jxo^l=ixo^l+1;
1144  lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1145  lrlt=.false.
1146  where(lts(ixo^s) > trac_delta)
1147  lrlt(ixo^s)=.true.
1148  end where
1149  if(any(lrlt(ixo^s))) then
1150  tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1151  end if
1152  case(2)
1153  !> iijima et al. 2021, LTRAC method
1154  ltrc=1.5d0
1155  ltrp=4.d0
1156  ixp^l=ixo^l^ladd1;
1157  hxo^l=ixo^l-1;
1158  jxo^l=ixo^l+1;
1159  hxp^l=ixp^l-1;
1160  jxp^l=ixp^l+1;
1161  lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1162  lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1163  lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
1164  block%wextra(ixo^s,tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
1165  case default
1166  call mpistop("mhd_trac_type not allowed for 1D simulation")
1167  end select
1168  }
1169  {^nooned
1170  select case(mhd_trac_type)
1171  case(0)
1172  !> test case, fixed cutoff temperature
1173  block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
1174  case(1,4,6)
1175  ! temperature gradient at cell centers
1176  do idims=1,ndim
1177  call gradient(te,ixi^l,ixo^l,idims,tmp1)
1178  gradt(ixo^s,idims)=tmp1(ixo^s)
1179  end do
1180  ! B vector
1181  if(b0field) then
1182  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
1183  else
1184  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1185  end if
1186  if(mhd_trac_type .gt. 1) then
1187  ! B direction at cell center
1188  bdir=zero
1189  {do ixa^d=0,1\}
1190  ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
1191  bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
1192  {end do\}
1193  if(sum(bdir(:)**2) .gt. zero) then
1194  bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1195  end if
1196  block%special_values(3:ndim+2)=bdir(1:ndim)
1197  end if
1198  tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1199  where(tmp1(ixo^s)/=0.d0)
1200  tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1201  elsewhere
1202  tmp1(ixo^s)=bigdouble
1203  end where
1204  ! b unit vector: magnetic field direction vector
1205  do idims=1,ndim
1206  bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1207  end do
1208  ! temperature length scale inversed
1209  lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1210  ! fraction of cells size to temperature length scale
1211  if(slab_uniform) then
1212  lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1213  else
1214  lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1215  end if
1216  lrlt=.false.
1217  where(lts(ixo^s) > trac_delta)
1218  lrlt(ixo^s)=.true.
1219  end where
1220  if(any(lrlt(ixo^s))) then
1221  block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1222  else
1223  block%special_values(1)=zero
1224  end if
1225  block%special_values(2)=tmax_local
1226  case(2)
1227  !> iijima et al. 2021, LTRAC method
1228  ltrc=1.5d0
1229  ltrp=4.d0
1230  ixp^l=ixo^l^ladd1;
1231  ! temperature gradient at cell centers
1232  do idims=1,ndim
1233  call gradient(te,ixi^l,ixp^l,idims,tmp1)
1234  gradt(ixp^s,idims)=tmp1(ixp^s)
1235  end do
1236  ! B vector
1237  if(b0field) then
1238  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
1239  else
1240  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
1241  end if
1242  tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
1243  where(tmp1(ixp^s)/=0.d0)
1244  tmp1(ixp^s)=1.d0/tmp1(ixp^s)
1245  elsewhere
1246  tmp1(ixp^s)=bigdouble
1247  end where
1248  ! b unit vector: magnetic field direction vector
1249  do idims=1,ndim
1250  bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
1251  end do
1252  ! temperature length scale inversed
1253  lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
1254  ! fraction of cells size to temperature length scale
1255  if(slab_uniform) then
1256  lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
1257  else
1258  lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
1259  end if
1260  lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1261 
1262  altr=zero
1263  do idims=1,ndim
1264  hxo^l=ixo^l-kr(idims,^d);
1265  jxo^l=ixo^l+kr(idims,^d);
1266  altr(ixo^s)=altr(ixo^s)+0.25d0*(lts(hxo^s)+two*lts(ixo^s)+lts(jxo^s))*bunitvec(ixo^s,idims)**2
1267  end do
1268  block%wextra(ixo^s,tcoff_)=te(ixo^s)*altr(ixo^s)**0.4d0
1269  ! need one ghost layer for thermal conductivity
1270  {block%wextra(ixomin^d-1^d%ixP^s,tcoff_)=block%wextra(ixomin^d^d%ixP^s,tcoff_) \}
1271  {block%wextra(ixomax^d+1^d%ixP^s,tcoff_)=block%wextra(ixomax^d^d%ixP^s,tcoff_) \}
1272  case(3,5)
1273  !> do nothing here
1274  case default
1275  call mpistop("unknown mhd_trac_type")
1276  end select
1277  }
1278  end subroutine mhd_get_tcutoff
1279 
1280  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
1281  subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
1283 
1284  integer, intent(in) :: ixI^L, ixO^L, idim
1285  double precision, intent(in) :: wprim(ixI^S, nw)
1286  double precision, intent(in) :: x(ixI^S,1:ndim)
1287  double precision, intent(out) :: Hspeed(ixI^S)
1288 
1289  double precision :: csound(ixI^S,ndim),tmp(ixI^S)
1290  integer :: jxC^L, ixC^L, ixA^L, id, ix^D
1291 
1292  hspeed=0.d0
1293  ixa^l=ixo^l^ladd1;
1294  do id=1,ndim
1295  call mhd_get_csound_prim(wprim,x,ixi^l,ixa^l,id,tmp)
1296  csound(ixa^s,id)=tmp(ixa^s)
1297  end do
1298  ixcmax^d=ixomax^d;
1299  ixcmin^d=ixomin^d+kr(idim,^d)-1;
1300  jxcmax^d=ixcmax^d+kr(idim,^d);
1301  jxcmin^d=ixcmin^d+kr(idim,^d);
1302  hspeed(ixc^s)=0.5d0*abs(wprim(jxc^s,mom(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom(idim))+csound(ixc^s,idim))
1303 
1304  do id=1,ndim
1305  if(id==idim) cycle
1306  ixamax^d=ixcmax^d+kr(id,^d);
1307  ixamin^d=ixcmin^d+kr(id,^d);
1308  hspeed(ixc^s)=max(hspeed(ixc^s),0.5d0*abs(wprim(ixa^s,mom(id))+csound(ixa^s,id)-wprim(ixc^s,mom(id))+csound(ixc^s,id)))
1309  ixamax^d=ixcmax^d-kr(id,^d);
1310  ixamin^d=ixcmin^d-kr(id,^d);
1311  hspeed(ixc^s)=max(hspeed(ixc^s),0.5d0*abs(wprim(ixc^s,mom(id))+csound(ixc^s,id)-wprim(ixa^s,mom(id))+csound(ixa^s,id)))
1312  end do
1313 
1314  do id=1,ndim
1315  if(id==idim) cycle
1316  ixamax^d=jxcmax^d+kr(id,^d);
1317  ixamin^d=jxcmin^d+kr(id,^d);
1318  hspeed(ixc^s)=max(hspeed(ixc^s),0.5d0*abs(wprim(ixa^s,mom(id))+csound(ixa^s,id)-wprim(jxc^s,mom(id))+csound(jxc^s,id)))
1319  ixamax^d=jxcmax^d-kr(id,^d);
1320  ixamin^d=jxcmin^d-kr(id,^d);
1321  hspeed(ixc^s)=max(hspeed(ixc^s),0.5d0*abs(wprim(jxc^s,mom(id))+csound(jxc^s,id)-wprim(ixa^s,mom(id))+csound(ixa^s,id)))
1322  end do
1323 
1324  end subroutine mhd_get_h_speed
1325 
1326  !> Estimating bounds for the minimum and maximum signal velocities
1327  subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
1329 
1330  integer, intent(in) :: ixI^L, ixO^L, idim
1331  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
1332  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
1333  double precision, intent(in) :: x(ixI^S,1:ndim)
1334  double precision, intent(in) :: Hspeed(ixI^S)
1335  double precision, intent(inout) :: cmax(ixI^S)
1336  double precision, intent(inout), optional :: cmin(ixI^S)
1337 
1338  double precision :: wmean(ixI^S,nw)
1339  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
1340  integer :: ix^D
1341 
1342  select case (boundspeed)
1343  case (1)
1344  ! This implements formula (10.52) from "Riemann Solvers and Numerical
1345  ! Methods for Fluid Dynamics" by Toro.
1346  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
1347  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
1348  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
1349  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
1350  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
1351  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
1352  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
1353  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
1354  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
1355  dmean(ixo^s)=sqrt(dmean(ixo^s))
1356  if(present(cmin)) then
1357  cmin(ixo^s)=umean(ixo^s)-dmean(ixo^s)
1358  cmax(ixo^s)=umean(ixo^s)+dmean(ixo^s)
1359  if(h_correction) then
1360  {do ix^db=ixomin^db,ixomax^db\}
1361  cmin(ix^d)=sign(one,cmin(ix^d))*max(abs(cmin(ix^d)),hspeed(ix^d))
1362  cmax(ix^d)=sign(one,cmax(ix^d))*max(abs(cmax(ix^d)),hspeed(ix^d))
1363  {end do\}
1364  end if
1365  else
1366  cmax(ixo^s)=abs(umean(ixo^s))+dmean(ixo^s)
1367  end if
1368  case (2)
1369  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1370  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
1371  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
1372  if(present(cmin)) then
1373  cmax(ixo^s)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
1374  cmin(ixo^s)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
1375  if(h_correction) then
1376  {do ix^db=ixomin^db,ixomax^db\}
1377  cmin(ix^d)=sign(one,cmin(ix^d))*max(abs(cmin(ix^d)),hspeed(ix^d))
1378  cmax(ix^d)=sign(one,cmax(ix^d))*max(abs(cmax(ix^d)),hspeed(ix^d))
1379  {end do\}
1380  end if
1381  else
1382  cmax(ixo^s)=abs(tmp1(ixo^s))+csoundr(ixo^s)
1383  end if
1384  case (3)
1385  ! Miyoshi 2005 JCP 208, 315 equation (67)
1386  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
1387  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
1388  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
1389  if(present(cmin)) then
1390  cmin(ixo^s)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
1391  cmax(ixo^s)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
1392  if(h_correction) then
1393  {do ix^db=ixomin^db,ixomax^db\}
1394  cmin(ix^d)=sign(one,cmin(ix^d))*max(abs(cmin(ix^d)),hspeed(ix^d))
1395  cmax(ix^d)=sign(one,cmax(ix^d))*max(abs(cmax(ix^d)),hspeed(ix^d))
1396  {end do\}
1397  end if
1398  else
1399  cmax(ixo^s)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
1400  end if
1401  end select
1402 
1403  end subroutine mhd_get_cbounds
1404 
1405  !> prepare velocities for ct methods
1406  subroutine mhd_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
1408 
1409  integer, intent(in) :: ixI^L, ixO^L, idim
1410  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
1411  double precision, intent(in) :: cmax(ixI^S)
1412  double precision, intent(in), optional :: cmin(ixI^S)
1413  type(ct_velocity), intent(inout):: vcts
1414 
1415  integer :: idimE,idimN
1416 
1417  ! calculate velocities related to different UCT schemes
1418  select case(type_ct)
1419  case('average')
1420  case('uct_contact')
1421  if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
1422  ! get average normal velocity at cell faces
1423  vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom(idim))+wrp(ixo^s,mom(idim)))
1424  case('uct_hll')
1425  if(.not.allocated(vcts%vbarC)) then
1426  allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
1427  allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
1428  end if
1429  ! Store magnitude of characteristics
1430  if(present(cmin)) then
1431  vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
1432  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
1433  else
1434  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
1435  vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
1436  end if
1437 
1438  idimn=mod(idim,ndir)+1 ! 'Next' direction
1439  idime=mod(idim+1,ndir)+1 ! Electric field direction
1440  ! Store velocities
1441  vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom(idimn))
1442  vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(idimn))
1443  vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
1444  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
1445  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
1446 
1447  vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom(idime))
1448  vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(idime))
1449  vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
1450  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
1451  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
1452  case default
1453  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
1454  end select
1455 
1456  end subroutine mhd_get_ct_velocity
1457 
1458  !> Calculate fast magnetosonic wave speed
1459  subroutine mhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
1461 
1462  integer, intent(in) :: ixI^L, ixO^L, idim
1463  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1464  double precision, intent(out):: csound(ixI^S)
1465  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
1466  double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1467 
1468  inv_rho=1.d0/w(ixo^s,rho_)
1469 
1470  if (mhd_boris_type == boris_reduced_force) then
1471  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma2)
1472  else
1473  gamma2 = 1.0d0
1474  end if
1475 
1476  call mhd_get_csound2(w,x,ixi^l,ixo^l,csound)
1477 
1478  ! store |B|^2 in v
1479  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l) * gamma2
1480 
1481  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1482  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1483  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
1484  * inv_rho * gamma2
1485 
1486  where(avmincs2(ixo^s)<zero)
1487  avmincs2(ixo^s)=zero
1488  end where
1489 
1490  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1491 
1492  if (.not. mhd_hall) then
1493  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1494  if (mhd_boris_type == boris_simplification) then
1495  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
1496  end if
1497  else
1498  ! take the Hall velocity into account:
1499  ! most simple estimate, high k limit:
1500  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
1501  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
1502  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
1503  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
1504  end if
1505 
1506  end subroutine mhd_get_csound
1507 
1508  !> Calculate fast magnetosonic wave speed
1509  subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
1511 
1512  integer, intent(in) :: ixI^L, ixO^L, idim
1513  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1514  double precision, intent(out):: csound(ixI^S)
1515  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
1516  double precision :: inv_rho(ixO^S), gamma_A2(ixO^S)
1517 
1518  inv_rho=1.d0/w(ixo^s,rho_)
1519 
1520  if (mhd_boris_type == boris_reduced_force) then
1521  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a2)
1522  else
1523  gamma_a2 = 1.0d0
1524  end if
1525 
1526  if(mhd_energy) then
1527  csound(ixo^s)=mhd_gamma*w(ixo^s,p_)*inv_rho
1528  else
1529  csound(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
1530  end if
1531  ! store |B|^2 in v
1532  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l) * gamma_a2
1533  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1534  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1535  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
1536  * inv_rho * gamma_a2
1537 
1538  where(avmincs2(ixo^s)<zero)
1539  avmincs2(ixo^s)=zero
1540  end where
1541 
1542  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1543 
1544  if (.not. mhd_hall) then
1545  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1546  if (mhd_boris_type == boris_simplification) then
1547  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
1548  end if
1549  else
1550  ! take the Hall velocity into account:
1551  ! most simple estimate, high k limit:
1552  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
1553  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
1554  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
1555  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
1556  end if
1557 
1558  end subroutine mhd_get_csound_prim
1559 
1560  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
1561  subroutine mhd_get_pthermal(w,x,ixI^L,ixO^L,pth)
1564 
1565  integer, intent(in) :: ixi^l, ixo^l
1566  double precision, intent(in) :: w(ixi^s,nw)
1567  double precision, intent(in) :: x(ixi^s,1:ndim)
1568  double precision, intent(out):: pth(ixi^s)
1569  integer :: iw, ix^d
1570 
1571  if(mhd_energy) then
1572  if(mhd_internal_e) then
1573  pth(ixo^s)=gamma_1*w(ixo^s,e_)
1574  else
1575  pth(ixo^s)=gamma_1*(w(ixo^s,e_)&
1576  - mhd_kin_en(w,ixi^l,ixo^l)&
1577  - mhd_mag_en(w,ixi^l,ixo^l))
1578  end if
1579  else
1580  pth(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
1581  end if
1582 
1583  if (fix_small_values) then
1584  {do ix^db= ixo^lim^db\}
1585  if(pth(ix^d)<small_pressure) then
1586  pth(ix^d)=small_pressure
1587  end if
1588  {enddo^d&\}
1589  end if
1590 
1591  if (check_small_values) then
1592  {do ix^db= ixo^lim^db\}
1593  if(pth(ix^d)<small_pressure) then
1594  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
1595  " encountered when call mhd_get_pthermal"
1596  write(*,*) "Iteration: ", it, " Time: ", global_time
1597  write(*,*) "Location: ", x(ix^d,:)
1598  write(*,*) "Cell number: ", ix^d
1599  do iw=1,nw
1600  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
1601  end do
1602  ! use erroneous arithmetic operation to crash the run
1603  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
1604  write(*,*) "Saving status at the previous time step"
1605  crash=.true.
1606  end if
1607  {enddo^d&\}
1608  end if
1609 
1610  end subroutine mhd_get_pthermal
1611 
1612  !> Calculate temperature=p/rho when in e_ the internal energy is stored
1613  subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
1615  integer, intent(in) :: ixI^L, ixO^L
1616  double precision, intent(in) :: w(ixI^S, 1:nw)
1617  double precision, intent(in) :: x(ixI^S, 1:ndim)
1618  double precision, intent(out):: res(ixI^S)
1619  res(ixo^s) = gamma_1 * w(ixo^s, e_) /w(ixo^s,rho_)
1620  end subroutine mhd_get_temperature_from_eint
1621 
1622  !> Calculate temperature=p/rho when in e_ the total energy is stored
1623  !> this does not check the values of mhd_energy and mhd_internal_e,
1624  !> mhd_energy = .true. and mhd_internal_e = .false.
1625  !> also check small_values is avoided
1626  subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
1628  integer, intent(in) :: ixI^L, ixO^L
1629  double precision, intent(in) :: w(ixI^S, 1:nw)
1630  double precision, intent(in) :: x(ixI^S, 1:ndim)
1631  double precision, intent(out):: res(ixI^S)
1632  res(ixo^s)=(gamma_1*(w(ixo^s,e_)&
1633  - mhd_kin_en(w,ixi^l,ixo^l)&
1634  - mhd_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,rho_)
1635  end subroutine mhd_get_temperature_from_etot
1636 
1637  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
1638  !> csound2=gamma*p/rho
1639  subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
1641  integer, intent(in) :: ixi^l, ixo^l
1642  double precision, intent(in) :: w(ixi^s,nw)
1643  double precision, intent(in) :: x(ixi^s,1:ndim)
1644  double precision, intent(out) :: csound2(ixi^s)
1645 
1646  if(mhd_energy) then
1647  call mhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
1648  csound2(ixo^s)=mhd_gamma*csound2(ixo^s)/w(ixo^s,rho_)
1649  else
1650  csound2(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
1651  end if
1652  end subroutine mhd_get_csound2
1653 
1654  !> Calculate total pressure within ixO^L including magnetic pressure
1655  subroutine mhd_get_p_total(w,x,ixI^L,ixO^L,p)
1657 
1658  integer, intent(in) :: ixI^L, ixO^L
1659  double precision, intent(in) :: w(ixI^S,nw)
1660  double precision, intent(in) :: x(ixI^S,1:ndim)
1661  double precision, intent(out) :: p(ixI^S)
1662 
1663  call mhd_get_pthermal(w,x,ixi^l,ixo^l,p)
1664 
1665  p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1666 
1667  end subroutine mhd_get_p_total
1668 
1669  !> Calculate fluxes within ixO^L.
1670  subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
1672  use mod_geometry
1673 
1674  integer, intent(in) :: ixI^L, ixO^L, idim
1675  ! conservative w
1676  double precision, intent(in) :: wC(ixI^S,nw)
1677  ! primitive w
1678  double precision, intent(in) :: w(ixI^S,nw)
1679  double precision, intent(in) :: x(ixI^S,1:ndim)
1680  double precision,intent(out) :: f(ixI^S,nwflux)
1681 
1682  double precision :: pgas(ixO^S), ptotal(ixO^S),tmp(ixI^S)
1683  double precision :: vHall(ixI^S,1:ndir)
1684  integer :: idirmin, iw, idir, jdir, kdir
1685  double precision, allocatable, dimension(:^D&,:) :: Jambi, btot
1686  double precision, allocatable, dimension(:^D&) :: tmp2, tmp3
1687  double precision :: tmp4(ixO^S)
1688 
1689 
1690  if (mhd_hall) then
1691  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
1692  end if
1693 
1694  if(b0field) tmp(ixo^s)=sum(block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
1695 
1696  if(mhd_energy) then
1697  pgas=w(ixo^s,p_)
1698  else
1699  pgas=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
1700  end if
1701 
1702  ptotal = pgas + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1703 
1704  ! Get flux of density
1705  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
1706 
1707  ! Get flux of tracer
1708  do iw=1,mhd_n_tracer
1709  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
1710  end do
1711 
1712  ! Get flux of momentum
1713  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
1714  if (mhd_boris_type == boris_reduced_force) then
1715  do idir=1,ndir
1716  if(idim==idir) then
1717  f(ixo^s,mom(idir)) = pgas(ixo^s)
1718  else
1719  f(ixo^s,mom(idir)) = 0.0d0
1720  end if
1721  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+w(ixo^s,mom(idim))*wc(ixo^s,mom(idir))
1722  end do
1723  else
1724  ! Normal case (no Boris approximation)
1725  do idir=1,ndir
1726  if(idim==idir) then
1727  f(ixo^s,mom(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
1728  if(b0field) f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+tmp(ixo^s)
1729  else
1730  f(ixo^s,mom(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
1731  end if
1732  if (b0field) then
1733  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))&
1734  -w(ixo^s,mag(idir))*block%B0(ixo^s,idim,idim)&
1735  -w(ixo^s,mag(idim))*block%B0(ixo^s,idir,idim)
1736  end if
1737  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+w(ixo^s,mom(idim))*wc(ixo^s,mom(idir))
1738  end do
1739  end if
1740 
1741  ! Get flux of energy
1742  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
1743  if(mhd_energy) then
1744  if (mhd_internal_e) then
1745  f(ixo^s,e_)=w(ixo^s,mom(idim))*w(ixo^s,p_)
1746  if (mhd_hall) then
1747  call mpistop("solve internal energy not implemented for Hall MHD")
1748  endif
1749  else
1750  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+ptotal(ixo^s))&
1751  -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
1752  if(mhd_solve_eaux) f(ixo^s,eaux_)=w(ixo^s,mom(idim))*wc(ixo^s,eaux_)
1753 
1754  if (b0field) then
1755  f(ixo^s,e_) = f(ixo^s,e_) &
1756  + w(ixo^s,mom(idim)) * tmp(ixo^s) &
1757  - sum(w(ixo^s,mom(:))*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
1758  end if
1759 
1760  if (mhd_hall) then
1761  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
1762  if (mhd_etah>zero) then
1763  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
1764  sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
1765  - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
1766  if (b0field) then
1767  f(ixo^s,e_) = f(ixo^s,e_) &
1768  + vhall(ixo^s,idim) * tmp(ixo^s) &
1769  - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
1770  end if
1771  end if
1772  end if
1773  end if
1774  end if
1775 
1776  ! compute flux of magnetic field
1777  ! f_i[b_k]=v_i*b_k-v_k*b_i
1778  do idir=1,ndir
1779  if (idim==idir) then
1780  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
1781  if (mhd_glm) then
1782  f(ixo^s,mag(idir))=w(ixo^s,psi_)
1783  else
1784  f(ixo^s,mag(idir))=zero
1785  end if
1786  else
1787  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))
1788 
1789  if (b0field) then
1790  f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
1791  +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
1792  -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
1793  end if
1794 
1795  if (mhd_hall) then
1796  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
1797  if (mhd_etah>zero) then
1798  if (b0field) then
1799  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1800  - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+block%B0(ixo^s,idim,idim)) &
1801  + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+block%B0(ixo^s,idir,idim))
1802  else
1803  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1804  - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
1805  + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
1806  end if
1807  end if
1808  end if
1809 
1810  end if
1811  end do
1812 
1813  if (mhd_glm) then
1814  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
1815  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
1816  end if
1817 
1818  ! Contributions of ambipolar term in explicit scheme
1819  if(mhd_ambipolar_exp.and. .not.stagger_grid) then
1820  ! ambipolar electric field
1821  ! E_ambi=-eta_ambi*JxBxB=-JaxBxB=B^2*Ja-(Ja dot B)*B
1822  !Ja=eta_ambi*J=J * mhd_eta_ambi/rho**2
1823  allocate(jambi(ixi^s,1:3))
1824  call mhd_get_jambi(w,x,ixi^l,ixo^l,jambi)
1825  allocate(btot(ixo^s,1:3))
1826  if(b0field) then
1827  do idir=1,3
1828  btot(ixo^s, idir) = w(ixo^s,mag(idir)) + block%B0(ixo^s,idir,idim)
1829  enddo
1830  else
1831  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
1832  endif
1833  allocate(tmp2(ixo^s),tmp3(ixo^s))
1834  !tmp2 = Btot^2
1835  tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
1836  !tmp3 = J_ambi dot Btot
1837  tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
1838 
1839  select case(idim)
1840  case(1)
1841  tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
1842  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(2)) * btot(ixo^s,3) - w(ixo^s,mag(3)) * btot(ixo^s,2)
1843  f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
1844  f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
1845  case(2)
1846  tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
1847  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(3)) * btot(ixo^s,1) - w(ixo^s,mag(1)) * btot(ixo^s,3)
1848  f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
1849  f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
1850  case(3)
1851  tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
1852  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(1)) * btot(ixo^s,2) - w(ixo^s,mag(2)) * btot(ixo^s,1)
1853  f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
1854  f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
1855  endselect
1856 
1857  if(mhd_energy .and. .not. mhd_internal_e) then
1858  f(ixo^s,e_) = f(ixo^s,e_) + tmp2(ixo^s) * tmp(ixo^s)
1859  if(b0field) f(ixo^s,e_) = f(ixo^s,e_) + tmp3(ixo^s) * tmp4(ixo^s)
1860  endif
1861 
1862  deallocate(jambi,btot,tmp2,tmp3)
1863  endif
1864 
1865  end subroutine mhd_get_flux
1866 
1867  !> Source terms J.E in internal energy.
1868  !> For the ambipolar term E = ambiCoef * JxBxB=ambiCoef * B^2(-J_perpB)
1869  !=> the source term J.E = ambiCoef * B^2 * J_perpB^2 = ambiCoef * JxBxB^2/B^2
1870  !> ambiCoef is calculated as mhd_ambi_coef/rho^2, see also the subroutine mhd_get_Jambi
1871  subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
1873  integer, intent(in) :: ixI^L, ixO^L,ie
1874  double precision, intent(in) :: qdt
1875  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
1876  double precision, intent(inout) :: w(ixI^S,1:nw)
1877  double precision :: tmp(ixI^S)
1878  double precision :: jxbxb(ixI^S,1:3)
1879 
1880  call mhd_get_jxbxb(wct,x,ixi^l,ixo^l,jxbxb)
1881  tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=ndim+1) / mhd_mag_en_all(wct, ixi^l, ixo^l)
1882  call multiplyambicoef(ixi^l,ixo^l,tmp,wct,x)
1883  w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
1884 
1886 
1887  subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
1889 
1890  integer, intent(in) :: ixI^L, ixO^L
1891  double precision, intent(in) :: w(ixI^S,nw)
1892  double precision, intent(in) :: x(ixI^S,1:ndim)
1893  double precision, intent(out) :: res(:^D&,:)
1894 
1895  double precision :: btot(ixI^S,1:3)
1896  integer :: idir, idirmin
1897  double precision :: current(ixI^S,7-2*ndir:3)
1898  double precision :: tmp(ixI^S),b2(ixI^S)
1899 
1900  res=0.d0
1901  ! Calculate current density and idirmin
1902  call get_current(w,ixi^l,ixo^l,idirmin,current)
1903  !!!here we know that current has nonzero values only for components in the range idirmin, 3
1904 
1905  if(b0field) then
1906  do idir=1,3
1907  btot(ixo^s, idir) = w(ixo^s,mag(idir)) + block%B0(ixo^s,idir,b0i)
1908  enddo
1909  else
1910  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
1911  endif
1912  tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=ndim+1) !J.B
1913  b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1) !B^2
1914  do idir=1,idirmin-1
1915  res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
1916  enddo
1917  do idir=idirmin,3
1918  res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
1919  enddo
1920  end subroutine mhd_get_jxbxb
1921 
1922  !> Sets the sources for the ambipolar
1923  !> this is used for the STS method
1924  ! The sources are added directly (instead of fluxes as in the explicit)
1925  !> at the corresponding indices
1926  !> store_flux_var is explicitly called for each of the fluxes one by one
1927  subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
1929  use mod_fix_conserve
1930 
1931  integer, intent(in) :: ixI^L, ixO^L,igrid,nflux
1932  double precision, intent(in) :: x(ixI^S,1:ndim)
1933  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
1934  double precision, intent(in) :: my_dt
1935  logical, intent(in) :: fix_conserve_at_step
1936 
1937  double precision, dimension(ixI^S,1:3) :: tmp,ff
1938  double precision :: fluxall(ixI^S,1:nflux,1:ndim)
1939  double precision :: fE(ixI^S,7-2*ndim:3)
1940  double precision :: btot(ixI^S,1:3),tmp2(ixI^S)
1941  integer :: i, ixA^L, ie_
1942 
1943  ixa^l=ixo^l^ladd1;
1944 
1945  fluxall=zero
1946 
1947  call mhd_get_jxbxb(w,x,ixi^l,ixa^l,tmp)
1948 
1949  !set electric field in tmp: E=nuA * jxbxb, where nuA=-etaA/rho^2
1950  do i=1,3
1951  !tmp(ixA^S,i) = -(mhd_eta_ambi/w(ixA^S, rho_)**2) * tmp(ixA^S,i)
1952  call multiplyambicoef(ixi^l,ixa^l,tmp(ixi^s,i),w,x)
1953  enddo
1954 
1955  if(mhd_energy .and. .not.mhd_internal_e) then
1956  !btot should be only mag. pert.
1957  btot(ixa^s,1:3)=0.d0
1958  !if(B0field) then
1959  ! do i=1,ndir
1960  ! btot(ixA^S, i) = w(ixA^S,mag(i)) + block%B0(ixA^S,i,0)
1961  ! enddo
1962  !else
1963  btot(ixa^s,1:ndir) = w(ixa^s,mag(1:ndir))
1964  !endif
1965  call cross_product(ixi^l,ixa^l,tmp,btot,ff)
1966  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
1967  if(fix_conserve_at_step) fluxall(ixi^s,1,1:ndim)=ff(ixi^s,1:ndim)
1968  !- sign comes from the fact that the flux divergence is a source now
1969  wres(ixo^s,e_)=-tmp2(ixo^s)
1970  endif
1971 
1972  if(stagger_grid) then
1973  if(ndir>ndim) then
1974  !!!Bz
1975  ff(ixa^s,1) = tmp(ixa^s,2)
1976  ff(ixa^s,2) = -tmp(ixa^s,1)
1977  ff(ixa^s,3) = 0.d0
1978  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
1979  if(fix_conserve_at_step) fluxall(ixi^s,1+ndir,1:ndim)=ff(ixi^s,1:ndim)
1980  wres(ixo^s,mag(ndir))=-tmp2(ixo^s)
1981  end if
1982  fe=0.d0
1983  call update_faces_ambipolar(ixi^l,ixo^l,w,x,tmp,fe,btot)
1984  ixamax^d=ixomax^d;
1985  ixamin^d=ixomin^d-1;
1986  wres(ixa^s,mag(1:ndim))=-btot(ixa^s,1:ndim)
1987  else
1988  !write curl(ele) as the divergence
1989  !m1={0,ele[[3]],-ele[[2]]}
1990  !m2={-ele[[3]],0,ele[[1]]}
1991  !m3={ele[[2]],-ele[[1]],0}
1992 
1993  !!!Bx
1994  ff(ixa^s,1) = 0.d0
1995  ff(ixa^s,2) = tmp(ixa^s,3)
1996  ff(ixa^s,3) = -tmp(ixa^s,2)
1997  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
1998  if(fix_conserve_at_step) fluxall(ixi^s,2,1:ndim)=ff(ixi^s,1:ndim)
1999  !flux divergence is a source now
2000  wres(ixo^s,mag(1))=-tmp2(ixo^s)
2001  !!!By
2002  ff(ixa^s,1) = -tmp(ixa^s,3)
2003  ff(ixa^s,2) = 0.d0
2004  ff(ixa^s,3) = tmp(ixa^s,1)
2005  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
2006  if(fix_conserve_at_step) fluxall(ixi^s,3,1:ndim)=ff(ixi^s,1:ndim)
2007  wres(ixo^s,mag(2))=-tmp2(ixo^s)
2008 
2009  if(ndir==3) then
2010  !!!Bz
2011  ff(ixa^s,1) = tmp(ixa^s,2)
2012  ff(ixa^s,2) = -tmp(ixa^s,1)
2013  ff(ixa^s,3) = 0.d0
2014  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
2015  if(fix_conserve_at_step) fluxall(ixi^s,1+ndir,1:ndim)=ff(ixi^s,1:ndim)
2016  wres(ixo^s,mag(ndir))=-tmp2(ixo^s)
2017  end if
2018 
2019  end if
2020 
2021  if(fix_conserve_at_step) then
2022  fluxall=my_dt*fluxall
2023  call store_flux(igrid,fluxall,1,ndim,nflux)
2024  if(stagger_grid) then
2025  call store_edge(igrid,ixi^l,my_dt*fe,1,ndim)
2026  end if
2027  end if
2028 
2029  end subroutine sts_set_source_ambipolar
2030 
2031  !> get ambipolar electric field and the integrals around cell faces
2032  subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
2034 
2035  integer, intent(in) :: ixI^L, ixO^L
2036  double precision, intent(in) :: w(ixI^S,1:nw)
2037  double precision, intent(in) :: x(ixI^S,1:ndim)
2038  ! amibipolar electric field at cell centers
2039  double precision, intent(in) :: ECC(ixI^S,1:3)
2040  double precision, intent(out) :: fE(ixI^S,7-2*ndim:3)
2041  double precision, intent(out) :: circ(ixI^S,1:ndim)
2042 
2043  integer :: hxC^L,ixC^L,ixA^L
2044  integer :: idim1,idim2,idir,ix^D
2045 
2046  fe=zero
2047  ! calcuate ambipolar electric field on cell edges from cell centers
2048  do idir=7-2*ndim,3
2049  ixcmax^d=ixomax^d;
2050  ixcmin^d=ixomin^d+kr(idir,^d)-1;
2051  {do ix^db=0,1\}
2052  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
2053  ixamin^d=ixcmin^d+ix^d;
2054  ixamax^d=ixcmax^d+ix^d;
2055  fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
2056  {end do\}
2057  fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
2058  end do
2059 
2060  ! Calculate circulation on each face to get value of line integral of
2061  ! electric field in the positive idir direction.
2062  ixcmax^d=ixomax^d;
2063  ixcmin^d=ixomin^d-1;
2064 
2065  circ=zero
2066 
2067  do idim1=1,ndim ! Coordinate perpendicular to face
2068  do idim2=1,ndim
2069  do idir=7-2*ndim,3 ! Direction of line integral
2070  ! Assemble indices
2071  hxc^l=ixc^l-kr(idim2,^d);
2072  ! Add line integrals in direction idir
2073  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2074  +lvc(idim1,idim2,idir)&
2075  *(fe(ixc^s,idir)&
2076  -fe(hxc^s,idir))
2077  end do
2078  end do
2079  circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
2080  end do
2081 
2082  end subroutine update_faces_ambipolar
2083 
2084  !> use cell-center flux to get cell-face flux
2085  !> and get the source term as the divergence of the flux
2086  subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
2088 
2089  integer, intent(in) :: ixI^L, ixO^L
2090  double precision, dimension(:^D&,:), intent(inout) :: ff
2091  double precision, intent(out) :: src(ixI^S)
2092 
2093  double precision :: ffc(ixI^S,1:ndim)
2094  double precision :: dxinv(ndim)
2095  integer :: idims, ix^D, ixA^L, ixB^L, ixC^L
2096 
2097  ixa^l=ixo^l^ladd1;
2098  dxinv=1.d0/dxlevel
2099  ! cell corner flux in ffc
2100  ffc=0.d0
2101  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
2102  {do ix^db=0,1\}
2103  ixbmin^d=ixcmin^d+ix^d;
2104  ixbmax^d=ixcmax^d+ix^d;
2105  ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
2106  {end do\}
2107  ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
2108  ! flux at cell face
2109  ff(ixi^s,1:ndim)=0.d0
2110  do idims=1,ndim
2111  ixb^l=ixo^l-kr(idims,^d);
2112  ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
2113  {do ix^db=0,1 \}
2114  if({ ix^d==0 .and. ^d==idims | .or.}) then
2115  ixbmin^d=ixcmin^d-ix^d;
2116  ixbmax^d=ixcmax^d-ix^d;
2117  ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
2118  end if
2119  {end do\}
2120  ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
2121  end do
2122  src=0.d0
2123  if(slab_uniform) then
2124  do idims=1,ndim
2125  ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
2126  ixb^l=ixo^l-kr(idims,^d);
2127  src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
2128  end do
2129  else
2130  do idims=1,ndim
2131  ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
2132  ixb^l=ixo^l-kr(idims,^d);
2133  src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
2134  end do
2135  src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
2136  end if
2137  end subroutine get_flux_on_cell_face
2138 
2139  !> Calculates the explicit dt for the ambipokar term
2140  !> This function is used by both explicit scheme and STS method
2141  function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
2143 
2144  integer, intent(in) :: ixi^l, ixo^l
2145  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
2146  double precision, intent(in) :: w(ixi^s,1:nw)
2147  double precision :: dtnew
2148 
2149  double precision :: coef
2150  double precision :: dxarr(ndim)
2151  double precision :: tmp(ixi^s)
2152 
2153  ^d&dxarr(^d)=dx^d;
2154  tmp(ixo^s) = mhd_mag_en_all(w, ixi^l, ixo^l)
2155  call multiplyambicoef(ixi^l,ixo^l,tmp,w,x)
2156  coef = maxval(abs(tmp(ixo^s)))
2157  if(coef/=0.d0) then
2158  coef=1.d0/coef
2159  else
2160  coef=bigdouble
2161  end if
2162  if(slab_uniform) then
2163  dtnew=minval(dxarr(1:ndim))**2.0d0*coef
2164  else
2165  dtnew=minval(block%ds(ixo^s,1:ndim))**2.0d0*coef
2166  end if
2167 
2168  end function get_ambipolar_dt
2169 
2170  !> multiply res by the ambipolar coefficient
2171  !> The ambipolar coefficient is calculated as -mhd_eta_ambi/rho^2
2172  !> The user may mask its value in the user file
2173  !> by implemneting usr_mask_ambipolar subroutine
2174  subroutine multiplyambicoef(ixI^L,ixO^L,res,w,x)
2176  integer, intent(in) :: ixi^l, ixo^l
2177  double precision, intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:ndim)
2178  double precision, intent(inout) :: res(ixi^s)
2179  double precision :: tmp(ixi^s)
2180 
2181  tmp=0.d0
2182  tmp(ixo^s)=-mhd_eta_ambi/w(ixo^s, rho_)**2
2183  if (associated(usr_mask_ambipolar)) then
2184  call usr_mask_ambipolar(ixi^l,ixo^l,w,x,tmp)
2185  end if
2186 
2187  res(ixo^s) = tmp(ixo^s) * res(ixo^s)
2188  end subroutine multiplyambicoef
2189 
2190  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
2191  subroutine mhd_add_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
2195  use mod_gravity, only: gravity_add_source
2196 
2197  integer, intent(in) :: ixI^L, ixO^L
2198  double precision, intent(in) :: qdt
2199  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2200  double precision, intent(inout) :: w(ixI^S,1:nw)
2201  logical, intent(in) :: qsourcesplit
2202  logical, intent(inout) :: active
2203 
2204  if (.not. qsourcesplit) then
2205  if(mhd_internal_e) then
2206  ! Source for solving internal energy
2207  active = .true.
2208  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x,e_)
2209  else if(mhd_solve_eaux) then
2210  ! Source for auxiliary internal energy equation
2211  active = .true.
2212  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x,eaux_)
2213  endif
2214 
2215  ! Source for B0 splitting
2216  if (b0field) then
2217  active = .true.
2218  call add_source_b0split(qdt,ixi^l,ixo^l,wct,w,x)
2219  end if
2220 
2221  ! Sources for resistivity in eqs. for e, B1, B2 and B3
2222  if (abs(mhd_eta)>smalldouble)then
2223  active = .true.
2224  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
2225  end if
2226 
2227  if (mhd_eta_hyper>0.d0)then
2228  active = .true.
2229  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
2230  end if
2231  end if
2232 
2233  {^nooned
2234  if(.not.source_split_divb .and. .not.qsourcesplit .and. istep==nstep) then
2235  ! Sources related to div B
2236  select case (type_divb)
2237  case (divb_none)
2238  ! Do nothing
2239  case (divb_glm)
2240  active = .true.
2241  call add_source_glm(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2242  case (divb_powel)
2243  active = .true.
2244  call add_source_powel(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2245  case (divb_janhunen)
2246  active = .true.
2247  call add_source_janhunen(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2248  case (divb_linde)
2249  active = .true.
2250  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2251  case (divb_lindejanhunen)
2252  active = .true.
2253  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2254  call add_source_janhunen(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2255  case (divb_lindepowel)
2256  active = .true.
2257  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2258  call add_source_powel(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2259  case (divb_lindeglm)
2260  active = .true.
2261  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2262  call add_source_glm(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
2263  case (divb_ct)
2264  continue ! Do nothing
2265  case (divb_multigrid)
2266  continue ! Do nothing
2267  case default
2268  call mpistop('Unknown divB fix')
2269  end select
2270  else if(source_split_divb .and. qsourcesplit) then
2271  ! Sources related to div B
2272  select case (type_divb)
2273  case (divb_none)
2274  ! Do nothing
2275  case (divb_glm)
2276  active = .true.
2277  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
2278  case (divb_powel)
2279  active = .true.
2280  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
2281  case (divb_janhunen)
2282  active = .true.
2283  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
2284  case (divb_linde)
2285  active = .true.
2286  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
2287  case (divb_lindejanhunen)
2288  active = .true.
2289  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
2290  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
2291  case (divb_lindepowel)
2292  active = .true.
2293  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
2294  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
2295  case (divb_lindeglm)
2296  active = .true.
2297  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
2298  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
2299  case (divb_ct)
2300  continue ! Do nothing
2301  case (divb_multigrid)
2302  continue ! Do nothing
2303  case default
2304  call mpistop('Unknown divB fix')
2305  end select
2306  end if
2307  }
2308 
2309  if(mhd_radiative_cooling) then
2310  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,&
2311  w,x,qsourcesplit,active)
2312  end if
2313 
2314  if(mhd_viscosity) then
2315  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
2316  w,x,mhd_energy,qsourcesplit,active)
2317  end if
2318 
2319  if(mhd_gravity) then
2320  call gravity_add_source(qdt,ixi^l,ixo^l,wct,&
2321  w,x,total_energy,qsourcesplit,active)
2322  end if
2323 
2324  if (mhd_boris_type == boris_reduced_force) then
2325  call boris_add_source(qdt,ixi^l,ixo^l,wct,&
2326  w,x,qsourcesplit,active)
2327  end if
2328 
2329  end subroutine mhd_add_source
2330 
2331  subroutine boris_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
2332  qsourcesplit,active)
2334 
2335  integer, intent(in) :: ixI^L, ixO^L
2336  double precision, intent(in) :: qdt, x(ixI^S,1:ndim)
2337  double precision, intent(in) :: wCT(ixI^S,1:nw)
2338  double precision, intent(inout) :: w(ixI^S,1:nw)
2339  logical, intent(in) :: qsourcesplit
2340  logical, intent(inout) :: active
2341 
2342  double precision :: JxB(ixI^S,3)
2343  double precision :: gamma_A2(ixO^S)
2344  integer :: idir
2345 
2346  ! Boris source term is always unsplit
2347  if (qsourcesplit) return
2348 
2349  call get_lorentz(ixi^l,ixo^l,w,jxb)
2350  call mhd_gamma2_alfven(ixi^l, ixo^l, wct, gamma_a2)
2351 
2352  do idir = 1, ndir
2353  w(ixo^s,mom(idir)) = w(ixo^s,mom(idir)) + &
2354  qdt * gamma_a2 * jxb(ixo^s, idir)
2355  end do
2356 
2357  end subroutine boris_add_source
2358 
2359  !> Compute the Lorentz force (JxB)
2360  subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
2362  integer, intent(in) :: ixI^L, ixO^L
2363  double precision, intent(in) :: w(ixI^S,1:nw)
2364  double precision, intent(inout) :: JxB(ixI^S,3)
2365  double precision :: a(ixI^S,3), b(ixI^S,3)
2366  integer :: idir, idirmin
2367  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
2368  double precision :: current(ixI^S,7-2*ndir:3)
2369 
2370  b=0.0d0
2371  do idir = 1, ndir
2372  b(ixo^s, idir) = mhd_mag_i_all(w, ixi^l, ixo^l,idir)
2373  end do
2374 
2375  ! store J current in a
2376  call get_current(w,ixi^l,ixo^l,idirmin,current)
2377 
2378  a=0.0d0
2379  do idir=7-2*ndir,3
2380  a(ixo^s,idir)=current(ixo^s,idir)
2381  end do
2382 
2383  call cross_product(ixi^l,ixo^l,a,b,jxb)
2384  end subroutine get_lorentz
2385 
2386  !> Compute 1/(1+v_A^2/c^2) for Boris' approximation, where v_A is the Alfven
2387  !> velocity
2388  subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
2390  integer, intent(in) :: ixI^L, ixO^L
2391  double precision, intent(in) :: w(ixI^S, nw)
2392  double precision, intent(out) :: gamma_A2(ixO^S)
2393 
2394  if (mhd_boris_c < 0.0d0) then
2395  ! Good for testing the non-conservative momentum treatment
2396  gamma_a2 = 1.0d0
2397  else
2398  ! Compute the inverse of 1 + B^2/(rho * c^2)
2399  gamma_a2 = 1.0d0 / (1.0d0 + mhd_mag_en_all(w, ixi^l, ixo^l) / &
2400  (w(ixo^s, rho_) * mhd_boris_c**2))
2401  end if
2402  end subroutine mhd_gamma2_alfven
2403 
2404  !> Compute 1/sqrt(1+v_A^2/c^2) for Boris simplification, where v_A is the
2405  !> Alfven velocity
2406  function mhd_gamma_alfven(w, ixI^L, ixO^L) result(gamma_A)
2408  integer, intent(in) :: ixi^l, ixo^l
2409  double precision, intent(in) :: w(ixi^s, nw)
2410  double precision :: gamma_a(ixo^s)
2411 
2412  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a)
2413  gamma_a = sqrt(gamma_a)
2414  end function mhd_gamma_alfven
2415 
2416  subroutine internal_energy_add_source(qdt,ixI^L,ixO^L,wCT,w,x,ie)
2418  use mod_geometry
2419 
2420  integer, intent(in) :: ixI^L, ixO^L,ie
2421  double precision, intent(in) :: qdt
2422  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2423  double precision, intent(inout) :: w(ixI^S,1:nw)
2424  double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
2425 
2426  call mhd_get_v(wct,x,ixi^l,ixi^l,v)
2427  if(slab_uniform) then
2428  if(nghostcells .gt. 2) then
2429  call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
2430  else
2431  call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
2432  end if
2433  else
2434  call divvector(v,ixi^l,ixo^l,divv)
2435  end if
2436  call mhd_get_pthermal(wct,x,ixi^l,ixo^l,pth)
2437  w(ixo^s,ie)=w(ixo^s,ie)-qdt*pth(ixo^s)*divv(ixo^s)
2438  if(mhd_ambipolar)then
2439  call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,ie)
2440  end if
2441  if(fix_small_values) then
2442  call mhd_handle_small_ei(w,x,ixi^l,ixo^l,ie,'internal_energy_add_source')
2443  end if
2444  end subroutine internal_energy_add_source
2445 
2446  !> handle small or negative internal energy
2447  subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
2449  use mod_small_values
2450  integer, intent(in) :: ixI^L,ixO^L, ie
2451  double precision, intent(inout) :: w(ixI^S,1:nw)
2452  double precision, intent(in) :: x(ixI^S,1:ndim)
2453  character(len=*), intent(in) :: subname
2454 
2455  integer :: idir
2456  logical :: flag(ixI^S,1:nw)
2457 
2458  flag=.false.
2459  where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
2460  if(any(flag(ixo^s,ie))) then
2461  select case (small_values_method)
2462  case ("replace")
2463  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
2464  case ("average")
2465  call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
2466  case default
2467  ! small values error shows primitive variables
2468  w(ixo^s,e_)=w(ixo^s,e_)*gamma_1
2469  do idir = 1, ndir
2470  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/w(ixo^s,rho_)
2471  end do
2472  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2473  end select
2474  end if
2475 
2476  end subroutine mhd_handle_small_ei
2477 
2478  !> Source terms after split off time-independent magnetic field
2479  subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
2481 
2482  integer, intent(in) :: ixI^L, ixO^L
2483  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2484  double precision, intent(inout) :: w(ixI^S,1:nw)
2485 
2486  double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
2487  integer :: idir
2488 
2489  a=0.d0
2490  b=0.d0
2491  ! for force-free field J0xB0 =0
2492  if(.not.b0field_forcefree) then
2493  ! store B0 magnetic field in b
2494  b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
2495 
2496  ! store J0 current in a
2497  do idir=7-2*ndir,3
2498  a(ixo^s,idir)=block%J0(ixo^s,idir)
2499  end do
2500  call cross_product(ixi^l,ixo^l,a,b,axb)
2501  axb(ixo^s,:)=axb(ixo^s,:)*qdt
2502  ! add J0xB0 source term in momentum equations
2503  w(ixo^s,mom(1:ndir))=w(ixo^s,mom(1:ndir))+axb(ixo^s,1:ndir)
2504  end if
2505 
2506  if(total_energy) then
2507  a=0.d0
2508  ! for free-free field -(vxB0) dot J0 =0
2509  b(ixo^s,:)=wct(ixo^s,mag(:))
2510  ! store full magnetic field B0+B1 in b
2511  if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
2512  ! store velocity in a
2513  do idir=1,ndir
2514  a(ixo^s,idir)=wct(ixo^s,mom(idir))/wct(ixo^s,rho_)
2515  end do
2516  call cross_product(ixi^l,ixo^l,a,b,axb)
2517  axb(ixo^s,:)=axb(ixo^s,:)*qdt
2518  ! add -(vxB) dot J0 source term in energy equation
2519  do idir=7-2*ndir,3
2520  w(ixo^s,e_)=w(ixo^s,e_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
2521  end do
2522  end if
2523 
2524  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
2525 
2526  end subroutine add_source_b0split
2527 
2528  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
2529  !> each direction, non-conservative. If the fourthorder precompiler flag is
2530  !> set, uses fourth order central difference for the laplacian. Then the
2531  !> stencil is 5 (2 neighbours).
2532  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
2534  use mod_usr_methods
2535  use mod_geometry
2536 
2537  integer, intent(in) :: ixI^L, ixO^L
2538  double precision, intent(in) :: qdt
2539  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2540  double precision, intent(inout) :: w(ixI^S,1:nw)
2541  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
2542  integer :: lxO^L, kxO^L
2543 
2544  double precision :: tmp(ixI^S),tmp2(ixI^S)
2545 
2546  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
2547  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
2548  double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
2549 
2550  ! Calculating resistive sources involve one extra layer
2551  if (mhd_4th_order) then
2552  ixa^l=ixo^l^ladd2;
2553  else
2554  ixa^l=ixo^l^ladd1;
2555  end if
2556 
2557  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
2558  call mpistop("Error in add_source_res1: Non-conforming input limits")
2559 
2560  ! Calculate current density and idirmin
2561  call get_current(wct,ixi^l,ixo^l,idirmin,current)
2562 
2563  if (mhd_eta>zero)then
2564  eta(ixa^s)=mhd_eta
2565  gradeta(ixo^s,1:ndim)=zero
2566  else
2567  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
2568  ! assumes that eta is not function of current?
2569  do idim=1,ndim
2570  call gradient(eta,ixi^l,ixo^l,idim,tmp)
2571  gradeta(ixo^s,idim)=tmp(ixo^s)
2572  end do
2573  end if
2574 
2575  if(b0field) then
2576  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
2577  else
2578  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
2579  end if
2580 
2581  do idir=1,ndir
2582  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
2583  if (mhd_4th_order) then
2584  tmp(ixo^s)=zero
2585  tmp2(ixi^s)=bf(ixi^s,idir)
2586  do idim=1,ndim
2587  lxo^l=ixo^l+2*kr(idim,^d);
2588  jxo^l=ixo^l+kr(idim,^d);
2589  hxo^l=ixo^l-kr(idim,^d);
2590  kxo^l=ixo^l-2*kr(idim,^d);
2591  tmp(ixo^s)=tmp(ixo^s)+&
2592  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
2593  /(12.0d0 * dxlevel(idim)**2)
2594  end do
2595  else
2596  tmp(ixo^s)=zero
2597  tmp2(ixi^s)=bf(ixi^s,idir)
2598  do idim=1,ndim
2599  jxo^l=ixo^l+kr(idim,^d);
2600  hxo^l=ixo^l-kr(idim,^d);
2601  tmp(ixo^s)=tmp(ixo^s)+&
2602  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
2603  end do
2604  end if
2605 
2606  ! Multiply by eta
2607  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
2608 
2609  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
2610  if (mhd_eta<zero)then
2611  do jdir=1,ndim; do kdir=idirmin,3
2612  if (lvc(idir,jdir,kdir)/=0)then
2613  if (lvc(idir,jdir,kdir)==1)then
2614  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
2615  else
2616  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
2617  end if
2618  end if
2619  end do; end do
2620  end if
2621 
2622  ! Add sources related to eta*laplB-grad(eta) x J to B and e
2623  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
2624  if (mhd_energy) then
2625  w(ixo^s,e_)=w(ixo^s,e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
2626  end if
2627  end do ! idir
2628 
2629  if(mhd_energy) then
2630  ! de/dt+=eta*J**2
2631  tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
2632  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)
2633  if(mhd_solve_eaux) then
2634  ! add eta*J**2 source term in the internal energy equation
2635  w(ixo^s,eaux_)=w(ixo^s,eaux_)+tmp(ixo^s)
2636  end if
2637  end if
2638 
2639  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
2640 
2641  end subroutine add_source_res1
2642 
2643  !> Add resistive source to w within ixO
2644  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
2645  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
2647  use mod_usr_methods
2648  use mod_geometry
2649 
2650  integer, intent(in) :: ixI^L, ixO^L
2651  double precision, intent(in) :: qdt
2652  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2653  double precision, intent(inout) :: w(ixI^S,1:nw)
2654 
2655  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
2656  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
2657  double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
2658  integer :: ixA^L,idir,idirmin,idirmin1
2659 
2660  ixa^l=ixo^l^ladd2;
2661 
2662  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
2663  call mpistop("Error in add_source_res2: Non-conforming input limits")
2664 
2665  ixa^l=ixo^l^ladd1;
2666  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
2667  ! Determine exact value of idirmin while doing the loop.
2668  call get_current(wct,ixi^l,ixa^l,idirmin,current)
2669 
2670  if (mhd_eta>zero)then
2671  eta(ixa^s)=mhd_eta
2672  else
2673  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
2674  end if
2675 
2676  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
2677  tmpvec(ixa^s,1:ndir)=zero
2678  do idir=idirmin,3
2679  tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
2680  end do
2681  curlj=0.d0
2682  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
2683  if(stagger_grid) then
2684  if(ndim==2.and.ndir==3) then
2685  ! if 2.5D
2686  w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
2687  end if
2688  else
2689  w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
2690  end if
2691 
2692  if(mhd_energy) then
2693  ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
2694  ! de1/dt= eta J^2 - B1 dot curl(eta J)
2695  tmp(ixo^s)=eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
2696  w(ixo^s,e_)=w(ixo^s,e_)+qdt*(tmp(ixo^s)-&
2697  sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
2698  if(mhd_solve_eaux) then
2699  ! add eta*J**2 source term in the internal energy equation
2700  w(ixo^s,eaux_)=w(ixo^s,eaux_)+tmp(ixo^s)
2701  end if
2702  end if
2703 
2704  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
2705  end subroutine add_source_res2
2706 
2707  !> Add Hyper-resistive source to w within ixO
2708  !> Uses 9 point stencil (4 neighbours) in each direction.
2709  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
2711  use mod_geometry
2712 
2713  integer, intent(in) :: ixI^L, ixO^L
2714  double precision, intent(in) :: qdt
2715  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2716  double precision, intent(inout) :: w(ixI^S,1:nw)
2717  !.. local ..
2718  double precision :: current(ixI^S,7-2*ndir:3)
2719  double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
2720  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
2721 
2722  ixa^l=ixo^l^ladd3;
2723  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
2724  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
2725 
2726  call get_current(wct,ixi^l,ixa^l,idirmin,current)
2727  tmpvec(ixa^s,1:ndir)=zero
2728  do jdir=idirmin,3
2729  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
2730  end do
2731 
2732  ixa^l=ixo^l^ladd2;
2733  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
2734 
2735  ixa^l=ixo^l^ladd1;
2736  tmpvec(ixa^s,1:ndir)=zero
2737  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
2738  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mhd_eta_hyper
2739 
2740  ixa^l=ixo^l;
2741  tmpvec2(ixa^s,1:ndir)=zero
2742  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
2743 
2744  do idir=1,ndir
2745  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
2746  end do
2747 
2748  if (mhd_energy) then
2749  ! de/dt= +div(B x Ehyper)
2750  ixa^l=ixo^l^ladd1;
2751  tmpvec2(ixa^s,1:ndir)=zero
2752  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
2753  tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
2754  + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
2755  end do; end do; end do
2756  tmp(ixo^s)=zero
2757  call divvector(tmpvec2,ixi^l,ixo^l,tmp)
2758  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
2759  end if
2760 
2761  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
2762 
2763  end subroutine add_source_hyperres
2764 
2765  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
2766  ! Add divB related sources to w within ixO
2767  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
2768  ! giving the EGLM-MHD scheme
2770  use mod_geometry
2771 
2772  integer, intent(in) :: ixI^L, ixO^L
2773  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2774  double precision, intent(inout) :: w(ixI^S,1:nw)
2775  double precision:: divb(ixI^S)
2776  integer :: idim,idir
2777  double precision :: gradPsi(ixI^S)
2778 
2779  ! We calculate now div B
2780  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
2781 
2782  ! dPsi/dt = - Ch^2/Cp^2 Psi
2783  if (mhd_glm_alpha < zero) then
2784  w(ixo^s,psi_) = abs(mhd_glm_alpha)*wct(ixo^s,psi_)
2785  else
2786  ! implicit update of Psi variable
2787  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
2788  if(slab_uniform) then
2789  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
2790  else
2791  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
2792  end if
2793  end if
2794 
2795  ! gradient of Psi
2796  do idim=1,ndim
2797  select case(typegrad)
2798  case("central")
2799  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
2800  case("limited")
2801  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
2802  end select
2803  if (total_energy) then
2804  ! e = e -qdt (b . grad(Psi))
2805  w(ixo^s,e_) = w(ixo^s,e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
2806  end if
2807  end do
2808 
2809  ! m = m - qdt b div b
2810  do idir=1,ndir
2811  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
2812  end do
2813 
2814  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
2815 
2816  end subroutine add_source_glm
2817 
2818  !> Add divB related sources to w within ixO corresponding to Powel
2819  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
2821 
2822  integer, intent(in) :: ixI^L, ixO^L
2823  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2824  double precision, intent(inout) :: w(ixI^S,1:nw)
2825  double precision :: divb(ixI^S),v(ixI^S,1:ndir)
2826  integer :: idir
2827 
2828  ! We calculate now div B
2829  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
2830 
2831  ! calculate velocity
2832  call mhd_get_v(wct,x,ixi^l,ixo^l,v)
2833 
2834  if (total_energy) then
2835  ! e = e - qdt (v . b) * div b
2836  w(ixo^s,e_)=w(ixo^s,e_)-&
2837  qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
2838  end if
2839 
2840  ! b = b - qdt v * div b
2841  do idir=1,ndir
2842  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
2843  end do
2844 
2845  ! m = m - qdt b div b
2846  do idir=1,ndir
2847  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
2848  end do
2849 
2850  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
2851 
2852  end subroutine add_source_powel
2853 
2854  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
2855  ! Add divB related sources to w within ixO
2856  ! corresponding to Janhunen, just the term in the induction equation.
2858 
2859  integer, intent(in) :: ixI^L, ixO^L
2860  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2861  double precision, intent(inout) :: w(ixI^S,1:nw)
2862  double precision :: divb(ixI^S)
2863  integer :: idir
2864 
2865  ! We calculate now div B
2866  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
2867 
2868  ! b = b - qdt v * div b
2869  do idir=1,ndir
2870  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))/wct(ixo^s,rho_)*divb(ixo^s)
2871  end do
2872 
2873  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
2874 
2875  end subroutine add_source_janhunen
2876 
2877  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
2878  ! Add Linde's divB related sources to wnew within ixO
2880  use mod_geometry
2881 
2882  integer, intent(in) :: ixI^L, ixO^L
2883  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
2884  double precision, intent(inout) :: w(ixI^S,1:nw)
2885  integer :: idim, idir, ixp^L, i^D, iside
2886  double precision :: divb(ixI^S),graddivb(ixI^S)
2887  logical, dimension(-1:1^D&) :: leveljump
2888 
2889  ! Calculate div B
2890  ixp^l=ixo^l^ladd1;
2891  call get_divb(wct,ixi^l,ixp^l,divb, mhd_divb_4thorder)
2892 
2893  ! for AMR stability, retreat one cell layer from the boarders of level jump
2894  {do i^db=-1,1\}
2895  if(i^d==0|.and.) cycle
2896  if(neighbor_type(i^d,block%igrid)==2 .or. neighbor_type(i^d,block%igrid)==4) then
2897  leveljump(i^d)=.true.
2898  else
2899  leveljump(i^d)=.false.
2900  end if
2901  {end do\}
2902 
2903  ixp^l=ixo^l;
2904  do idim=1,ndim
2905  select case(idim)
2906  {case(^d)
2907  do iside=1,2
2908  i^dd=kr(^dd,^d)*(2*iside-3);
2909  if (leveljump(i^dd)) then
2910  if (iside==1) then
2911  ixpmin^d=ixomin^d-i^d
2912  else
2913  ixpmax^d=ixomax^d-i^d
2914  end if
2915  end if
2916  end do
2917  \}
2918  end select
2919  end do
2920 
2921  ! Add Linde's diffusive terms
2922  do idim=1,ndim
2923  ! Calculate grad_idim(divb)
2924  select case(typegrad)
2925  case("central")
2926  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
2927  case("limited")
2928  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
2929  end select
2930 
2931  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
2932  if (slab_uniform) then
2933  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
2934  else
2935  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
2936  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
2937  end if
2938 
2939  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
2940 
2941  if (typedivbdiff=='all' .and. total_energy) then
2942  ! e += B_idim*eta*grad_idim(divb)
2943  w(ixp^s,e_)=w(ixp^s,e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
2944  end if
2945  end do
2946 
2947  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
2948 
2949  end subroutine add_source_linde
2950 
2951  !> Calculate div B within ixO
2952  subroutine get_divb(w,ixI^L,ixO^L,divb, fourthorder)
2953 
2955  use mod_geometry
2956 
2957  integer, intent(in) :: ixi^l, ixo^l
2958  double precision, intent(in) :: w(ixi^s,1:nw)
2959  double precision, intent(inout) :: divb(ixi^s)
2960  logical, intent(in), optional :: fourthorder
2961 
2962  double precision :: bvec(ixi^s,1:ndir)
2963  double precision :: divb_corner(ixi^s), sign
2964  double precision :: aux_vol(ixi^s)
2965  integer :: ixc^l, idir, ic^d, ix^l
2966 
2967  if(stagger_grid) then
2968  divb=0.d0
2969  do idir=1,ndim
2970  ixc^l=ixo^l-kr(idir,^d);
2971  divb(ixo^s)=divb(ixo^s)+block%ws(ixo^s,idir)*block%surfaceC(ixo^s,idir)-&
2972  block%ws(ixc^s,idir)*block%surfaceC(ixc^s,idir)
2973  end do
2974  divb(ixo^s)=divb(ixo^s)/block%dvolume(ixo^s)
2975  else
2976  bvec(ixi^s,:)=w(ixi^s,mag(:))
2977  select case(typediv)
2978  case("central")
2979  call divvector(bvec,ixi^l,ixo^l,divb,fourthorder)
2980  case("limited")
2981  call divvectors(bvec,ixi^l,ixo^l,divb)
2982  end select
2983  end if
2984 
2985  end subroutine get_divb
2986 
2987  !> get dimensionless div B = |divB| * volume / area / |B|
2988  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
2989 
2991 
2992  integer, intent(in) :: ixi^l, ixo^l
2993  double precision, intent(in) :: w(ixi^s,1:nw)
2994  double precision :: divb(ixi^s), dsurface(ixi^s)
2995 
2996  double precision :: invb(ixo^s)
2997  integer :: ixa^l,idims
2998 
2999  call get_divb(w,ixi^l,ixo^l,divb)
3000  invb(ixo^s)=sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))
3001  where(invb(ixo^s)/=0.d0)
3002  invb(ixo^s)=1.d0/invb(ixo^s)
3003  end where
3004  if(slab_uniform) then
3005  divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/dxlevel(:))
3006  else
3007  ixamin^d=ixomin^d-1;
3008  ixamax^d=ixomax^d-1;
3009  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
3010  do idims=1,ndim
3011  ixa^l=ixo^l-kr(idims,^d);
3012  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
3013  end do
3014  divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
3015  block%dvolume(ixo^s)/dsurface(ixo^s)
3016  end if
3017 
3018  end subroutine get_normalized_divb
3019 
3020  !> Calculate idirmin and the idirmin:3 components of the common current array
3021  !> make sure that dxlevel(^D) is set correctly.
3022  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
3024  use mod_geometry
3025 
3026  integer, intent(in) :: ixo^l, ixi^l
3027  double precision, intent(in) :: w(ixi^s,1:nw)
3028  integer, intent(out) :: idirmin
3029  integer :: idir, idirmin0
3030 
3031  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
3032  double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)
3033 
3034  idirmin0 = 7-2*ndir
3035 
3036  bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
3037 
3038  call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
3039 
3040  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
3041  block%J0(ixo^s,idirmin0:3)
3042 
3043  end subroutine get_current
3044 
3045  !> If resistivity is not zero, check diffusion time limit for dt
3046  subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
3048  use mod_usr_methods
3050  use mod_viscosity, only: viscosity_get_dt
3051  use mod_gravity, only: gravity_get_dt
3052 
3053  integer, intent(in) :: ixI^L, ixO^L
3054  double precision, intent(inout) :: dtnew
3055  double precision, intent(in) :: dx^D
3056  double precision, intent(in) :: w(ixI^S,1:nw)
3057  double precision, intent(in) :: x(ixI^S,1:ndim)
3058 
3059  integer :: idirmin,idim
3060  double precision :: dxarr(ndim)
3061  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
3062 
3063  dtnew = bigdouble
3064 
3065  ^d&dxarr(^d)=dx^d;
3066  if (mhd_eta>zero)then
3067  dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mhd_eta
3068  else if (mhd_eta<zero)then
3069  call get_current(w,ixi^l,ixo^l,idirmin,current)
3070  call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
3071  dtnew=bigdouble
3072  do idim=1,ndim
3073  if(slab_uniform) then
3074  dtnew=min(dtnew,&
3075  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
3076  else
3077  dtnew=min(dtnew,&
3078  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
3079  end if
3080  end do
3081  end if
3082 
3083  if(mhd_eta_hyper>zero) then
3084  if(slab_uniform) then
3085  dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mhd_eta_hyper,dtnew)
3086  else
3087  dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mhd_eta_hyper,dtnew)
3088  end if
3089  end if
3090 
3091  if(mhd_radiative_cooling) then
3092  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
3093  end if
3094 
3095  if(mhd_viscosity) then
3096  call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
3097  end if
3098 
3099  if(mhd_gravity) then
3100  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
3101  end if
3102 
3103  if(mhd_ambipolar_exp) then
3104  dtnew=min(dtdiffpar*get_ambipolar_dt(w,ixi^l,ixo^l,dx^d,x),dtnew)
3105  endif
3106 
3107  end subroutine mhd_get_dt
3108 
3109  ! Add geometrical source terms to w
3110  subroutine mhd_add_source_geom(qdt,ixI^L,ixO^L,wCT,w,x)
3112  use mod_geometry
3113 
3114  integer, intent(in) :: ixI^L, ixO^L
3115  double precision, intent(in) :: qdt, x(ixI^S,1:ndim)
3116  double precision, intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
3117 
3118  integer :: iw,idir, h1x^L{^NOONED, h2x^L}
3119  double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S)
3120 
3121  integer :: mr_,mphi_ ! Polar var. names
3122  integer :: br_,bphi_
3123 
3124  mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
3125  br_=mag(1); bphi_=mag(1)-1+phi_
3126 
3127  select case (coordinate)
3128  case (cylindrical)
3129  if (angmomfix) then
3130  call mpistop("angmomfix not implemented yet in MHD")
3131  endif
3132  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp)
3133  if(phi_>0) then
3134  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
3135  wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/wct(ixo^s,rho_))
3136  w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
3137  -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/wct(ixo^s,rho_) &
3138  +wct(ixo^s,bphi_)*wct(ixo^s,br_))
3139  if(.not.stagger_grid) then
3140  w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
3141  (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
3142  -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
3143  /wct(ixo^s,rho_)
3144  end if
3145  else
3146  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
3147  end if
3148  if(mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
3149  case (spherical)
3150  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
3151  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp1)
3152  tmp(ixo^s)=tmp1(ixo^s)
3153  if(b0field) then
3154  tmp2(ixo^s)=sum(block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
3155  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
3156  end if
3157  ! m1
3158  tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
3159  *(block%surfaceC(ixo^s,1)-block%surfaceC(h1x^s,1))/block%dvolume(ixo^s)
3160  if(ndir>1) then
3161  do idir=2,ndir
3162  tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(idir))**2/wct(ixo^s,rho_)-wct(ixo^s,mag(idir))**2
3163  if(b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
3164  end do
3165  end if
3166  w(ixo^s,mom(1))=w(ixo^s,mom(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
3167  ! b1
3168  if(mhd_glm) then
3169  w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,psi_)
3170  end if
3171 
3172  {^nooned
3173  ! m2
3174  tmp(ixo^s)=tmp1(ixo^s)
3175  if(b0field) then
3176  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
3177  end if
3178  ! This will make hydrostatic p=const an exact solution
3179  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s) &
3180  *(block%surfaceC(ixo^s,2)-block%surfaceC(h2x^s,2)) &
3181  /block%dvolume(ixo^s)
3182  tmp(ixo^s)=-(wct(ixo^s,mom(1))*wct(ixo^s,mom(2))/wct(ixo^s,rho_) &
3183  -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
3184  if (b0field) then
3185  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
3186  +wct(ixo^s,mag(1))*block%B0(ixo^s,2,0)
3187  end if
3188  if(ndir==3) then
3189  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(3))**2/wct(ixo^s,rho_) &
3190  -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
3191  if (b0field) then
3192  tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
3193  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
3194  end if
3195  end if
3196  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
3197  ! b2
3198  if(.not.stagger_grid) then
3199  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
3200  -wct(ixo^s,mom(2))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_)
3201  if(b0field) then
3202  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,2,0) &
3203  -wct(ixo^s,mom(2))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_)
3204  end if
3205  if(mhd_glm) then
3206  tmp(ixo^s)=tmp(ixo^s) &
3207  + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,psi_)
3208  end if
3209  w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
3210  end if
3211  }
3212 
3213  if(ndir==3) then
3214  ! m3
3215  if(.not.angmomfix) then
3216  tmp(ixo^s)=-(wct(ixo^s,mom(3))*wct(ixo^s,mom(1))/wct(ixo^s,rho_) &
3217  -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
3218  -(wct(ixo^s,mom(2))*wct(ixo^s,mom(3))/wct(ixo^s,rho_) &
3219  -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
3220  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
3221  if (b0field) then
3222  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
3223  +wct(ixo^s,mag(1))*block%B0(ixo^s,3,0) {^nooned &
3224  +(block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
3225  +wct(ixo^s,mag(2))*block%B0(ixo^s,3,0)) &
3226  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
3227  end if
3228  w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
3229  else
3230  call mpistop("angmomfix not implemented yet in MHD")
3231  end if
3232  ! b3
3233  if(.not.stagger_grid) then
3234  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
3235  -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_) {^nooned &
3236  -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
3237  -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
3238  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
3239  if (b0field) then
3240  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,3,0) &
3241  -wct(ixo^s,mom(3))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_){^nooned &
3242  -(wct(ixo^s,mom(3))*block%B0(ixo^s,2,0) &
3243  -wct(ixo^s,mom(2))*block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
3244  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
3245  end if
3246  w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
3247  end if
3248  end if
3249  end select
3250  end subroutine mhd_add_source_geom
3251 
3252  !> Compute 2 times total magnetic energy
3253  function mhd_mag_en_all(w, ixI^L, ixO^L) result(mge)
3255  integer, intent(in) :: ixi^l, ixo^l
3256  double precision, intent(in) :: w(ixi^s, nw)
3257  double precision :: mge(ixo^s)
3258 
3259  if (b0field) then
3260  mge = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,b0i))**2, dim=ndim+1)
3261  else
3262  mge = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3263  end if
3264  end function mhd_mag_en_all
3265 
3266  !> Compute full magnetic field by direction
3267  function mhd_mag_i_all(w, ixI^L, ixO^L,idir) result(mgf)
3269  integer, intent(in) :: ixi^l, ixo^l, idir
3270  double precision, intent(in) :: w(ixi^s, nw)
3271  double precision :: mgf(ixo^s)
3272 
3273  if (b0field) then
3274  mgf = w(ixo^s, mag(idir))+block%B0(ixo^s,idir,b0i)
3275  else
3276  mgf = w(ixo^s, mag(idir))
3277  end if
3278  end function mhd_mag_i_all
3279 
3280  !> Compute evolving magnetic energy
3281  function mhd_mag_en(w, ixI^L, ixO^L) result(mge)
3282  use mod_global_parameters, only: nw, ndim
3283  integer, intent(in) :: ixi^l, ixo^l
3284  double precision, intent(in) :: w(ixi^s, nw)
3285  double precision :: mge(ixo^s)
3286 
3287  mge = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3288  end function mhd_mag_en
3289 
3290  !> compute kinetic energy
3291  function mhd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
3292  use mod_global_parameters, only: nw, ndim
3293  integer, intent(in) :: ixi^l, ixo^l
3294  double precision, intent(in) :: w(ixi^s, nw)
3295  double precision :: ke(ixo^s)
3296  double precision, intent(in), optional :: inv_rho(ixo^s)
3297 
3298  if (present(inv_rho)) then
3299  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
3300  else
3301  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
3302  end if
3303  end function mhd_kin_en
3304 
3305  subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
3307 
3308  integer, intent(in) :: ixI^L, ixO^L
3309  double precision, intent(in) :: w(ixI^S,nw)
3310  double precision, intent(in) :: x(ixI^S,1:ndim)
3311  double precision, intent(inout) :: vHall(ixI^S,1:3)
3312 
3313  integer :: idir, idirmin
3314  double precision :: current(ixI^S,7-2*ndir:3)
3315 
3316  ! Calculate current density and idirmin
3317  call get_current(w,ixi^l,ixo^l,idirmin,current)
3318  vhall(ixo^s,1:3) = zero
3319  vhall(ixo^s,idirmin:3) = - mhd_etah*current(ixo^s,idirmin:3)
3320  do idir = idirmin, 3
3321  vhall(ixo^s,idir) = vhall(ixo^s,idir)/w(ixo^s,rho_)
3322  end do
3323 
3324  end subroutine mhd_getv_hall
3325 
3326  subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
3328 
3329  integer, intent(in) :: ixI^L, ixO^L
3330  double precision, intent(in) :: w(ixI^S,nw)
3331  double precision, intent(in) :: x(ixI^S,1:ndim)
3332  double precision, allocatable, intent(inout) :: res(:^D&,:)
3333 
3334 
3335  integer :: idir, idirmin
3336  double precision :: current(ixI^S,7-2*ndir:3)
3337 
3338  res = 0d0
3339 
3340  ! Calculate current density and idirmin
3341  call get_current(w,ixi^l,ixo^l,idirmin,current)
3342 
3343  res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
3344  do idir = idirmin, 3
3345  call multiplyambicoef(ixi^l,ixo^l,res(ixi^s,idir),w,x)
3346  enddo
3347 
3348  end subroutine mhd_get_jambi
3349 
3350  subroutine mhd_getdt_hall(w,x,ixI^L,ixO^L,dx^D,dthall)
3352 
3353  integer, intent(in) :: ixI^L, ixO^L
3354  double precision, intent(in) :: dx^D
3355  double precision, intent(in) :: w(ixI^S,1:nw)
3356  double precision, intent(in) :: x(ixI^S,1:ndim)
3357  double precision, intent(out) :: dthall
3358  !.. local ..
3359  double precision :: dxarr(ndim)
3360  double precision :: bmag(ixI^S)
3361 
3362  dthall=bigdouble
3363 
3364  ! because we have that in cmax now:
3365  return
3366 
3367  ^d&dxarr(^d)=dx^d;
3368 
3369  if (.not. b0field) then
3370  bmag(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2, dim=ndim+1))
3371  bmag(ixo^s)=sqrt(sum((w(ixo^s,mag(:)) + block%B0(ixo^s,1:ndir,b0i))**2))
3372  end if
3373 
3374  if(slab_uniform) then
3375  dthall=dtdiffpar*minval(dxarr(1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
3376  else
3377  dthall=dtdiffpar*minval(block%ds(ixo^s,1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
3378  end if
3379 
3380  end subroutine mhd_getdt_hall
3381 
3382  subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
3384  use mod_usr_methods
3385  integer, intent(in) :: ixI^L, ixO^L, idir
3386  double precision, intent(in) :: qt
3387  double precision, intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
3388  double precision, intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
3389  type(state) :: s
3390  double precision :: dB(ixI^S), dPsi(ixI^S)
3391 
3392  if(stagger_grid) then
3393  wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
3394  wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
3395  wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
3396  wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
3397  else
3398  ! Solve the Riemann problem for the linear 2x2 system for normal
3399  ! B-field and GLM_Psi according to Dedner 2002:
3400  ! This implements eq. (42) in Dedner et al. 2002 JcP 175
3401  ! Gives the Riemann solution on the interface
3402  ! for the normal B component and Psi in the GLM-MHD system.
3403  ! 23/04/2013 Oliver Porth
3404  db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
3405  dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
3406 
3407  wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
3408  - 0.5d0/cmax_global * dpsi(ixo^s)
3409  wlp(ixo^s,psi_) = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
3410  - 0.5d0*cmax_global * db(ixo^s)
3411 
3412  wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
3413  wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
3414 
3415  if(total_energy) then
3416  wrc(ixo^s,e_)=wrc(ixo^s,e_)-half*wrc(ixo^s,mag(idir))**2
3417  wlc(ixo^s,e_)=wlc(ixo^s,e_)-half*wlc(ixo^s,mag(idir))**2
3418  end if
3419  wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
3420  wrc(ixo^s,psi_) = wlp(ixo^s,psi_)
3421  wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
3422  wlc(ixo^s,psi_) = wlp(ixo^s,psi_)
3423  ! modify total energy according to the change of magnetic field
3424  if(total_energy) then
3425  wrc(ixo^s,e_)=wrc(ixo^s,e_)+half*wrc(ixo^s,mag(idir))**2
3426  wlc(ixo^s,e_)=wlc(ixo^s,e_)+half*wlc(ixo^s,mag(idir))**2
3427  end if
3428  end if
3429 
3430  if(associated(usr_set_wlr)) call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
3431 
3432  end subroutine mhd_modify_wlr
3433 
3434  subroutine mhd_boundary_adjust(igrid,psb)
3436  integer, intent(in) :: igrid
3437  type(state), target :: psb(max_blocks)
3438 
3439  integer :: iB, idims, iside, ixO^L, i^D
3440 
3441  block=>ps(igrid)
3442  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
3443  do idims=1,ndim
3444  ! to avoid using as yet unknown corner info in more than 1D, we
3445  ! fill only interior mesh ranges of the ghost cell ranges at first,
3446  ! and progressively enlarge the ranges to include corners later
3447  do iside=1,2
3448  i^d=kr(^d,idims)*(2*iside-3);
3449  if (neighbor_type(i^d,igrid)/=1) cycle
3450  ib=(idims-1)*2+iside
3451  if(.not.boundary_divbfix(ib)) cycle
3452  if(any(typeboundary(:,ib)==bc_special)) then
3453  ! MF nonlinear force-free B field extrapolation and data driven
3454  ! require normal B of the first ghost cell layer to be untouched by
3455  ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
3456  select case (idims)
3457  {case (^d)
3458  if (iside==2) then
3459  ! maximal boundary
3460  ixomin^dd=ixghi^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixglo^dd;
3461  ixomax^dd=ixghi^dd;
3462  else
3463  ! minimal boundary
3464  ixomin^dd=ixglo^dd;
3465  ixomax^dd=ixglo^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixghi^dd;
3466  end if \}
3467  end select
3468  call fixdivb_boundary(ixg^ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
3469  end if
3470  end do
3471  end do
3472 
3473  end subroutine mhd_boundary_adjust
3474 
3475  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
3477 
3478  integer, intent(in) :: ixG^L,ixO^L,iB
3479  double precision, intent(inout) :: w(ixG^S,1:nw)
3480  double precision, intent(in) :: x(ixG^S,1:ndim)
3481 
3482  double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
3483  integer :: ix^D,ixF^L
3484 
3485  select case(ib)
3486  case(1)
3487  ! 2nd order CD for divB=0 to set normal B component better
3488  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
3489  {^iftwod
3490  ixfmin1=ixomin1+1
3491  ixfmax1=ixomax1+1
3492  ixfmin2=ixomin2+1
3493  ixfmax2=ixomax2-1
3494  if(slab_uniform) then
3495  dx1x2=dxlevel(1)/dxlevel(2)
3496  do ix1=ixfmax1,ixfmin1,-1
3497  w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
3498  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
3499  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
3500  enddo
3501  else
3502  do ix1=ixfmax1,ixfmin1,-1
3503  w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
3504  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
3505  +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
3506  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
3507  -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
3508  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
3509  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
3510  end do
3511  end if
3512  }
3513  {^ifthreed
3514  ixfmin1=ixomin1+1
3515  ixfmax1=ixomax1+1
3516  ixfmin2=ixomin2+1
3517  ixfmax2=ixomax2-1
3518  ixfmin3=ixomin3+1
3519  ixfmax3=ixomax3-1
3520  if(slab_uniform) then
3521  dx1x2=dxlevel(1)/dxlevel(2)
3522  dx1x3=dxlevel(1)/dxlevel(3)
3523  do ix1=ixfmax1,ixfmin1,-1
3524  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
3525  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
3526  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
3527  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
3528  +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
3529  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
3530  end do
3531  else
3532  do ix1=ixfmax1,ixfmin1,-1
3533  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
3534  ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
3535  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
3536  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
3537  +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
3538  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
3539  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
3540  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
3541  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
3542  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
3543  +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
3544  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
3545  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
3546  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
3547  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
3548  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
3549  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
3550  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
3551  end do
3552  end if
3553  }
3554  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
3555  case(2)
3556  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
3557  {^iftwod
3558  ixfmin1=ixomin1-1
3559  ixfmax1=ixomax1-1
3560  ixfmin2=ixomin2+1
3561  ixfmax2=ixomax2-1
3562  if(slab_uniform) then
3563  dx1x2=dxlevel(1)/dxlevel(2)
3564  do ix1=ixfmin1,ixfmax1
3565  w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
3566  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
3567  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
3568  enddo
3569  else
3570  do ix1=ixfmin1,ixfmax1
3571  w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
3572  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
3573  -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
3574  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
3575  +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
3576  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
3577  /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
3578  end do
3579  end if
3580  }
3581  {^ifthreed
3582  ixfmin1=ixomin1-1
3583  ixfmax1=ixomax1-1
3584  ixfmin2=ixomin2+1
3585  ixfmax2=ixomax2-1
3586  ixfmin3=ixomin3+1
3587  ixfmax3=ixomax3-1
3588  if(slab_uniform) then
3589  dx1x2=dxlevel(1)/dxlevel(2)
3590  dx1x3=dxlevel(1)/dxlevel(3)
3591  do ix1=ixfmin1,ixfmax1
3592  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
3593  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
3594  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
3595  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
3596  -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
3597  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
3598  end do
3599  else
3600  do ix1=ixfmin1,ixfmax1
3601  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
3602  ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
3603  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
3604  block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
3605  -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
3606  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
3607  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
3608  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
3609  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
3610  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
3611  -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
3612  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
3613  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
3614  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
3615  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
3616  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
3617  /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
3618  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
3619  end do
3620  end if
3621  }
3622  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
3623  case(3)
3624  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
3625  {^iftwod
3626  ixfmin1=ixomin1+1
3627  ixfmax1=ixomax1-1
3628  ixfmin2=ixomin2+1
3629  ixfmax2=ixomax2+1
3630  if(slab_uniform) then
3631  dx2x1=dxlevel(2)/dxlevel(1)
3632  do ix2=ixfmax2,ixfmin2,-1
3633  w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
3634  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
3635  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
3636  enddo
3637  else
3638  do ix2=ixfmax2,ixfmin2,-1
3639  w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
3640  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
3641  +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
3642  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
3643  -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
3644  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
3645  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
3646  end do
3647  end if
3648  }
3649  {^ifthreed
3650  ixfmin1=ixomin1+1
3651  ixfmax1=ixomax1-1
3652  ixfmin3=ixomin3+1
3653  ixfmax3=ixomax3-1
3654  ixfmin2=ixomin2+1
3655  ixfmax2=ixomax2+1
3656  if(slab_uniform) then
3657  dx2x1=dxlevel(2)/dxlevel(1)
3658  dx2x3=dxlevel(2)/dxlevel(3)
3659  do ix2=ixfmax2,ixfmin2,-1
3660  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
3661  ix2+1,ixfmin3:ixfmax3,mag(2)) &
3662  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
3663  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
3664  +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
3665  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
3666  end do
3667  else
3668  do ix2=ixfmax2,ixfmin2,-1
3669  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
3670  ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
3671  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
3672  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
3673  +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
3674  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
3675  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
3676  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
3677  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
3678  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
3679  +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
3680  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
3681  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
3682  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
3683  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
3684  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
3685  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
3686  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
3687  end do
3688  end if
3689  }
3690  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
3691  case(4)
3692  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
3693  {^iftwod
3694  ixfmin1=ixomin1+1
3695  ixfmax1=ixomax1-1
3696  ixfmin2=ixomin2-1
3697  ixfmax2=ixomax2-1
3698  if(slab_uniform) then
3699  dx2x1=dxlevel(2)/dxlevel(1)
3700  do ix2=ixfmin2,ixfmax2
3701  w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
3702  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
3703  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
3704  end do
3705  else
3706  do ix2=ixfmin2,ixfmax2
3707  w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
3708  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
3709  -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
3710  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
3711  +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
3712  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
3713  /block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
3714  end do
3715  end if
3716  }
3717  {^ifthreed
3718  ixfmin1=ixomin1+1
3719  ixfmax1=ixomax1-1
3720  ixfmin3=ixomin3+1
3721  ixfmax3=ixomax3-1
3722  ixfmin2=ixomin2-1
3723  ixfmax2=ixomax2-1
3724  if(slab_uniform) then
3725  dx2x1=dxlevel(2)/dxlevel(1)
3726  dx2x3=dxlevel(2)/dxlevel(3)
3727  do ix2=ixfmin2,ixfmax2
3728  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
3729  ix2-1,ixfmin3:ixfmax3,mag(2)) &
3730  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
3731  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
3732  -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
3733  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
3734  end do
3735  else
3736  do ix2=ixfmin2,ixfmax2
3737  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
3738  ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
3739  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
3740  block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
3741  -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
3742  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
3743  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
3744  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
3745  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
3746  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
3747  -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
3748  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
3749  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
3750  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
3751  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
3752  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
3753  /block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
3754  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
3755  end do
3756  end if
3757  }
3758  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
3759  {^ifthreed
3760  case(5)
3761  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
3762  ixfmin1=ixomin1+1
3763  ixfmax1=ixomax1-1
3764  ixfmin2=ixomin2+1
3765  ixfmax2=ixomax2-1
3766  ixfmin3=ixomin3+1
3767  ixfmax3=ixomax3+1
3768  if(slab_uniform) then
3769  dx3x1=dxlevel(3)/dxlevel(1)
3770  dx3x2=dxlevel(3)/dxlevel(2)
3771  do ix3=ixfmax3,ixfmin3,-1
3772  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
3773  ixfmin2:ixfmax2,ix3+1,mag(3)) &
3774  +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
3775  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
3776  +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
3777  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
3778  end do
3779  else
3780  do ix3=ixfmax3,ixfmin3,-1
3781  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
3782  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
3783  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
3784  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
3785  +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
3786  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
3787  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
3788  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
3789  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
3790  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
3791  +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
3792  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
3793  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
3794  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
3795  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
3796  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
3797  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
3798  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
3799  end do
3800  end if
3801  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
3802  case(6)
3803  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
3804  ixfmin1=ixomin1+1
3805  ixfmax1=ixomax1-1
3806  ixfmin2=ixomin2+1
3807  ixfmax2=ixomax2-1
3808  ixfmin3=ixomin3-1
3809  ixfmax3=ixomax3-1
3810  if(slab_uniform) then
3811  dx3x1=dxlevel(3)/dxlevel(1)
3812  dx3x2=dxlevel(3)/dxlevel(2)
3813  do ix3=ixfmin3,ixfmax3
3814  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
3815  ixfmin2:ixfmax2,ix3-1,mag(3)) &
3816  -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
3817  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
3818  -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
3819  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
3820  end do
3821  else
3822  do ix3=ixfmin3,ixfmax3
3823  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
3824  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
3825  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
3826  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
3827  -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
3828  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
3829  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
3830  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
3831  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
3832  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
3833  -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
3834  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
3835  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
3836  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
3837  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
3838  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
3839  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
3840  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
3841  end do
3842  end if
3843  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
3844  }
3845  case default
3846  call mpistop("Special boundary is not defined for this region")
3847  end select
3848 
3849  end subroutine fixdivb_boundary
3850 
3851  {^nooned
3852  subroutine mhd_clean_divb_multigrid(qdt, qt, active)
3853  use mod_forest
3856  use mod_geometry
3857 
3858  double precision, intent(in) :: qdt !< Current time step
3859  double precision, intent(in) :: qt !< Current time
3860  logical, intent(inout) :: active !< Output if the source is active
3861  integer :: iigrid, igrid, id
3862  integer :: n, nc, lvl, ix^l, ixc^l, idim
3863  type(tree_node), pointer :: pnode
3864  double precision :: tmp(ixg^t), grad(ixg^t, ndim)
3865  double precision :: res
3866  double precision, parameter :: max_residual = 1d-3
3867  double precision, parameter :: residual_reduction = 1d-10
3868  integer, parameter :: max_its = 50
3869  double precision :: residual_it(max_its), max_divb
3870 
3871  mg%operator_type = mg_laplacian
3872 
3873  ! Set boundary conditions
3874  do n = 1, 2*ndim
3875  idim = (n+1)/2
3876  select case (typeboundary(mag(idim), n))
3877  case (bc_symm)
3878  ! d/dx B = 0, take phi = 0
3879  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
3880  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3881  case (bc_asymm)
3882  ! B = 0, so grad(phi) = 0
3883  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
3884  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3885  case (bc_cont)
3886  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
3887  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3888  case (bc_special)
3889  ! Assume Dirichlet boundary conditions, derivative zero
3890  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
3891  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3892  case (bc_periodic)
3893  ! Nothing to do here
3894  case default
3895  write(*,*) "mhd_clean_divb_multigrid warning: unknown boundary type"
3896  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
3897  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3898  end select
3899  end do
3900 
3901  ix^l=ixm^ll^ladd1;
3902  max_divb = 0.0d0
3903 
3904  ! Store divergence of B as right-hand side
3905  do iigrid = 1, igridstail
3906  igrid = igrids(iigrid);
3907  pnode => igrid_to_node(igrid, mype)%node
3908  id = pnode%id
3909  lvl = mg%boxes(id)%lvl
3910  nc = mg%box_size_lvl(lvl)
3911 
3912  ! Geometry subroutines expect this to be set
3913  block => ps(igrid)
3914  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
3915 
3916  call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^ll, ixm^ll, tmp, &
3918  mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(ixm^t)
3919  max_divb = max(max_divb, maxval(abs(tmp(ixm^t))))
3920  end do
3921 
3922  ! Solve laplacian(phi) = divB
3923  if(stagger_grid) then
3924  call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
3925  mpi_max, icomm, ierrmpi)
3926 
3927  if (mype == 0) print *, "Performing multigrid divB cleaning"
3928  if (mype == 0) print *, "iteration vs residual"
3929  ! Solve laplacian(phi) = divB
3930  do n = 1, max_its
3931  call mg_fas_fmg(mg, n>1, max_res=residual_it(n))
3932  if (mype == 0) write(*, "(I4,E11.3)") n, residual_it(n)
3933  if (residual_it(n) < residual_reduction * max_divb) exit
3934  end do
3935  if (mype == 0 .and. n > max_its) then
3936  print *, "divb_multigrid warning: not fully converged"
3937  print *, "current amplitude of divb: ", residual_it(max_its)
3938  print *, "multigrid smallest grid: ", &
3939  mg%domain_size_lvl(:, mg%lowest_lvl)
3940  print *, "note: smallest grid ideally has <= 8 cells"
3941  print *, "multigrid dx/dy/dz ratio: ", mg%dr(:, 1)/mg%dr(1, 1)
3942  print *, "note: dx/dy/dz should be similar"
3943  end if
3944  else
3945  do n = 1, max_its
3946  call mg_fas_vcycle(mg, max_res=res)
3947  if (res < max_residual) exit
3948  end do
3949  if (res > max_residual) call mpistop("divb_multigrid: no convergence")
3950  end if
3951 
3952 
3953  ! Correct the magnetic field
3954  do iigrid = 1, igridstail
3955  igrid = igrids(iigrid);
3956  pnode => igrid_to_node(igrid, mype)%node
3957  id = pnode%id
3958 
3959  ! Geometry subroutines expect this to be set
3960  block => ps(igrid)
3961  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
3962 
3963  ! Compute the gradient of phi
3964  tmp(ix^s) = mg%boxes(id)%cc({:,}, mg_iphi)
3965 
3966  if(stagger_grid) then
3967  do idim =1, ndim
3968  ixcmin^d=ixmlo^d-kr(idim,^d);
3969  ixcmax^d=ixmhi^d;
3970  call gradientx(tmp,ps(igrid)%x,ixg^ll,ixc^l,idim,grad(ixg^t,idim),.false.)
3971  ! Apply the correction B* = B - gradient(phi)
3972  ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
3973  end do
3974  ! store cell-center magnetic energy
3975  tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
3976  ! change cell-center magnetic field
3977  call mhd_face_to_center(ixm^ll,ps(igrid))
3978  else
3979  do idim = 1, ndim
3980  call gradient(tmp,ixg^ll,ixm^ll,idim,grad(ixg^t, idim))
3981  end do
3982  ! store cell-center magnetic energy
3983  tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
3984  ! Apply the correction B* = B - gradient(phi)
3985  ps(igrid)%w(ixm^t, mag(1:ndim)) = &
3986  ps(igrid)%w(ixm^t, mag(1:ndim)) - grad(ixm^t, :)
3987  end if
3988 
3989  if(total_energy) then
3990  ! Determine magnetic energy difference
3991  tmp(ixm^t) = 0.5_dp * (sum(ps(igrid)%w(ixm^t, &
3992  mag(1:ndim))**2, dim=ndim+1) - tmp(ixm^t))
3993  ! Keep thermal pressure the same
3994  ps(igrid)%w(ixm^t, e_) = ps(igrid)%w(ixm^t, e_) + tmp(ixm^t)
3995  end if
3996  end do
3997 
3998  active = .true.
3999 
4000  end subroutine mhd_clean_divb_multigrid
4001  }
4002 
4003  subroutine mhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
4005 
4006  integer, intent(in) :: ixI^L, ixO^L
4007  double precision, intent(in) :: qt,qdt
4008  ! cell-center primitive variables
4009  double precision, intent(in) :: wprim(ixI^S,1:nw)
4010  type(state) :: sCT, s
4011  type(ct_velocity) :: vcts
4012  double precision, intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
4013  double precision, intent(inout) :: fE(ixI^S,7-2*ndim:3)
4014 
4015  select case(type_ct)
4016  case('average')
4017  call update_faces_average(ixi^l,ixo^l,qt,qdt,fc,fe,sct,s)
4018  case('uct_contact')
4019  call update_faces_contact(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,s,vcts)
4020  case('uct_hll')
4021  call update_faces_hll(ixi^l,ixo^l,qt,qdt,fe,sct,s,vcts)
4022  case default
4023  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
4024  end select
4025 
4026  end subroutine mhd_update_faces
4027 
4028  !> get electric field though averaging neighors to update faces in CT
4029  subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
4031  use mod_usr_methods
4032 
4033  integer, intent(in) :: ixI^L, ixO^L
4034  double precision, intent(in) :: qt, qdt
4035  type(state) :: sCT, s
4036  double precision, intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
4037  double precision, intent(inout) :: fE(ixI^S,7-2*ndim:3)
4038 
4039  integer :: hxC^L,ixC^L,jxC^L,ixCm^L
4040  integer :: idim1,idim2,idir,iwdim1,iwdim2
4041  double precision :: circ(ixI^S,1:ndim)
4042  ! non-ideal electric field on cell edges
4043  double precision, dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
4044 
4045  associate(bfaces=>s%ws,x=>s%x)
4046 
4047  ! Calculate contribution to FEM of each edge,
4048  ! that is, estimate value of line integral of
4049  ! electric field in the positive idir direction.
4050  ixcmax^d=ixomax^d;
4051  ixcmin^d=ixomin^d-1;
4052 
4053  ! if there is resistivity, get eta J
4054  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
4055 
4056  ! if there is ambipolar diffusion, get E_ambi
4057  if(mhd_ambipolar_exp) call get_ambipolar_electric_field(ixi^l,ixo^l,sct%w,x,e_ambi)
4058 
4059  fe=zero
4060 
4061  do idim1=1,ndim
4062  iwdim1 = mag(idim1)
4063  do idim2=1,ndim
4064  iwdim2 = mag(idim2)
4065  do idir=7-2*ndim,3! Direction of line integral
4066  ! Allow only even permutations
4067  if (lvc(idim1,idim2,idir)==1) then
4068  ! Assemble indices
4069  jxc^l=ixc^l+kr(idim1,^d);
4070  hxc^l=ixc^l+kr(idim2,^d);
4071  ! Interpolate to edges
4072  fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
4073  -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
4074 
4075  ! add resistive electric field at cell edges E=-vxB+eta J
4076  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
4077  ! add ambipolar electric field
4078  if(mhd_ambipolar_exp) fe(ixc^s,idir)=fe(ixc^s,idir)+e_ambi(ixc^s,idir)
4079 
4080  fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
4081 
4082  if (.not.slab) then
4083  where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
4084  fe(ixc^s,idir)=zero
4085  end where
4086  end if
4087  end if
4088  end do
4089  end do
4090  end do
4091 
4092  ! allow user to change inductive electric field, especially for boundary driven applications
4093  if(associated(usr_set_electric_field)) &
4094  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
4095 
4096  circ(ixi^s,1:ndim)=zero
4097 
4098  ! Calculate circulation on each face
4099 
4100  do idim1=1,ndim ! Coordinate perpendicular to face
4101  do idim2=1,ndim
4102  do idir=7-2*ndim,3 ! Direction of line integral
4103  ! Assemble indices
4104  hxc^l=ixc^l-kr(idim2,^d);
4105  ! Add line integrals in direction idir
4106  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4107  +lvc(idim1,idim2,idir)&
4108  *(fe(ixc^s,idir)&
4109  -fe(hxc^s,idir))
4110  end do
4111  end do
4112  end do
4113 
4114  ! Divide by the area of the face to get dB/dt
4115  do idim1=1,ndim
4116  ixcmax^d=ixomax^d;
4117  ixcmin^d=ixomin^d-kr(idim1,^d);
4118  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
4119  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
4120  elsewhere
4121  circ(ixc^s,idim1)=zero
4122  end where
4123  ! Time update
4124  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
4125  end do
4126 
4127  end associate
4128 
4129  end subroutine update_faces_average
4130 
4131  !> update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
4132  subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
4134  use mod_usr_methods
4135 
4136  integer, intent(in) :: ixI^L, ixO^L
4137  double precision, intent(in) :: qt, qdt
4138  ! cell-center primitive variables
4139  double precision, intent(in) :: wp(ixI^S,1:nw)
4140  type(state) :: sCT, s
4141  type(ct_velocity) :: vcts
4142  double precision, intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
4143  double precision, intent(inout) :: fE(ixI^S,7-2*ndim:3)
4144 
4145  double precision :: circ(ixI^S,1:ndim)
4146  ! electric field at cell centers
4147  double precision :: ECC(ixI^S,7-2*ndim:3)
4148  ! gradient of E at left and right side of a cell face
4149  double precision :: EL(ixI^S),ER(ixI^S)
4150  ! gradient of E at left and right side of a cell corner
4151  double precision :: ELC(ixI^S),ERC(ixI^S)
4152  ! non-ideal electric field on cell edges
4153  double precision, dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
4154  ! total magnetic field at cell centers
4155  double precision :: Btot(ixI^S,1:ndim)
4156  integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
4157  integer :: idim1,idim2,idir,iwdim1,iwdim2
4158 
4159  associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
4160 
4161  if(b0field) then
4162  btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))+block%B0(ixi^s,1:ndim,0)
4163  else
4164  btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))
4165  end if
4166  ecc=0.d0
4167  ! Calculate electric field at cell centers
4168  do idim1=1,ndim; do idim2=1,ndim; do idir=7-2*ndim,3
4169  if(lvc(idim1,idim2,idir)==1)then
4170  ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,mom(idim2))
4171  else if(lvc(idim1,idim2,idir)==-1) then
4172  ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,mom(idim2))
4173  endif
4174  enddo; enddo; enddo
4175 
4176  ! if there is resistivity, get eta J
4177  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
4178 
4179  ! if there is ambipolar diffusion, get E_ambi
4180  if(mhd_ambipolar_exp) call get_ambipolar_electric_field(ixi^l,ixo^l,sct%w,x,e_ambi)
4181 
4182  ! Calculate contribution to FEM of each edge,
4183  ! that is, estimate value of line integral of
4184  ! electric field in the positive idir direction.
4185  fe=zero
4186  ! evaluate electric field along cell edges according to equation (41)
4187  do idim1=1,ndim
4188  iwdim1 = mag(idim1)
4189  do idim2=1,ndim
4190  iwdim2 = mag(idim2)
4191  do idir=7-2*ndim,3 ! Direction of line integral
4192  ! Allow only even permutations
4193  if (lvc(idim1,idim2,idir)==1) then
4194  ixcmax^d=ixomax^d;
4195  ixcmin^d=ixomin^d+kr(idir,^d)-1;
4196  ! Assemble indices
4197  jxc^l=ixc^l+kr(idim1,^d);
4198  hxc^l=ixc^l+kr(idim2,^d);
4199  ! average cell-face electric field to cell edges
4200  fe(ixc^s,idir)=quarter*&
4201  (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
4202  -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
4203 
4204  ! add slope in idim2 direction from equation (50)
4205  ixamin^d=ixcmin^d;
4206  ixamax^d=ixcmax^d+kr(idim1,^d);
4207  el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
4208  hxc^l=ixa^l+kr(idim2,^d);
4209  er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
4210  where(vnorm(ixc^s,idim1)>0.d0)
4211  elc(ixc^s)=el(ixc^s)
4212  else where(vnorm(ixc^s,idim1)<0.d0)
4213  elc(ixc^s)=el(jxc^s)
4214  else where
4215  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
4216  end where
4217  hxc^l=ixc^l+kr(idim2,^d);
4218  where(vnorm(hxc^s,idim1)>0.d0)
4219  erc(ixc^s)=er(ixc^s)
4220  else where(vnorm(hxc^s,idim1)<0.d0)
4221  erc(ixc^s)=er(jxc^s)
4222  else where
4223  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
4224  end where
4225  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
4226 
4227  ! add slope in idim1 direction from equation (50)
4228  jxc^l=ixc^l+kr(idim2,^d);
4229  ixamin^d=ixcmin^d;
4230  ixamax^d=ixcmax^d+kr(idim2,^d);
4231  el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
4232  hxc^l=ixa^l+kr(idim1,^d);
4233  er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
4234  where(vnorm(ixc^s,idim2)>0.d0)
4235  elc(ixc^s)=el(ixc^s)
4236  else where(vnorm(ixc^s,idim2)<0.d0)
4237  elc(ixc^s)=el(jxc^s)
4238  else where
4239  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
4240  end where
4241  hxc^l=ixc^l+kr(idim1,^d);
4242  where(vnorm(hxc^s,idim2)>0.d0)
4243  erc(ixc^s)=er(ixc^s)
4244  else where(vnorm(hxc^s,idim2)<0.d0)
4245  erc(ixc^s)=er(jxc^s)
4246  else where
4247  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
4248  end where
4249  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
4250 
4251  ! add resistive electric field at cell edges E=-vxB+eta J
4252  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
4253  ! add ambipolar electric field
4254  if(mhd_ambipolar_exp) fe(ixc^s,idir)=fe(ixc^s,idir)+e_ambi(ixc^s,idir)
4255 
4256  ! times time step and edge length
4257  fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
4258  if (.not.slab) then
4259  where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
4260  fe(ixc^s,idir)=zero
4261  end where
4262  end if
4263  end if
4264  end do
4265  end do
4266  end do
4267 
4268  ! allow user to change inductive electric field, especially for boundary driven applications
4269  if(associated(usr_set_electric_field)) &
4270  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
4271 
4272  circ(ixi^s,1:ndim)=zero
4273 
4274  ! Calculate circulation on each face
4275  do idim1=1,ndim ! Coordinate perpendicular to face
4276  ixcmax^d=ixomax^d;
4277  ixcmin^d=ixomin^d-kr(idim1,^d);
4278  do idim2=1,ndim
4279  do idir=7-2*ndim,3 ! Direction of line integral
4280  ! Assemble indices
4281  hxc^l=ixc^l-kr(idim2,^d);
4282  ! Add line integrals in direction idir
4283  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4284  +lvc(idim1,idim2,idir)&
4285  *(fe(ixc^s,idir)&
4286  -fe(hxc^s,idir))
4287  end do
4288  end do
4289  ! Divide by the area of the face to get dB/dt
4290  ixcmax^d=ixomax^d;
4291  ixcmin^d=ixomin^d-kr(idim1,^d);
4292  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
4293  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
4294  elsewhere
4295  circ(ixc^s,idim1)=zero
4296  end where
4297  ! Time update cell-face magnetic field component
4298  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
4299  end do
4300 
4301  end associate
4302 
4303  end subroutine update_faces_contact
4304 
4305  !> update faces
4306  subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
4309  use mod_usr_methods
4310 
4311  integer, intent(in) :: ixI^L, ixO^L
4312  double precision, intent(in) :: qt, qdt
4313  double precision, intent(inout) :: fE(ixI^S,7-2*ndim:3)
4314  type(state) :: sCT, s
4315  type(ct_velocity) :: vcts
4316 
4317  double precision :: vtilL(ixI^S,2)
4318  double precision :: vtilR(ixI^S,2)
4319  double precision :: bfacetot(ixI^S,ndim)
4320  double precision :: btilL(ixI^S,ndim)
4321  double precision :: btilR(ixI^S,ndim)
4322  double precision :: cp(ixI^S,2)
4323  double precision :: cm(ixI^S,2)
4324  double precision :: circ(ixI^S,1:ndim)
4325  ! non-ideal electric field on cell edges
4326  double precision, dimension(ixI^S,7-2*ndim:3) :: E_resi, E_ambi
4327  integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
4328  integer :: idim1,idim2,idir
4329 
4330  associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
4331  cbarmax=>vcts%cbarmax)
4332 
4333  ! Calculate contribution to FEM of each edge,
4334  ! that is, estimate value of line integral of
4335  ! electric field in the positive idir direction.
4336 
4337  ! Loop over components of electric field
4338 
4339  ! idir: electric field component we need to calculate
4340  ! idim1: directions in which we already performed the reconstruction
4341  ! idim2: directions in which we perform the reconstruction
4342 
4343  ! if there is resistivity, get eta J
4344  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
4345 
4346  ! if there is ambipolar diffusion, get E_ambi
4347  if(mhd_ambipolar_exp) call get_ambipolar_electric_field(ixi^l,ixo^l,sct%w,x,e_ambi)
4348 
4349  fe=zero
4350 
4351  do idir=7-2*ndim,3
4352  ! Indices
4353  ! idir: electric field component
4354  ! idim1: one surface
4355  ! idim2: the other surface
4356  ! cyclic permutation: idim1,idim2,idir=1,2,3
4357  ! Velocity components on the surface
4358  ! follow cyclic premutations:
4359  ! Sx(1),Sx(2)=y,z ; Sy(1),Sy(2)=z,x ; Sz(1),Sz(2)=x,y
4360 
4361  ixcmax^d=ixomax^d;
4362  ixcmin^d=ixomin^d-1+kr(idir,^d);
4363 
4364  ! Set indices and directions
4365  idim1=mod(idir,3)+1
4366  idim2=mod(idir+1,3)+1
4367 
4368  jxc^l=ixc^l+kr(idim1,^d);
4369  ixcp^l=ixc^l+kr(idim2,^d);
4370 
4371  ! Reconstruct transverse transport velocities
4372  call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
4373  vtill(ixi^s,2),vtilr(ixi^s,2))
4374 
4375  call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
4376  vtill(ixi^s,1),vtilr(ixi^s,1))
4377 
4378  ! Reconstruct magnetic fields
4379  ! Eventhough the arrays are larger, reconstruct works with
4380  ! the limits ixG.
4381  if(b0field) then
4382  bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+block%B0(ixi^s,idim1,idim1)
4383  bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+block%B0(ixi^s,idim2,idim2)
4384  else
4385  bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
4386  bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
4387  end if
4388  call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
4389  btill(ixi^s,idim1),btilr(ixi^s,idim1))
4390 
4391  call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
4392  btill(ixi^s,idim2),btilr(ixi^s,idim2))
4393 
4394  ! Take the maximum characteristic
4395 
4396  cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
4397  cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
4398 
4399  cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
4400  cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
4401 
4402 
4403  ! Calculate eletric field
4404  fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
4405  + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
4406  - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
4407  /(cp(ixc^s,1)+cm(ixc^s,1)) &
4408  +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
4409  + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
4410  - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
4411  /(cp(ixc^s,2)+cm(ixc^s,2))
4412 
4413  ! add resistive electric field at cell edges E=-vxB+eta J
4414  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
4415  ! add ambipolar electric field
4416  if(mhd_ambipolar_exp) fe(ixc^s,idir)=fe(ixc^s,idir)+e_ambi(ixc^s,idir)
4417 
4418  fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
4419 
4420  if (.not.slab) then
4421  where(abs(x(ixc^s,r_)+half*dxlevel(r_)).lt.1.0d-9)
4422  fe(ixc^s,idir)=zero
4423  end where
4424  end if
4425 
4426  end do
4427 
4428  ! allow user to change inductive electric field, especially for boundary driven applications
4429  if(associated(usr_set_electric_field)) &
4430  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
4431 
4432  circ(ixi^s,1:ndim)=zero
4433 
4434  ! Calculate circulation on each face: interal(fE dot dl)
4435 
4436  do idim1=1,ndim ! Coordinate perpendicular to face
4437  ixcmax^d=ixomax^d;
4438  ixcmin^d=ixomin^d-kr(idim1,^d);
4439  do idim2=1,ndim
4440  do idir=7-2*ndim,3 ! Direction of line integral
4441  ! Assemble indices
4442  hxc^l=ixc^l-kr(idim2,^d);
4443  ! Add line integrals in direction idir
4444  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4445  +lvc(idim1,idim2,idir)&
4446  *(fe(ixc^s,idir)&
4447  -fe(hxc^s,idir))
4448  end do
4449  end do
4450  end do
4451 
4452  ! Divide by the area of the face to get dB/dt
4453  do idim1=1,ndim
4454  ixcmax^d=ixomax^d;
4455  ixcmin^d=ixomin^d-kr(idim1,^d);
4456  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
4457  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
4458  elsewhere
4459  circ(ixc^s,idim1)=zero
4460  end where
4461  ! Time update
4462  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
4463  end do
4464 
4465  end associate
4466  end subroutine update_faces_hll
4467 
4468  !> calculate eta J at cell edges
4469  subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
4471  use mod_usr_methods
4472  use mod_geometry
4473 
4474  integer, intent(in) :: ixI^L, ixO^L
4475  type(state), intent(in) :: sCT, s
4476  ! current on cell edges
4477  double precision :: jce(ixI^S,7-2*ndim:3)
4478 
4479  ! current on cell centers
4480  double precision :: jcc(ixI^S,7-2*ndir:3)
4481  ! location at cell faces
4482  double precision :: xs(ixGs^T,1:ndim)
4483  ! resistivity
4484  double precision :: eta(ixI^S)
4485  double precision :: gradi(ixGs^T)
4486  integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
4487 
4488  associate(x=>s%x,dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
4489  ! calculate current density at cell edges
4490  jce=0.d0
4491  do idim1=1,ndim
4492  do idim2=1,ndim
4493  do idir=7-2*ndim,3
4494  if (lvc(idim1,idim2,idir)==0) cycle
4495  ixcmax^d=ixomax^d;
4496  ixcmin^d=ixomin^d+kr(idir,^d)-1;
4497  ixbmax^d=ixcmax^d-kr(idir,^d)+1;
4498  ixbmin^d=ixcmin^d;
4499  ! current at transverse faces
4500  xs(ixb^s,:)=x(ixb^s,:)
4501  xs(ixb^s,idim2)=x(ixb^s,idim2)+half*dx(ixb^s,idim2)
4502  call gradientx(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi,.true.)
4503  if (lvc(idim1,idim2,idir)==1) then
4504  jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
4505  else
4506  jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
4507  end if
4508  end do
4509  end do
4510  end do
4511  ! get resistivity
4512  if(mhd_eta>zero)then
4513  jce(ixi^s,:)=jce(ixi^s,:)*mhd_eta
4514  else
4515  ixa^l=ixo^l^ladd1;
4516  call get_current(wct,ixi^l,ixa^l,idirmin,jcc)
4517  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,jcc,eta)
4518  ! calcuate eta on cell edges
4519  do idir=7-2*ndim,3
4520  ixcmax^d=ixomax^d;
4521  ixcmin^d=ixomin^d+kr(idir,^d)-1;
4522  jcc(ixc^s,idir)=0.d0
4523  {do ix^db=0,1\}
4524  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
4525  ixamin^d=ixcmin^d+ix^d;
4526  ixamax^d=ixcmax^d+ix^d;
4527  jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
4528  {end do\}
4529  jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
4530  jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
4531  enddo
4532  end if
4533 
4534  end associate
4535  end subroutine get_resistive_electric_field
4536 
4537  !> get ambipolar electric field on cell edges
4538  subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
4540 
4541  integer, intent(in) :: ixI^L, ixO^L
4542  double precision, intent(in) :: w(ixI^S,1:nw)
4543  double precision, intent(in) :: x(ixI^S,1:ndim)
4544  double precision, intent(out) :: fE(ixI^S,7-2*ndim:3)
4545 
4546  double precision :: jxbxb(ixI^S,1:3)
4547  integer :: idir,ixA^L,ixC^L,ix^D
4548 
4549  ixa^l=ixo^l^ladd1;
4550  call mhd_get_jxbxb(w,x,ixi^l,ixa^l,jxbxb)
4551  ! calcuate electric field on cell edges from cell centers
4552  do idir=7-2*ndim,3
4553  !set electric field in jxbxb: E=nuA * jxbxb, where nuA=-etaA/rho^2
4554  !jxbxb(ixA^S,i) = -(mhd_eta_ambi/w(ixA^S, rho_)**2) * jxbxb(ixA^S,i)
4555  call multiplyambicoef(ixi^l,ixa^l,jxbxb(ixi^s,idir),w,x)
4556  ixcmax^d=ixomax^d;
4557  ixcmin^d=ixomin^d+kr(idir,^d)-1;
4558  fe(ixc^s,idir)=0.d0
4559  {do ix^db=0,1\}
4560  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
4561  ixamin^d=ixcmin^d+ix^d;
4562  ixamax^d=ixcmax^d+ix^d;
4563  fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
4564  {end do\}
4565  fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
4566  end do
4567 
4568  end subroutine get_ambipolar_electric_field
4569 
4570  !> calculate cell-center values from face-center values
4571  subroutine mhd_face_to_center(ixO^L,s)
4573  ! Non-staggered interpolation range
4574  integer, intent(in) :: ixo^l
4575  type(state) :: s
4576 
4577  integer :: fxo^l, gxo^l, hxo^l, jxo^l, kxo^l, idim
4578 
4579  associate(w=>s%w, ws=>s%ws)
4580 
4581  ! calculate cell-center values from face-center values in 2nd order
4582  do idim=1,ndim
4583  ! Displace index to the left
4584  ! Even if ixI^L is the full size of the w arrays, this is ok
4585  ! because the staggered arrays have an additional place to the left.
4586  hxo^l=ixo^l-kr(idim,^d);
4587  ! Interpolate to cell barycentre using arithmetic average
4588  ! This might be done better later, to make the method less diffusive.
4589  w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
4590  (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
4591  +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
4592  end do
4593 
4594  ! calculate cell-center values from face-center values in 4th order
4595  !do idim=1,ndim
4596  ! gxO^L=ixO^L-2*kr(idim,^D);
4597  ! hxO^L=ixO^L-kr(idim,^D);
4598  ! jxO^L=ixO^L+kr(idim,^D);
4599 
4600  ! ! Interpolate to cell barycentre using fourth order central formula
4601  ! w(ixO^S,mag(idim))=(0.0625d0/s%surface(ixO^S,idim))*&
4602  ! ( -ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
4603  ! +9.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
4604  ! +9.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
4605  ! -ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) )
4606  !end do
4607 
4608  ! calculate cell-center values from face-center values in 6th order
4609  !do idim=1,ndim
4610  ! fxO^L=ixO^L-3*kr(idim,^D);
4611  ! gxO^L=ixO^L-2*kr(idim,^D);
4612  ! hxO^L=ixO^L-kr(idim,^D);
4613  ! jxO^L=ixO^L+kr(idim,^D);
4614  ! kxO^L=ixO^L+2*kr(idim,^D);
4615 
4616  ! ! Interpolate to cell barycentre using sixth order central formula
4617  ! w(ixO^S,mag(idim))=(0.00390625d0/s%surface(ixO^S,idim))* &
4618  ! ( +3.0d0*ws(fxO^S,idim)*s%surfaceC(fxO^S,idim) &
4619  ! -25.0d0*ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
4620  ! +150.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
4621  ! +150.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
4622  ! -25.0d0*ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) &
4623  ! +3.0d0*ws(kxO^S,idim)*s%surfaceC(kxO^S,idim) )
4624  !end do
4625 
4626  end associate
4627 
4628  end subroutine mhd_face_to_center
4629 
4630  !> calculate magnetic field from vector potential
4631  subroutine b_from_vector_potential(ixIs^L, ixI^L, ixO^L, ws, x)
4634 
4635  integer, intent(in) :: ixis^l, ixi^l, ixo^l
4636  double precision, intent(inout) :: ws(ixis^s,1:nws)
4637  double precision, intent(in) :: x(ixi^s,1:ndim)
4638 
4639  double precision :: adummy(ixis^s,1:3)
4640 
4641  call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, adummy)
4642 
4643  end subroutine b_from_vector_potential
4644 
4645 end module mod_mhd_phys
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
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 flux conservation near refinement boundaries.
subroutine, public store_edge(igrid, ixIL, fE, idimLIM)
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
Module with basic grid data structures.
Definition: mod_forest.t:2
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
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges.
Definition: mod_geometry.t:571
integer coordinate
Definition: mod_geometry.t:6
integer, parameter spherical
Definition: mod_geometry.t:10
integer, parameter cylindrical
Definition: mod_geometry.t:9
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:320
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
Definition: mod_geometry.t:626
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:421
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:479
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:364
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer nstep
How many sub-steps the time integrator takes.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
integer istep
Index of the sub-step in a multi-step time integrator.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
character(len=std_len) typediv
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range (within a block with ghost cells)
integer ierrmpi
A global MPI error return code.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical crash
Save a snapshot before crash a run met unphysical values.
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
double precision, dimension(ndim) dxlevel
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for including gravity in (magneto)hydrodynamics simulations.
Definition: mod_gravity.t:2
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_gravity.t:74
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 gravity_init()
Initialize the module.
Definition: mod_gravity.t:27
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
Definition: mod_mhd_phys.t:2
subroutine mhd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities.
logical, public, protected mhd_gravity
Whether gravity is added.
Definition: mod_mhd_phys.t:25
subroutine update_faces_hll(ixIL, ixOL, qt, qdt, fE, sCT, s, vcts)
update faces
subroutine mhd_e_to_ei_aux(ixIL, ixOL, w, x)
Transform total energy to internal energy via eaux as internal energy.
Definition: mod_mhd_phys.t:887
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,...
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
Definition: mod_mhd_phys.t:64
subroutine internal_energy_add_source(qdt, ixIL, ixOL, wCT, w, x, ie)
character(len=std_len), public, protected type_ct
Method type of constrained transport.
Definition: mod_mhd_phys.t:148
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_mhd_phys.t:95
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
subroutine mhd_check_params
Definition: mod_mhd_phys.t:675
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.
subroutine mhd_get_tcutoff(ixIL, ixOL, w, x, Tco_local, Tmax_local)
get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
subroutine update_faces_contact(ixIL, ixOL, qt, qdt, wp, fC, fE, sCT, s, vcts)
update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
subroutine mhd_get_jambi(w, x, ixIL, ixOL, res)
double precision function, dimension(ixo^s) mhd_mag_i_all(w, ixIL, ixOL, idir)
Compute full magnetic field by direction.
subroutine boris_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
Definition: mod_mhd_phys.t:19
double precision, public mhd_adiab
The adiabatic constant.
Definition: mod_mhd_phys.t:124
double precision function, dimension(ixo^s), public mhd_kin_en(w, ixIL, ixOL, inv_rho)
compute kinetic energy
subroutine mhd_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
subroutine, public mhd_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
Definition: mod_mhd_phys.t:776
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
Definition: mod_mhd_phys.t:130
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 add_source_b0split(qdt, ixIL, ixOL, wCT, w, x)
Source terms after split off time-independent magnetic field.
subroutine mhd_get_a2max(w, x, ixIL, ixOL, a2max)
subroutine mhd_get_temperature_from_etot(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the total energy is stored this does not check the values of m...
logical, public, protected mhd_divb_4thorder
Whether divB is computed with a fourth order approximation.
Definition: mod_mhd_phys.t:151
double precision, public mhd_gamma
The adiabatic index.
Definition: mod_mhd_phys.t:121
integer, dimension(:), allocatable, public, protected mag
Indices of the magnetic field.
Definition: mod_mhd_phys.t:104
integer, public, protected mhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
Definition: mod_mhd_phys.t:58
subroutine mhd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L.
subroutine mhd_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
subroutine, public mhd_face_to_center(ixOL, s)
calculate cell-center values from face-center values
subroutine, public mhd_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
Definition: mod_mhd_phys.t:806
subroutine mhd_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
subroutine update_faces_ambipolar(ixIL, ixOL, w, x, ECC, fE, circ)
get ambipolar electric field and the integrals around cell faces
subroutine mhd_ei_to_e_aux(ixIL, ixOL, w, x)
Update eaux and transform internal energy to total energy.
Definition: mod_mhd_phys.t:872
subroutine mhd_get_csound(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
integer, public, protected mhd_n_tracer
Number of tracer species.
Definition: mod_mhd_phys.t:89
subroutine mhd_angmomfix(fC, x, wnew, ixIL, ixOL, idim)
Definition: mod_mhd_phys.t:287
integer, public, protected mhd_trac_type
Which TRAC method is used
Definition: mod_mhd_phys.t:52
subroutine mhd_check_w(primitive, ixIL, ixOL, w, flag)
Definition: mod_mhd_phys.t:747
subroutine add_source_ambipolar_internal_energy(qdt, ixIL, ixOL, wCT, w, x, ie)
Source terms J.E in internal energy. For the ambipolar term E = ambiCoef * JxBxB=ambiCoef * B^2(-J_pe...
subroutine mhd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Definition: mod_mhd_phys.t:67
integer, public, protected paux_
Definition: mod_mhd_phys.t:111
logical, public, protected mhd_hall
Whether Hall-MHD is used.
Definition: mod_mhd_phys.t:28
subroutine rhos_to_e(ixIL, ixOL, w, x)
Convert entropy to energy.
subroutine mhd_energy_synchro(ixIL, ixOL, w, x)
Definition: mod_mhd_phys.t:897
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
subroutine mhd_gamma2_alfven(ixIL, ixOL, w, gamma_A2)
Compute 1/(1+v_A^2/c^2) for Boris' approximation, where v_A is the Alfven velocity.
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
logical, public, protected mhd_ambipolar
Whether Ambipolar term is used.
Definition: mod_mhd_phys.t:31
subroutine mhd_get_ct_velocity(vcts, wLp, wRp, ixIL, ixOL, idim, cmax, cmin)
prepare velocities for ct methods
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.
subroutine mhd_ei_to_e(ixIL, ixOL, w, x)
Transform internal energy to total energy.
Definition: mod_mhd_phys.t:840
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
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:71
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
subroutine mhd_get_jxbxb(w, x, ixIL, ixOL, res)
subroutine mhd_get_h_speed(wprim, x, ixIL, ixOL, idim, Hspeed)
get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
subroutine get_flux_on_cell_face(ixIL, ixOL, ff, src)
use cell-center flux to get cell-face flux and get the source term as the divergence of the flux
logical, public, protected mhd_viscosity
Whether viscosity is added.
Definition: mod_mhd_phys.t:22
subroutine mhd_get_p_total(w, x, ixIL, ixOL, p)
Calculate total pressure within ixO^L including magnetic pressure.
logical, public, protected mhd_energy
Whether an energy equation is used.
Definition: mod_mhd_phys.t:8
logical, public, protected mhd_ambipolar_exp
Whether Ambipolar term is implemented explicitly.
Definition: mod_mhd_phys.t:37
subroutine get_ambipolar_electric_field(ixIL, ixOL, w, x, fE)
get ambipolar electric field on cell edges
logical, public, protected mhd_solve_eaux
Whether auxiliary internal energy is solved.
Definition: mod_mhd_phys.t:61
logical, public, protected mhd_glm
Whether GLM-MHD is used.
Definition: mod_mhd_phys.t:46
logical, public clean_initial_divb
clean initial divB
Definition: mod_mhd_phys.t:169
subroutine update_faces_average(ixIL, ixOL, qt, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
subroutine mhd_update_faces(ixIL, ixOL, qt, qdt, wprim, fC, fE, sCT, s, vcts)
subroutine mhd_get_csound_prim(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
double precision, public mhd_eta
The MHD resistivity.
Definition: mod_mhd_phys.t:127
logical, public divbwave
Add divB wave in Roe solver.
Definition: mod_mhd_phys.t:166
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.
double precision function, dimension(ixo^s) mhd_mag_en(w, ixIL, ixOL)
Compute evolving magnetic energy.
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
Definition: mod_mhd_phys.t:43
double precision, public, protected mhd_trac_mask
Height of the mask used in the TRAC method.
Definition: mod_mhd_phys.t:55
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
Definition: mod_mhd_phys.t:145
subroutine sts_set_source_ambipolar(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
Sets the sources for the ambipolar this is used for the STS method.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
Definition: mod_mhd_phys.t:11
subroutine mhd_getdt_hall(w, x, ixIL, ixOL, dxD, dthall)
subroutine mhd_getv_hall(w, x, ixIL, ixOL, vHall)
double precision function get_ambipolar_dt(w, ixIL, ixOL, dxD, x)
Calculates the explicit dt for the ambipokar term This function is used by both explicit scheme and S...
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_mhd_phys.t:101
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
subroutine mhd_get_temperature_from_eint(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the internal energy is stored.
integer, public, protected eaux_
Indices of auxiliary internal energy.
Definition: mod_mhd_phys.t:110
logical, public, protected mhd_particles
Whether particles module is added.
Definition: mod_mhd_phys.t:40
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,...
subroutine mhd_handle_small_values(primitive, w, x, ixIL, ixOL, subname)
Definition: mod_mhd_phys.t:936
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
Definition: mod_mhd_phys.t:175
double precision, public mhd_etah
TODO: what is this?
Definition: mod_mhd_phys.t:133
double precision, public mhd_eta_ambi
The MHD ambipolar coefficient.
Definition: mod_mhd_phys.t:136
subroutine, public mhd_get_v(w, x, ixIL, ixOL, v)
Calculate v vector.
subroutine mhd_handle_small_ei(w, x, ixIL, ixOL, ie, subname)
handle small or negative internal energy
subroutine, public mhd_phys_init()
Definition: mod_mhd_phys.t:300
logical, public, protected mhd_trac
Whether TRAC method is used.
Definition: mod_mhd_phys.t:49
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
Definition: mod_mhd_phys.t:118
procedure(mask_subroutine), pointer, public usr_mask_ambipolar
Definition: mod_mhd_phys.t:237
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_mhd_phys.t:92
subroutine, public mhd_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate v component.
subroutine, public multiplyambicoef(ixIL, ixOL, res, w, x)
multiply res by the ambipolar coefficient The ambipolar coefficient is calculated as -mhd_eta_ambi/rh...
logical, public, protected b0field_forcefree
B0 field is force-free.
Definition: mod_mhd_phys.t:181
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:178
subroutine mhd_write_info(fh)
Write this module's parameters to a snapsoht.
Definition: mod_mhd_phys.t:270
integer, public, protected tweight_
Definition: mod_mhd_phys.t:115
subroutine mhd_physical_units()
Definition: mod_mhd_phys.t:697
logical, public, protected mhd_ambipolar_sts
Whether Ambipolar term is implemented using supertimestepping.
Definition: mod_mhd_phys.t:34
subroutine mhd_e_to_ei(ixIL, ixOL, w, x)
Transform total energy to internal energy.
Definition: mod_mhd_phys.t:854
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixIL, ixOL)
Compute 2 times total magnetic energy.
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition: mod_mhd_phys.t:98
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition: mod_mhd_phys.t:172
logical, public, protected mhd_4th_order
MHD fourth order.
Definition: mod_mhd_phys.t:86
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
Definition: mod_mhd_phys.t:114
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.
subroutine get_lorentz(ixIL, ixOL, w, JxB)
Compute the Lorentz force (JxB)
integer, public, protected psi_
Indices of the GLM psi.
Definition: mod_mhd_phys.t:107
subroutine mhd_boundary_adjust(igrid, psb)
subroutine e_to_rhos(ixIL, ixOL, w, x)
Convert energy to entropy.
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
Definition: mod_particles.t:2
subroutine particles_init()
Initialize particle data and parameters.
Definition: mod_particles.t:15
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_get_a2max), pointer phys_get_a2max
Definition: mod_physics.t:68
procedure(sub_get_ct_velocity), pointer phys_get_ct_velocity
Definition: mod_physics.t:87
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:63
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
Definition: mod_physics.t:43
procedure(sub_small_values), pointer phys_handle_small_values
Definition: mod_physics.t:86
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:84
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:72
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
Definition: mod_physics.t:22
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:71
logical phys_solve_eaux
Solve internal energy and total energy equations.
Definition: mod_physics.t:34
procedure(sub_convert), pointer phys_e_to_ei
Definition: mod_physics.t:65
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:75
logical phys_total_energy
Solve total energy equation or not.
Definition: mod_physics.t:28