22 double precision,
public ::
mhd_eta = 0.0d0
43 double precision,
public,
protected ::
h_ion_fr=1d0
46 double precision,
public,
protected ::
he_ion_fr=1d0
53 double precision,
public,
protected ::
rr=1d0
55 double precision :: inv_squared_c0=0.d0, inv_squared_c=0.d0
62 integer,
public,
protected ::
rho_
64 integer,
allocatable,
public,
protected ::
mom(:)
66 integer,
public,
protected :: ^
c&m^C_
68 integer,
public,
protected ::
e_
70 integer,
public,
protected :: ^
c&b^C_
72 integer,
public,
protected ::
p_
74 integer,
public,
protected ::
ne_
76 integer,
public,
protected ::
qpar_
80 integer,
public,
protected ::
psi_
82 integer,
public,
protected ::
r_e
84 integer,
public,
protected ::
te_
86 integer,
public,
protected ::
fip_ = -1
88 logical,
public,
protected ::
mhd_fip = .false.
93 integer,
allocatable,
public,
protected ::
tracer(:)
101 integer,
parameter :: divb_none = 0
102 integer,
parameter :: divb_multigrid = -1
103 integer,
parameter :: divb_glm = 1
104 integer,
parameter :: divb_powel = 2
105 integer,
parameter :: divb_janhunen = 3
106 integer,
parameter :: divb_linde = 4
107 integer,
parameter :: divb_lindejanhunen = 5
108 integer,
parameter :: divb_lindepowel = 6
109 integer,
parameter :: divb_lindeglm = 7
110 integer,
parameter :: divb_ct = 8
151 logical,
public,
protected ::
mhd_glm = .false.
191 logical :: total_energy = .true.
195 logical :: gravity_energy
197 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
199 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
201 character(len=std_len) :: typedivbdiff =
'all'
212 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
214 integer,
intent(in) :: ixi^
l, ixo^
l
215 double precision,
intent(in) :: x(ixi^s,1:
ndim)
216 double precision,
intent(in) :: w(ixi^s,1:nw)
217 double precision,
intent(inout) :: res(ixi^s)
218 end subroutine mask_subroutine
261 subroutine mhd_read_params(files)
264 character(len=*),
intent(in) :: files(:)
283 do n = 1,
size(files)
284 open(
unitpar, file=trim(files(n)), status=
"old")
285 read(
unitpar, mhd_list,
end=111)
291 end subroutine mhd_read_params
294 subroutine mhd_write_info(fh)
296 integer,
intent(in) :: fh
299 integer,
parameter :: n_par = 1
300 double precision :: values(n_par)
301 integer,
dimension(MPI_STATUS_SIZE) :: st
302 character(len=name_len) :: names(n_par)
304 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
307 values(1) = eos%gamma
308 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
309 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
310 end subroutine mhd_write_info
338 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
342 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_internal_e=T'
349 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_hydrodynamic_e=T'
353 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_hydrodynamic_e=T'
357 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_hydrodynamic_e=T'
364 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_semirelativistic=T'
368 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_semirelativistic=T'
372 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_semirelativistic=T'
379 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
383 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
387 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
391 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_tc=F when mhd_energy=F'
395 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
399 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
403 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
407 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_energy=F'
413 if(
mype==0)
write(*,*)
'WARNING: set either parabolic TC or hyperbolic TC to F'
414 if(
mype==0)
write(*,*)
'WARNING: defaulting to only mhd_hyperbolic_tc=T'
418 call mpistop(
"mhd_hyperbolic_tc_use_perp is not supported in 1D")
428 phys_gamma = eos%gamma
441 phys_total_energy=total_energy
444 gravity_energy=.false.
446 gravity_energy=.true.
449 gravity_energy=.false.
455 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
460 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
468 type_divb = divb_none
471 if(
mhd_radiation_fld)
call mpistop(
'To verify whether mg usage for FLD versus divB can be combined')
472 type_divb = divb_multigrid
474 mg%operator_type = mg_laplacian
481 case (
'powel',
'powell')
482 type_divb = divb_powel
484 type_divb = divb_janhunen
486 type_divb = divb_linde
487 case (
'lindejanhunen')
488 type_divb = divb_lindejanhunen
490 type_divb = divb_lindepowel
494 type_divb = divb_lindeglm
499 call mpistop(
'Unknown divB fix')
502 allocate(start_indices(number_species),stop_indices(number_species))
509 mom(:) = var_set_momentum(
ndir)
515 e_ = var_set_energy()
524 mag(:) = var_set_bfield(
ndir)
528 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
534 qpar_ = var_set_fluxvar(
'q',
'q', need_bc=.false.)
536 qperp_ = var_set_fluxvar(
'qperp',
'qperp', need_bc=.false.)
547 fip_ = var_set_fluxvar(
'rho_fip',
'fip', need_bc=.false.)
552 if (eos%eos_type ==
'LTE')
then
555 else if (eos%eos_type ==
'PI')
then
566 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
572 write(*,*)
'Warning: CAK force addition together with FLD radiation'
577 write(*,*)
'Warning: Optically thin cooling together with FLD radiation'
581 call mpistop(
'using FLD implies the use of an energy equation, set mhd_energy=T')
584 call mpistop(
'using FLD not yet with semirelativistic energy formalism')
587 call mpistop(
'using FLD not yet with hydrodynamic or internal energy formalism')
590 call mpistop(
'using FLD not yet with split off rho and p')
594 r_e = var_set_radiation_energy()
603 phys_implicit_update => mhd_fld_implicit_update
604 phys_evaluate_implicit => mhd_fld_evaluate_implicit
615 stop_indices(1)=nwflux
643 allocate(iw_vector(nvector))
644 iw_vector(1) =
mom(1) - 1
645 iw_vector(2) = mag(1) - 1
648 if (.not.
allocated(flux_type))
then
649 allocate(flux_type(
ndir, nwflux))
650 flux_type = flux_default
651 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
652 call mpistop(
"phys_check error: flux_type has wrong shape")
655 if(nwflux>mag(
ndir))
then
657 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
662 flux_type(:,
psi_)=flux_special
664 flux_type(idir,mag(idir))=flux_special
668 flux_type(idir,mag(idir))=flux_tvdlf
674 phys_get_dt => mhd_get_dt
677 phys_get_cmax => mhd_get_cmax_semirelati
679 phys_get_cmax => mhd_get_cmax_semirelati_noe
683 phys_get_cmax => mhd_get_cmax_origin
685 phys_get_cmax => mhd_get_cmax_origin_noe
688 phys_get_tcutoff => mhd_get_tcutoff
689 phys_get_h_speed => mhd_get_h_speed
691 phys_get_cbounds => mhd_get_cbounds_split_rho
693 phys_get_cbounds => mhd_get_cbounds_semirelati
695 phys_get_cbounds => mhd_get_cbounds
705 phys_get_flux => mhd_get_flux_hde
708 phys_get_flux => mhd_get_flux_semirelati
710 phys_get_flux => mhd_get_flux_semirelati_noe
714 phys_get_flux => mhd_get_flux_split
716 phys_get_flux => mhd_get_flux
718 phys_get_flux => mhd_get_flux_noe
723 phys_add_source_geom => mhd_add_source_geom_semirelati
725 phys_add_source_geom => mhd_add_source_geom_split
727 phys_add_source_geom => mhd_add_source_geom
729 phys_add_source => mhd_add_source
730 phys_check_params => mhd_check_params
731 phys_write_info => mhd_write_info
734 phys_handle_small_values => mhd_handle_small_values_inte
735 mhd_handle_small_values => mhd_handle_small_values_inte
736 phys_check_w => mhd_check_w_inte
738 phys_handle_small_values => mhd_handle_small_values_hde
739 mhd_handle_small_values => mhd_handle_small_values_hde
740 phys_check_w => mhd_check_w_hde
742 phys_handle_small_values => mhd_handle_small_values_semirelati
743 mhd_handle_small_values => mhd_handle_small_values_semirelati
744 phys_check_w => mhd_check_w_semirelati
746 phys_handle_small_values => mhd_handle_small_values_split
747 mhd_handle_small_values => mhd_handle_small_values_split
748 phys_check_w => mhd_check_w_split
750 phys_handle_small_values => mhd_handle_small_values_origin
751 mhd_handle_small_values => mhd_handle_small_values_origin
752 phys_check_w => mhd_check_w_origin
754 phys_handle_small_values => mhd_handle_small_values_noe
755 mhd_handle_small_values => mhd_handle_small_values_noe
756 phys_check_w => mhd_check_w_noe
762 phys_set_equi_vars => set_equi_vars_grid
765 if(type_divb==divb_glm)
then
766 phys_modify_wlr => mhd_modify_wlr
777 transverse_ghost_cells = 1
778 phys_get_ct_velocity => mhd_get_ct_velocity_average
779 phys_update_faces => mhd_update_faces_average
781 transverse_ghost_cells = 1
782 phys_get_ct_velocity => mhd_get_ct_velocity_contact
783 phys_update_faces => mhd_update_faces_contact
785 transverse_ghost_cells = 2
786 phys_get_ct_velocity => mhd_get_ct_velocity_hll
787 phys_update_faces => mhd_update_faces_hll
789 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
792 phys_modify_wlr => mhd_modify_wlr
794 phys_boundary_adjust => mhd_boundary_adjust
800 call mpistop(
'To verify whether mg usage for FLD versus divB can be combined')
805 call mhd_physical_units()
838 call mpistop(
"thermal conduction needs mhd_energy=T")
841 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
844 call mpistop(
"radiative cooling needs mhd_energy=T")
849 iw_log_nh = var_set_wextra()
855 if(
mype==0)
write(*,*)
'WARNING: turning mhd_equi_thermal=F as no splitting or total e in use'
858 if(
mype==0)
write(*,*)
'Will subtract thermal balance in TC or RC with mhd_equi_thermal=T'
861 if(
mype==0)
write(*,*)
'WARNING: turning mhd_equi_thermal=F as no TC or RC in use'
885 phys_e_to_ei => mhd_e_to_ei_hde
886 phys_ei_to_e => mhd_ei_to_e_hde
889 phys_e_to_ei => mhd_e_to_ei_semirelati
890 phys_ei_to_e => mhd_ei_to_e_semirelati
892 if (iw_log_nh > 0)
then
909 phys_e_to_ei => mhd_e_to_ei_hde
910 phys_ei_to_e => mhd_ei_to_e_hde
912 phys_e_to_ei => mhd_e_to_ei_semirelati
913 phys_ei_to_e => mhd_ei_to_e_semirelati
936 phys_te_images => mhd_te_images
942 write(*,*)
'*****Using hyperresistivity: with mhd_eta_hyper :',
mhd_eta_hyper
946 call mpistop(
"Must have B0field=F when using hyperresistivity")
950 call mpistop(
"Must have mhd_eta_hyper positive when using hyperresistivity")
967 call mpistop(
"Must have has_equi_rho_and_p=F when mhd_rotating_frame=T")
981 call mpistop(
"Must have mhd_hall=F when mhd_semirelativistic=T")
985 call mpistop(
"Must have Cartesian coordinates for Hall")
989 phys_wider_stencil = 1
996 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
1007 phys_wider_stencil = 1
1017 call mpistop(
"CAK implementation not available in internal or semirelativistic variants")
1020 call mpistop(
"CAK force implementation not available for split off pressure and density")
1028 subroutine mhd_te_images
1033 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
1035 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
1037 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
1039 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
1042 call mpistop(
"Error in synthesize emission: Unknown convert_type")
1044 end subroutine mhd_te_images
1050 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
1054 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
1055 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1056 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
1057 double precision,
intent(in) :: my_dt
1058 logical,
intent(in) :: fix_conserve_at_step
1060 end subroutine mhd_sts_set_source_tc_mhd
1062 subroutine mhd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
1066 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
1067 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1068 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
1069 double precision,
intent(in) :: my_dt
1070 logical,
intent(in) :: fix_conserve_at_step
1072 end subroutine mhd_sts_set_source_tc_hd
1074 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1081 integer,
intent(in) :: ixi^
l, ixo^
l
1082 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1083 double precision,
intent(in) :: w(ixi^s,1:nw)
1084 double precision :: dtnew
1087 end function mhd_get_tc_dt_mhd
1089 function mhd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1096 integer,
intent(in) :: ixi^
l, ixo^
l
1097 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1098 double precision,
intent(in) :: w(ixi^s,1:nw)
1099 double precision :: dtnew
1102 end function mhd_get_tc_dt_hd
1104 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
1107 integer,
intent(in) :: ixi^
l,ixo^
l
1108 double precision,
intent(inout) :: w(ixi^s,1:nw)
1109 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1110 integer,
intent(in) :: step
1111 character(len=140) :: error_msg
1117 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1118 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1119 end subroutine mhd_tc_handle_small_e
1122 subroutine tc_params_read_mhd(fl)
1124 type(tc_fluid),
intent(inout) :: fl
1126 double precision :: tc_k_para=0d0
1127 double precision :: tc_k_perp=0d0
1130 logical :: tc_perpendicular=.false.
1131 logical :: tc_saturate=.false.
1132 logical :: tc_patch_eint=.false.
1133 double precision :: trac_t_floor=0.d0
1134 character(len=std_len) :: tc_slope_limiter=
"MC"
1136 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp, tc_patch_eint, trac_t_floor
1140 read(
unitpar, tc_list,
end=111)
1144 fl%tc_perpendicular = tc_perpendicular
1145 fl%tc_saturate = tc_saturate
1146 fl%tc_patch_eint = tc_patch_eint
1147 fl%tc_k_para = tc_k_para
1148 fl%tc_k_perp = tc_k_perp
1149 fl%trac_T_floor = trac_t_floor / unit_temperature
1150 select case(tc_slope_limiter)
1152 fl%tc_slope_limiter = 0
1155 fl%tc_slope_limiter = 1
1158 fl%tc_slope_limiter = 2
1161 fl%tc_slope_limiter = 3
1164 fl%tc_slope_limiter = 4
1167 fl%tc_slope_limiter = 5
1169 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod, superbee, koren, vanleer")
1171 end subroutine tc_params_read_mhd
1175 subroutine rc_params_read(fl)
1178 type(rc_fluid),
intent(inout) :: fl
1182 double precision :: rad_damp_height=0.5d0
1183 double precision :: rad_damp_scale=0.15d0
1186 integer :: ncool = 4000
1188 logical :: tfix=.false.
1190 logical :: rc_split=.false.
1191 logical :: rad_damp=.false.
1193 character(len=std_len) :: coolcurve=
'JCcorona'
1194 logical :: rad_newton = .false.
1195 double precision :: rad_newton_trad = 0.006d0
1196 double precision :: rad_newton_rhosurf = 1.d4
1197 double precision :: rad_newton_pthick = 25.d0
1199 double precision :: cfrac=0.1d0
1200 double precision :: rad_cut_hgt=0.5d0
1201 double precision :: rad_cut_dey=0.15d0
1203 character(len=8) :: rc_y_mod_quadrature=
'boole'
1204 integer :: rc_y_mod_n_sub=16
1206 namelist /rc_list/ coolcurve, ncool, cfrac, tlow, tfix, rc_split, &
1207 rad_cut_hgt, rad_cut_dey, &
1208 rc_y_mod_quadrature, rc_y_mod_n_sub, &
1209 rad_newton, rad_newton_trad, rad_newton_rhosurf, &
1210 rad_newton_pthick, rad_damp, rad_damp_height, rad_damp_scale
1214 read(
unitpar, rc_list,
end=111)
1219 fl%coolcurve=coolcurve
1222 fl%rc_split=rc_split
1224 fl%rad_cut_hgt=rad_cut_hgt
1225 fl%rad_cut_dey=rad_cut_dey
1226 fl%Y_mod_quadrature=rc_y_mod_quadrature
1227 fl%Y_mod_N_sub=rc_y_mod_n_sub
1228 fl%rad_damp=rad_damp
1229 fl%rad_damp_height=rad_damp_height
1230 fl%rad_damp_scale=rad_damp_scale
1231 fl%rad_newton=rad_newton
1232 fl%rad_newton_trad=rad_newton_trad
1233 fl%rad_newton_rhosurf=rad_newton_rhosurf
1234 fl%rad_newton_pthick=rad_newton_pthick
1235 end subroutine rc_params_read
1238 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1241 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1242 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1244 double precision :: delx(ixi^s,1:
ndim)
1245 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1246 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1252 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1256 hxo^
l=ixo^
l-
kr(idims,^
d);
1262 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1265 xshift^
d=half*(one-
kr(^
d,idims));
1272 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1276 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1278 end subroutine set_equi_vars_grid_faces
1281 subroutine set_equi_vars_grid(igrid)
1285 integer,
intent(in) :: igrid
1291 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1293 end subroutine set_equi_vars_grid
1296 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1298 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1299 double precision,
intent(in) :: w(ixi^s, 1:nw)
1300 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1301 double precision :: wnew(ixo^s, 1:nwc)
1308 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1314 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1318 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1320 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+
block%equi_vars(ixo^s,
equi_pe0_,0)*eos%inv_gamma_minus_1
1322 if(
b0field .and. total_energy)
then
1323 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1324 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1328 end function convert_vars_splitting
1330 subroutine mhd_check_params
1337 ngridvars,num_particles,physics_type_particles
1340 double precision :: a,
b,xfrac,yfrac
1345 if (particles_eta < zero) particles_eta =
mhd_eta
1346 if (particles_etah < zero) particles_eta =
mhd_etah
1351 if (eos%gamma <= 0.0d0)
call mpistop (
"Error: gamma <= 0")
1352 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1355 if (eos%gamma <= 0.0d0 .or. eos%gamma == 1.0d0) &
1356 call mpistop (
"Error: gamma <= 0 or gamma == 1")
1362 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1367 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1375 call mpistop(
'select IMEX scheme for FLD radiation use')
1378 call phys_set_mg_bounds()
1380 if(.not.
fld_no_mg)
call mpistop(
'multigrid must have BCs for IMEX and FLD radiation use')
1383 write(*,*)
'==FLD SETUP======================'
1384 write(*,*)
'Using FLD with settings:'
1389 write(*,*)
'Using FLD with settings: fld_kappa0=',
fld_kappa0
1390 write(*,*)
'Using FLD with settings: fld_opal_table=',
fld_opal_table
1392 write(*,*)
'Using FLD with settings: fld_bisect_tol=',
fld_bisect_tol
1393 write(*,*)
'Using FLD with settings: fld_diff_tol=',
fld_diff_tol
1397 print *,
'NORMALIZED arad_norm=',
arad_norm
1398 print *,
'NORMALIZED c_norm=',
c_norm
1405 print *,
'physical fld_kappa (in cgs or SI) =',
fld_kappa0
1408 write(*,*)
'===FLD SETUP====================='
1413 write(*,*)
'====MHD run with settings===================='
1414 write(*,*)
'Using mod_mhd_phys with settings:'
1416 write(*,*)
'Dimensionality :',
ndim
1417 write(*,*)
'vector components:',
ndir
1419 write(*,*)
'number of variables nw=',nw
1420 write(*,*)
' start index iwstart=',iwstart
1421 write(*,*)
'number of vector variables=',nvector
1422 write(*,*)
'number of stagger variables nws=',nws
1423 write(*,*)
'number of variables with BCs=',nwgc
1424 write(*,*)
'number of vars with fluxes=',nwflux
1425 write(*,*)
'number of vars with flux + BC=',nwfluxbc
1426 write(*,*)
'number of auxiliary variables=',nwaux
1427 write(*,*)
'number of extra vars without flux=',nwextra
1428 write(*,*)
'number of extra vars for wextra=',nw_extra
1429 write(*,*)
'number of auxiliary I/O variables=',
nwauxio
1431 write(*,*)
' mhd_energy=',
mhd_energy,
' with total_energy=',total_energy
1436 write(*,*)
' mhd_eta=',
mhd_eta,
' nonzero implies resistivity'
1450 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
1451 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
1452 write(*,*)
'*****Using particles: npayload,ngridvars :', npayload,ngridvars
1453 write(*,*)
'*****Using particles: nusrpayload :', nusrpayload
1454 write(*,*)
'*****Using particles: num_particles :', num_particles
1455 write(*,*)
'*****Using particles: physics_type_particles=',physics_type_particles
1458 write(*,*)
'number due to phys_wider_stencil=',phys_wider_stencil
1459 write(*,*)
'==========================================='
1460 print *,
'========EOS and UNITS==========='
1462 print *,
'gamma=',eos%gamma
1463 print *,
'He_abundance =',eos%He_abundance
1465 print *,
'========EOS and UNITS==========='
1487 print *,
' compare this to ',mp_si*(1.d0+4.d0*eos%He_abundance)
1489 print *,
' compare this to ',mp_cgs*(1.d0+4.d0*eos%He_abundance)
1493 print *,
' compare this to ',kb_si*(2.d0+3.d0*eos%He_abundance)
1497 print *,
' compare this to ',kb_cgs*(2.d0+3.d0*eos%He_abundance)
1501 if(eos%eos_type /=
'LTE')
then
1502 print *,
'mean molecular weight mu is =',a/
b,
' = ', (1.d0+4.d0*eos%He_abundance)/(2.d0+3.d0*eos%He_abundance)
1504 yfrac=4.d0*eos%He_abundance/(1.d0+4.d0*eos%He_abundance)
1505 print *,
'mass fraction hydrogen X is =',1/a,
' and this equals ', 1.d0/(1.d0+4.d0*eos%He_abundance)
1506 print *,
'mass fraction helium Y is =',yfrac
1507 print *,
' check that 1/mu',
b/a,
' is equal to 2X+3Y/4=',2.d0*xfrac+3.d0*yfrac/4.d0
1508 print *,
' ratio n_e/n_p=',1.d0+2.0d0*eos%He_abundance
1510 print *,
'========UNITS==========='
1513 end subroutine mhd_check_params
1515 subroutine mhd_physical_units()
1517 double precision :: mp,kb,miu0,c_lightspeed,xfrac,sigma_telectron
1518 double precision :: a,
b
1526 sigma_telectron=sigma_te_si
1532 c_lightspeed=const_c
1533 sigma_telectron=sigma_te_cgs
1538 if (eos%eos_type ==
'LTE')
then
1542 eos%nH2rhoFactor = 1d0+4d0*eos%He_abundance
1543 rr=(2d0+3d0*eos%He_abundance) / (1d0+4d0*eos%He_abundance)
1544 xfrac=1.d0/(1.d0+4.d0*eos%He_abundance)
1548 a=1d0+4d0*eos%He_abundance
1549 if(eos%eos_type==
'PI')
then
1552 b=2d0+3d0*eos%He_abundance
1684 eos%inv_squared_c0 = inv_squared_c0
1685 eos%inv_squared_c = inv_squared_c
1701 end subroutine mhd_physical_units
1703 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1706 logical,
intent(in) :: primitive
1707 logical,
intent(inout) :: flag(ixi^s,1:nw)
1708 integer,
intent(in) :: ixi^
l, ixo^
l
1709 double precision,
intent(in) :: w(ixi^s,nw)
1711 double precision :: tmp,
b(1:
ndir),v(1:
ndir),factor
1722 {
do ix^db=ixomin^db,ixomax^db \}
1726 {
do ix^db=ixomin^db,ixomax^db \}
1728 tmp=(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)*inv_squared_c
1729 factor=1.0d0/(w(ix^d,
rho_)*(w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)*inv_squared_c))
1730 ^
c&v(^
c)=factor*(w(ix^d,
m^
c_)*w(ix^d,
rho_)+w(ix^d,
b^
c_)*tmp)\
1733 b(1)=w(ix^d,b2_)*v(3)-w(ix^d,b3_)*v(2)
1734 b(2)=w(ix^d,b3_)*v(1)-w(ix^d,b1_)*v(3)
1735 b(3)=w(ix^d,b1_)*v(2)-w(ix^d,b2_)*v(1)
1740 b(2)=w(ix^d,b1_)*v(2)-w(ix^d,b2_)*v(1)
1746 tmp=w(ix^d,
e_)-half*((^
c&v(^
c)**2+)*w(ix^d,
rho_)&
1747 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(^
c)**2+)*inv_squared_c)
1748 if(tmp<small_e) flag(ix^d,
e_)=.true.
1754 end subroutine mhd_check_w_semirelati
1756 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1759 logical,
intent(in) :: primitive
1760 integer,
intent(in) :: ixi^
l, ixo^
l
1761 double precision,
intent(in) :: w(ixi^s,nw)
1762 logical,
intent(inout) :: flag(ixi^s,1:nw)
1767 {
do ix^db=ixomin^db,ixomax^db\}
1780 end subroutine mhd_check_w_origin
1782 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1785 logical,
intent(in) :: primitive
1786 integer,
intent(in) :: ixi^
l, ixo^
l
1787 double precision,
intent(in) :: w(ixi^s,nw)
1788 logical,
intent(inout) :: flag(ixi^s,1:nw)
1790 double precision :: tmp
1794 {
do ix^db=ixomin^db,ixomax^db\}
1800 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1805 end subroutine mhd_check_w_split
1807 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1810 logical,
intent(in) :: primitive
1811 integer,
intent(in) :: ixi^
l, ixo^
l
1812 double precision,
intent(in) :: w(ixi^s,nw)
1813 logical,
intent(inout) :: flag(ixi^s,1:nw)
1818 {
do ix^db=ixomin^db,ixomax^db\}
1822 end subroutine mhd_check_w_noe
1824 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1827 logical,
intent(in) :: primitive
1828 integer,
intent(in) :: ixi^
l, ixo^
l
1829 double precision,
intent(in) :: w(ixi^s,nw)
1830 logical,
intent(inout) :: flag(ixi^s,1:nw)
1835 {
do ix^db=ixomin^db,ixomax^db\}
1844 end subroutine mhd_check_w_inte
1846 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1849 logical,
intent(in) :: primitive
1850 integer,
intent(in) :: ixi^
l, ixo^
l
1851 double precision,
intent(in) :: w(ixi^s,nw)
1852 logical,
intent(inout) :: flag(ixi^s,1:nw)
1857 {
do ix^db=ixomin^db,ixomax^db\}
1866 end subroutine mhd_check_w_hde
1868 subroutine mhd_bound_fip(primitive, ixI^L, ixO^L, w)
1870 logical,
intent(in) :: primitive
1871 integer,
intent(in) :: ixi^
l, ixo^
l
1872 double precision,
intent(inout) :: w(ixi^s,1:nw)
1874 double precision :: rho_safe(ixi^s), fip_prim(ixi^s)
1886 fip_prim(ixo^s) = w(ixo^s,
fip_) / rho_safe(ixo^s)
1887 fip_prim(ixo^s) = min(
maxfip, max(
minfip, fip_prim(ixo^s)))
1888 w(ixo^s,
fip_) = rho_safe(ixo^s) * fip_prim(ixo^s)
1890 end subroutine mhd_bound_fip
1895 integer,
intent(in) :: ixi^
l, ixo^
l
1896 double precision,
intent(inout) :: w(ixi^s, nw)
1897 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1902 {
do ix^db=ixomin^db,ixomax^db\}
1905 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1907 +(^
c&w(ix^
d,
b^
c_)**2+))
1910 {
do ix^db=ixomin^db,ixomax^db\}
1912 w(ix^d,
e_)=w(ix^d,
e_)&
1913 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1914 +(^
c&w(ix^d,
b^
c_)**2+))
1920 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
1922 integer,
intent(in) :: ixi^
l, ixo^
l
1923 double precision,
intent(inout) :: w(ixi^s, nw)
1924 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1928 {
do ix^db=ixomin^db,ixomax^db\}
1934 end subroutine mhd_ei_to_e_hde
1937 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
1939 integer,
intent(in) :: ixi^
l, ixo^
l
1940 double precision,
intent(inout) :: w(ixi^s, nw)
1941 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1943 w(ixo^s,
p_)=w(ixo^s,
e_)*eos%gamma_minus_1
1945 call eos%to_conserved(ixi^
l,ixo^
l,w,x)
1947 end subroutine mhd_ei_to_e_semirelati
1952 integer,
intent(in) :: ixi^
l, ixo^
l
1953 double precision,
intent(inout) :: w(ixi^s, nw)
1954 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1959 {
do ix^db=ixomin^db,ixomax^db\}
1962 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1964 +(^
c&w(ix^
d,
b^
c_)**2+))
1967 {
do ix^db=ixomin^db,ixomax^db\}
1969 w(ix^d,
e_)=w(ix^d,
e_)&
1970 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1971 +(^
c&w(ix^d,
b^
c_)**2+))
1975 if(fix_small_values)
then
1976 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
1984 subroutine mhd_e_to_ei_and_cache_log_nh(ixI^L,ixO^L,w,x)
1986 integer,
intent(in) :: ixi^
l, ixo^
l
1987 double precision,
intent(inout) :: w(ixi^s, nw)
1988 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1991 block%wextra(ixo^s, iw_log_nh) = dlog10(w(ixo^s,
rho_) / eos%nH2rhoFactor)
1992 end subroutine mhd_e_to_ei_and_cache_log_nh
1995 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
1997 integer,
intent(in) :: ixi^
l, ixo^
l
1998 double precision,
intent(inout) :: w(ixi^s, nw)
1999 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2003 {
do ix^db=ixomin^db,ixomax^db\}
2009 if(fix_small_values)
then
2010 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2013 end subroutine mhd_e_to_ei_hde
2016 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2018 integer,
intent(in) :: ixi^
l, ixo^
l
2019 double precision,
intent(inout) :: w(ixi^s, nw)
2020 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2022 call eos%to_primitive(ixi^
l,ixo^
l,w,x)
2023 w(ixo^s,
e_)=w(ixo^s,
p_)*eos%inv_gamma_minus_1
2025 end subroutine mhd_e_to_ei_semirelati
2027 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2030 logical,
intent(in) :: primitive
2031 integer,
intent(in) :: ixi^
l,ixo^
l
2032 double precision,
intent(inout) :: w(ixi^s,1:nw)
2033 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2034 character(len=*),
intent(in) :: subname
2036 double precision :: e(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2037 double precision :: tmp, factor
2039 logical :: flag(ixi^s,1:nw)
2048 {
do ix^db=ixomin^db,ixomax^db\}
2050 tmp=(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)*inv_squared_c
2051 factor=1.0d0/(w(ix^
d,
rho_)*(w(ix^
d,
rho_)+(^
c&w(ix^
d,
b^
c_)**2+)*inv_squared_c))
2055 e(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2056 e(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2057 e(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2061 e(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2067 pressure(ix^
d)=eos%gamma_minus_1*(w(ix^
d,
e_)&
2068 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2069 +(^
c&w(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
2076 select case (small_values_method)
2078 {
do ix^db=ixomin^db,ixomax^db\}
2079 if(flag(ix^d,
rho_))
then
2080 w(ix^d,
rho_) = small_density
2081 ^
c&w(ix^d,
m^
c_)=0.d0\
2085 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2087 if(flag(ix^d,
e_))
then
2088 w(ix^d,
e_)=small_pressure*eos%inv_gamma_minus_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2089 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
2096 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2099 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2101 w(ixo^s,
e_)=pressure(ixo^s)
2102 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2103 {
do ix^db=ixomin^db,ixomax^db\}
2104 w(ix^d,
e_)=w(ix^d,
p_)*eos%inv_gamma_minus_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2105 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
2110 if(.not.primitive)
then
2112 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2115 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2118 if (
mhd_fip)
call mhd_bound_fip(primitive, ixi^l, ixo^l, w)
2119 end subroutine mhd_handle_small_values_semirelati
2121 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2124 logical,
intent(in) :: primitive
2125 integer,
intent(in) :: ixi^
l,ixo^
l
2126 double precision,
intent(inout) :: w(ixi^s,1:nw)
2127 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2128 character(len=*),
intent(in) :: subname
2131 logical :: flag(ixi^s,1:nw)
2133 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2138 {
do ix^db=ixomin^db,ixomax^db\}
2142 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2159 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2161 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2164 {
do ix^db=iximin^db,iximax^db\}
2165 w(ix^d,
e_)=w(ix^d,
e_)&
2166 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2168 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2170 {
do ix^db=iximin^db,iximax^db\}
2171 w(ix^d,
e_)=w(ix^d,
e_)&
2172 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2176 call small_values_average(ixi^l, ixo^l, w, x, flag,
r_e)
2179 if(.not.primitive)
then
2181 {
do ix^db=ixomin^db,ixomax^db\}
2183 w(ix^d,
p_)=eos%gamma_minus_1*(w(ix^d,
e_)&
2184 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2187 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2190 if (
mhd_fip)
call mhd_bound_fip(primitive, ixi^l, ixo^l, w)
2191 end subroutine mhd_handle_small_values_origin
2193 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2196 logical,
intent(in) :: primitive
2197 integer,
intent(in) :: ixi^
l,ixo^
l
2198 double precision,
intent(inout) :: w(ixi^s,1:nw)
2199 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2200 character(len=*),
intent(in) :: subname
2202 double precision :: rho
2204 logical :: flag(ixi^s,1:nw)
2206 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2211 {
do ix^db=ixomin^db,ixomax^db\}
2216 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2229 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2231 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2234 {
do ix^db=iximin^db,iximax^db\}
2236 w(ix^d,
e_)=w(ix^d,
e_)&
2237 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2239 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2241 {
do ix^db=iximin^db,iximax^db\}
2243 w(ix^d,
e_)=w(ix^d,
e_)&
2244 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2248 if(.not.primitive)
then
2250 {
do ix^db=ixomin^db,ixomax^db\}
2252 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2253 w(ix^d,
p_)=eos%gamma_minus_1*(w(ix^d,
e_)&
2254 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2257 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2260 if (
mhd_fip)
call mhd_bound_fip(primitive, ixi^l, ixo^l, w)
2261 end subroutine mhd_handle_small_values_split
2263 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2266 logical,
intent(in) :: primitive
2267 integer,
intent(in) :: ixi^
l,ixo^
l
2268 double precision,
intent(inout) :: w(ixi^s,1:nw)
2269 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2270 character(len=*),
intent(in) :: subname
2273 logical :: flag(ixi^s,1:nw)
2275 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2280 {
do ix^db=ixomin^db,ixomax^db\}
2281 if(flag(ix^
d,
rho_))
then
2283 ^
c&w(ix^
d,
m^
c_)=0.d0\
2293 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2295 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2297 if(.not.primitive)
then
2299 {
do ix^db=ixomin^db,ixomax^db\}
2301 w(ix^d,
p_)=eos%gamma_minus_1*w(ix^d,
e_)
2304 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2307 if (
mhd_fip)
call mhd_bound_fip(primitive, ixi^l, ixo^l, w)
2308 end subroutine mhd_handle_small_values_inte
2310 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2313 logical,
intent(in) :: primitive
2314 integer,
intent(in) :: ixi^
l,ixo^
l
2315 double precision,
intent(inout) :: w(ixi^s,1:nw)
2316 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2317 character(len=*),
intent(in) :: subname
2320 logical :: flag(ixi^s,1:nw)
2322 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2327 {
do ix^db=ixomin^db,ixomax^db\}
2331 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2337 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2339 if(.not.primitive)
then
2341 {
do ix^db=ixomin^db,ixomax^db\}
2345 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2348 if (
mhd_fip)
call mhd_bound_fip(primitive, ixi^l, ixo^l, w)
2349 end subroutine mhd_handle_small_values_noe
2351 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2354 logical,
intent(in) :: primitive
2355 integer,
intent(in) :: ixi^
l,ixo^
l
2356 double precision,
intent(inout) :: w(ixi^s,1:nw)
2357 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2358 character(len=*),
intent(in) :: subname
2361 logical :: flag(ixi^s,1:nw)
2363 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2368 {
do ix^db=ixomin^db,ixomax^db\}
2369 if(flag(ix^
d,
rho_))
then
2371 ^
c&w(ix^
d,
m^
c_)=0.d0\
2381 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2383 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2385 if(.not.primitive)
then
2387 {
do ix^db=ixomin^db,ixomax^db\}
2389 w(ix^d,
p_)=eos%gamma_minus_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2392 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2395 if (
mhd_fip)
call mhd_bound_fip(primitive, ixi^l, ixo^l, w)
2396 end subroutine mhd_handle_small_values_hde
2402 integer,
intent(in) :: ixi^
l, ixo^
l
2403 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2404 double precision,
intent(out) :: v(ixi^s,
ndir)
2406 double precision :: rho(ixi^s)
2411 rho(ixo^s)=1.d0/rho(ixo^s)
2414 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2420 subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,cs2)
2423 integer,
intent(in) :: ixi^
l, ixo^
l
2424 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2425 double precision,
intent(inout) :: cs2(ixi^s)
2427 double precision :: rho, inv_rho, ploc
2430 {
do ix^db=ixomin^db,ixomax^db \}
2440 cs2(ix^
d)=eos%gamma*ploc*inv_rho
2442 end subroutine mhd_get_csound2
2445 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2448 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2449 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2450 double precision,
intent(inout) :: cmax(ixi^s)
2452 double precision :: rho, inv_rho, ploc, cfast2, avmincs2, b2, kmax
2453 double precision :: cs2(ixi^s)
2454 double precision,
allocatable :: w_eos(:^
d&,:)
2464 allocate(w_eos(ixi^s,nw))
2465 w_eos(ixo^s,:) = w(ixo^s,:)
2470 call eos%get_csound2(w_eos, x, ixi^
l, ixo^
l, cs2)
2473 call eos%get_csound2(w, x, ixi^
l, ixo^
l, cs2)
2477 {
do ix^db=ixomin^db,ixomax^db \}
2486 cmax(ix^
d)=cs2(ix^
d)
2489 cfast2=b2*inv_rho+cmax(ix^
d)
2490 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2491 if(avmincs2<zero) avmincs2=zero
2492 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2496 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2498 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2501 {
do ix^db=ixomin^db,ixomax^db \}
2504 ploc=(w(ix^d,
p_)+block%equi_vars(ix^d,
equi_pe0_,b0i))
2510 cmax(ix^d)=cs2(ix^d)
2512 b2=(^
c&w(ix^d,
b^
c_)**2+)
2513 cfast2=b2*inv_rho+cmax(ix^d)
2514 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2515 if(avmincs2<zero) avmincs2=zero
2516 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2520 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2522 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2526 end subroutine mhd_get_cmax_origin
2529 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2533 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2534 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2535 double precision,
intent(inout) :: cmax(ixi^s)
2537 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2538 double precision :: adiabs(ixi^s), gammas(ixi^s)
2553 {
do ix^db=ixomin^db,ixomax^db \}
2557 cmax(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*rho**(gammas(ix^
d)-1.d0)
2559 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2560 cfast2=b2*inv_rho+cmax(ix^
d)
2561 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*w(ix^
d,mag(idim))**2*inv_rho
2562 if(avmincs2<zero) avmincs2=zero
2563 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2567 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2569 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2572 end subroutine mhd_get_cmax_origin_noe
2575 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2578 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2579 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2580 double precision,
intent(inout):: cmax(ixi^s)
2582 double precision :: csound, avmincs2, idim_alfven_speed2
2583 double precision :: inv_rho, alfven_speed2, gamma2
2586 {
do ix^db=ixomin^db,ixomax^db \}
2587 inv_rho=1.d0/w(ix^
d,
rho_)
2588 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2589 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2590 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2592 csound=eos%gamma*w(ix^
d,
p_)*inv_rho
2593 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2596 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2597 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2598 if(avmincs2<zero) avmincs2=zero
2600 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2601 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2604 end subroutine mhd_get_cmax_semirelati
2607 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2611 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2612 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2613 double precision,
intent(inout):: cmax(ixi^s)
2615 double precision :: adiabs(ixi^s), gammas(ixi^s)
2616 double precision :: csound, avmincs2, idim_alfven_speed2
2617 double precision :: inv_rho, alfven_speed2, gamma2
2631 {
do ix^db=ixomin^db,ixomax^db \}
2632 inv_rho=1.d0/w(ix^
d,
rho_)
2633 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2634 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2635 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2636 csound=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
2637 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2640 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2641 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2642 if(avmincs2<zero) avmincs2=zero
2644 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2645 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2648 end subroutine mhd_get_cmax_semirelati_noe
2651 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2654 integer,
intent(in) :: ixi^
l,ixo^
l
2655 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2657 double precision,
intent(inout) :: w(ixi^s,1:nw)
2658 double precision,
intent(out) :: tco_local,tmax_local
2660 double precision,
parameter :: trac_delta=0.25d0
2661 double precision :: te(ixi^s),lts(ixi^s)
2662 double precision,
dimension(1:ndim) :: bdir, bunitvec
2663 double precision,
dimension(ixI^S,1:ndim) :: gradt
2664 double precision :: ltrc,ltrp,altr
2665 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2666 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2668 if (eos%eos_type ==
'LTE' .or. eos%eos_type ==
'PI')
then
2669 te(ixi^s) = w(ixi^s,
te_)
2671 call eos%get_Rfactor(w,x,ixi^
l,ixi^
l,te)
2672 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2675 tmax_local=maxval(te(ixo^s))
2683 do ix1=ixomin1,ixomax1
2684 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2685 if(lts(ix1)>trac_delta)
then
2686 tco_local=max(tco_local,te(ix1))
2698 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2699 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2700 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2701 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2703 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2714 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2721 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2726 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2727 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2731 if(bdir(1)/=0.d0)
then
2732 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2734 block%special_values(3)=0.d0
2736 if(bdir(2)/=0.d0)
then
2737 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2739 block%special_values(4)=0.d0
2743 if(bdir(1)/=0.d0)
then
2744 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2745 (bdir(3)/bdir(1))**2)
2747 block%special_values(3)=0.d0
2749 if(bdir(2)/=0.d0)
then
2750 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2751 (bdir(3)/bdir(2))**2)
2753 block%special_values(4)=0.d0
2755 if(bdir(3)/=0.d0)
then
2756 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2757 (bdir(2)/bdir(3))**2)
2759 block%special_values(5)=0.d0
2764 block%special_values(1)=zero
2765 {
do ix^db=ixomin^db,ixomax^db\}
2767 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2769 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2772 if(bdir(1)/=0.d0)
then
2773 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2777 if(bdir(2)/=0.d0)
then
2778 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2783 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2784 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2787 if(bdir(1)/=0.d0)
then
2788 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2792 if(bdir(2)/=0.d0)
then
2793 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2797 if(bdir(3)/=0.d0)
then
2798 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2803 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2804 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2806 if(lts(ix^d)>trac_delta)
then
2807 block%special_values(1)=max(block%special_values(1),te(ix^d))
2810 block%special_values(2)=tmax_local
2829 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2830 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2831 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2835 {
do ix^db=ixpmin^db,ixpmax^db\}
2836 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2838 if(bdir(1)/=0.d0)
then
2839 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2843 if(bdir(2)/=0.d0)
then
2844 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2850 if(bdir(1)/=0.d0)
then
2851 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2855 if(bdir(2)/=0.d0)
then
2856 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2860 if(bdir(3)/=0.d0)
then
2861 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2867 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2869 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2870 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2873 {
do ix^db=ixpmin^db,ixpmax^db\}
2875 if(w(ix^d,iw_mag(1))/=0.d0)
then
2876 bunitvec(1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2)
2880 if(w(ix^d,iw_mag(2))/=0.d0)
then
2881 bunitvec(2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2)
2887 if(w(ix^d,iw_mag(1))/=0.d0)
then
2888 bunitvec(1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2+&
2889 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
2893 if(w(ix^d,iw_mag(2))/=0.d0)
then
2894 bunitvec(2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2+&
2895 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
2899 if(w(ix^d,iw_mag(3))/=0.d0)
then
2900 bunitvec(3)=sign(1.d0,w(ix^d,iw_mag(3)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(3)))**2+&
2901 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
2907 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2909 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2910 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2916 {
do ix^db=ixpmin^db,ixpmax^db\}
2918 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
2919 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
2920 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2923 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
2924 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
2925 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
2926 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2932 call mpistop(
"unknown mhd_trac_type")
2935 end subroutine mhd_get_tcutoff
2938 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2941 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2942 double precision,
intent(in) :: wprim(ixi^s, nw)
2943 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2944 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2946 double precision :: csound(ixi^s,
ndim)
2947 double precision,
allocatable :: tmp(:^
d&)
2948 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2952 allocate(tmp(ixa^s))
2955 call mhd_get_csound_prim_split(wprim,x,ixi^
l,ixa^
l,id,tmp)
2957 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2959 csound(ixa^s,id)=tmp(ixa^s)
2962 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2963 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2964 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2965 hspeed(ixc^s,1)=0.5d0*abs(wprim(jxc^s,
mom(idim))+csound(jxc^s,idim)-wprim(ixc^s,
mom(idim))+csound(ixc^s,idim))
2969 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2970 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2971 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,
mom(id))+csound(ixa^s,id)-wprim(ixc^s,
mom(id))+csound(ixc^s,id)))
2972 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2973 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2974 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixc^s,
mom(id))+csound(ixc^s,id)-wprim(ixa^s,
mom(id))+csound(ixa^s,id)))
2979 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2980 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2981 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,
mom(id))+csound(ixa^s,id)-wprim(jxc^s,
mom(id))+csound(jxc^s,id)))
2982 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2983 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2984 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(jxc^s,
mom(id))+csound(jxc^s,id)-wprim(ixa^s,
mom(id))+csound(ixa^s,id)))
2988 end subroutine mhd_get_h_speed
2991 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2994 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2995 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2996 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2997 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2998 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
2999 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3000 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3002 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3003 double precision :: umean, dmean, tmp1, tmp2, tmp3
3010 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3011 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3012 if(
present(cmin))
then
3013 {
do ix^db=ixomin^db,ixomax^db\}
3014 tmp1=sqrt(wlp(ix^
d,
rho_))
3015 tmp2=sqrt(wrp(ix^
d,
rho_))
3016 tmp3=1.d0/(tmp1+tmp2)
3017 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3018 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3019 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3020 cmin(ix^
d,1)=umean-dmean
3021 cmax(ix^
d,1)=umean+dmean
3023 if(h_correction)
then
3024 {
do ix^db=ixomin^db,ixomax^db\}
3025 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3026 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3030 {
do ix^db=ixomin^db,ixomax^db\}
3031 tmp1=sqrt(wlp(ix^d,
rho_))
3032 tmp2=sqrt(wrp(ix^d,
rho_))
3033 tmp3=1.d0/(tmp1+tmp2)
3034 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3035 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3036 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3037 cmax(ix^d,1)=abs(umean)+dmean
3041 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3042 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3043 if(
present(cmin))
then
3044 {
do ix^db=ixomin^db,ixomax^db\}
3045 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3046 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3048 if(h_correction)
then
3049 {
do ix^db=ixomin^db,ixomax^db\}
3050 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3051 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3055 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3059 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3060 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3061 if(
present(cmin))
then
3062 {
do ix^db=ixomin^db,ixomax^db\}
3063 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3064 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3065 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3067 if(h_correction)
then
3068 {
do ix^db=ixomin^db,ixomax^db\}
3069 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3070 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3074 {
do ix^db=ixomin^db,ixomax^db\}
3075 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3076 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3081 end subroutine mhd_get_cbounds
3084 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3087 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3088 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3089 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3090 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3091 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3092 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3093 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3095 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3100 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3101 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3103 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3104 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3106 if(
present(cmin))
then
3107 {
do ix^db=ixomin^db,ixomax^db\}
3108 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3109 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3110 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3113 {
do ix^db=ixomin^db,ixomax^db\}
3114 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3115 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3119 end subroutine mhd_get_cbounds_semirelati
3122 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3125 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3126 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3127 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3128 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3129 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3130 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3131 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3133 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3134 double precision :: umean, dmean, tmp1, tmp2, tmp3
3141 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3142 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3143 if(
present(cmin))
then
3144 {
do ix^db=ixomin^db,ixomax^db\}
3147 tmp3=1.d0/(tmp1+tmp2)
3148 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3149 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3150 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3151 cmin(ix^
d,1)=umean-dmean
3152 cmax(ix^
d,1)=umean+dmean
3154 if(h_correction)
then
3155 {
do ix^db=ixomin^db,ixomax^db\}
3156 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3157 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3161 {
do ix^db=ixomin^db,ixomax^db\}
3164 tmp3=1.d0/(tmp1+tmp2)
3165 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3166 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3167 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3168 cmax(ix^d,1)=abs(umean)+dmean
3172 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3173 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3174 if(
present(cmin))
then
3175 {
do ix^db=ixomin^db,ixomax^db\}
3176 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3177 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3179 if(h_correction)
then
3180 {
do ix^db=ixomin^db,ixomax^db\}
3181 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3182 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3186 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3190 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3191 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3192 if(
present(cmin))
then
3193 {
do ix^db=ixomin^db,ixomax^db\}
3194 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3195 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3196 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3198 if(h_correction)
then
3199 {
do ix^db=ixomin^db,ixomax^db\}
3200 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3201 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3205 {
do ix^db=ixomin^db,ixomax^db\}
3206 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3207 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3212 end subroutine mhd_get_cbounds_split_rho
3215 subroutine mhd_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3218 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3219 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3220 double precision,
intent(in) :: cmax(ixi^s)
3221 double precision,
intent(in),
optional :: cmin(ixi^s)
3222 type(ct_velocity),
intent(inout):: vcts
3224 end subroutine mhd_get_ct_velocity_average
3226 subroutine mhd_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3229 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3230 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3231 double precision,
intent(in) :: cmax(ixi^s)
3232 double precision,
intent(in),
optional :: cmin(ixi^s)
3233 type(ct_velocity),
intent(inout):: vcts
3235 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3237 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3239 end subroutine mhd_get_ct_velocity_contact
3241 subroutine mhd_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3244 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3245 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3246 double precision,
intent(in) :: cmax(ixi^s)
3247 double precision,
intent(in),
optional :: cmin(ixi^s)
3248 type(ct_velocity),
intent(inout):: vcts
3250 integer :: idime,idimn
3252 if(.not.
allocated(vcts%vbarC))
then
3253 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3254 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3257 if(
present(cmin))
then
3258 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3259 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3261 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3262 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3265 idimn=mod(idim,
ndir)+1
3266 idime=mod(idim+1,
ndir)+1
3268 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3269 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3270 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3271 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3272 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3274 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3275 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3276 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3277 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3278 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3280 end subroutine mhd_get_ct_velocity_hll
3287 integer,
intent(in) :: ixi^
l, ixo^
l
3288 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3289 double precision,
intent(out):: csound(ixi^s)
3291 double precision :: wprim(ixi^s, nw)
3293 wprim(ixi^s,1:nw)=w(ixi^s,1:nw)
3294 call eos%to_primitive(ixi^
l,ixo^
l,wprim,x)
3306 integer,
intent(in) :: ixi^
l, ixo^
l
3307 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3308 double precision,
intent(out):: csound(ixi^s)
3310 double precision :: inv_rho, b2
3311 double precision :: prad_tensor(ixi^s, 1:
ndim, 1:
ndim)
3312 double precision :: prad_max(ixi^s)
3318 {
do ix^db=ixomin^db,ixomax^db \}
3319 inv_rho=1.d0/w(ix^
d,
rho_)
3320 prad_max(ix^
d) = maxval(prad_tensor(ix^
d,:,:))
3322 csound(ix^
d)=(eos%gamma*w(ix^
d,
p_)+b2+prad_max(ix^
d))*inv_rho
3325 {
do ix^db=ixomin^db,ixomax^db \}
3326 inv_rho=1.d0/w(ix^d,
rho_)
3327 prad_max(ix^d) = maxval(prad_tensor(ix^d,:,:))
3328 b2=(^
c&w(ix^d,
b^
c_)**2+)
3329 csound(ix^d)=(eos%gamma*w(ix^d,
p_)+b2+prad_max(ix^d))*inv_rho
3333 if(minval(csound(ixo^s))<smalldouble)
then
3334 print *,
'issue with squared speed and rad pressure'
3335 print *,minval(csound(ixo^s))
3336 print *,minval(prad_max(ixo^s))
3337 call mpistop(
"negative squared speed in get_csrad2 for dt")
3343 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3347 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3348 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3349 double precision,
intent(out):: csound(ixo^s)
3351 double precision :: adiabs(ixi^s), gammas(ixi^s)
3352 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3353 double precision :: cs2(ixi^s)
3373 call eos%get_csound2(w, x, ixi^
l, ixo^
l, cs2)
3378 {
do ix^db=ixomin^db,ixomax^db \}
3379 inv_rho=1.d0/w(ix^
d,
rho_)
3381 csound(ix^
d)=cs2(ix^
d)
3383 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3386 cfast2=b2*inv_rho+csound(ix^
d)
3387 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3389 if(avmincs2<zero) avmincs2=zero
3390 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3392 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3396 {
do ix^db=ixomin^db,ixomax^db \}
3397 inv_rho=1.d0/w(ix^d,
rho_)
3399 csound(ix^d)=cs2(ix^d)
3401 csound(ix^d)=gammas(ix^d)*adiabs(ix^d)*w(ix^d,
rho_)**(gammas(ix^d)-1.d0)
3403 b2=(^
c&w(ix^d,
b^
c_)**2+)
3404 cfast2=b2*inv_rho+csound(ix^d)
3405 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3406 if(avmincs2<zero) avmincs2=zero
3407 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3409 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3414 end subroutine mhd_get_csound_prim
3418 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3421 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3422 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3423 double precision,
intent(out):: csound(ixo^s)
3425 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3432 {
do ix^db=ixomin^db,ixomax^db \}
3437 cfast2=b2*inv_rho+csound(ix^
d)
3438 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3440 if(avmincs2<zero) avmincs2=zero
3441 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3443 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3447 {
do ix^db=ixomin^db,ixomax^db \}
3450 csound(ix^d)=eos%gamma*(w(ix^d,
p_)+block%equi_vars(ix^d,
equi_pe0_,b0i))*inv_rho
3451 b2=(^
c&w(ix^d,
b^
c_)**2+)
3452 cfast2=b2*inv_rho+csound(ix^d)
3453 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3454 if(avmincs2<zero) avmincs2=zero
3455 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3457 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3462 end subroutine mhd_get_csound_prim_split
3465 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3468 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3470 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3471 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3473 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3476 {
do ix^db=ixomin^db,ixomax^db\}
3477 inv_rho = 1.d0/w(ix^
d,
rho_)
3479 csound(ix^
d)=eos%gamma*w(ix^
d,
p_)*inv_rho
3480 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3481 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3482 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3483 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3486 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3487 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3488 if(avmincs2<zero) avmincs2=zero
3490 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3493 end subroutine mhd_get_csound_semirelati
3496 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3500 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3502 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3503 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3505 double precision :: adiabs(ixi^s), gammas(ixi^s)
3506 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3519 {
do ix^db=ixomin^db,ixomax^db\}
3520 inv_rho = 1.d0/w(ix^
d,
rho_)
3522 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3523 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3524 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3525 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3526 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3529 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3530 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3531 if(avmincs2<zero) avmincs2=zero
3533 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3536 end subroutine mhd_get_csound_semirelati_noe
3549 integer,
intent(in) :: ixi^
l, ixo^
l
3550 double precision,
intent(in) :: w(ixi^s, 1:nw)
3551 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3552 double precision,
intent(out):: prad(ixi^s, 1:
ndim, 1:
ndim)
3561 integer,
intent(in) :: ixi^
l, ixo^
l
3562 double precision,
intent(in) :: w(ixi^s, 1:nw)
3563 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3564 double precision,
intent(out) :: pth_plus_prad(ixi^s)
3566 double precision :: wprim(ixi^s, 1:nw)
3567 double precision :: prad_tensor(ixi^s, 1:
ndim, 1:
ndim)
3568 double precision :: prad_max(ixi^s)
3571 wprim(ixi^s,1:nw)=w(ixi^s,1:nw)
3572 call eos%to_primitive(ixi^
l,ixo^
l,wprim,x)
3574 {
do ix^
d = ixomin^
d,ixomax^
d\}
3575 prad_max(ix^
d) = maxval(prad_tensor(ix^
d,:,:))
3577 pth_plus_prad(ixo^s) = wprim(ixo^s,
p_) + prad_max(ixo^s)
3585 integer,
intent(in) :: ixi^
l, ixo^
l
3586 double precision,
intent(in) :: w(ixi^s, 1:nw)
3587 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3588 double precision,
intent(out):: trad(ixi^s)
3595 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3599 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3601 double precision,
intent(in) :: wc(ixi^s,nw)
3603 double precision,
intent(in) :: w(ixi^s,nw)
3604 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3605 double precision,
intent(out) :: f(ixi^s,nwflux)
3607 double precision :: vhall(ixi^s,1:
ndir)
3608 double precision :: ptotal
3609 double precision :: r(ixi^s), te(ixi^s), rho_loc(ixi^s)
3610 double precision :: bvec(ixi^s,1:
ndir)
3611 double precision :: bgradt(ixi^s), gradtperp_mag(ixi^s)
3612 double precision :: nperp(ixi^s,1:
ndir)
3613 logical :: use_perp_flux
3614 integer :: iw, ix^
d, idir
3617 {
do ix^db=ixomin^db,ixomax^db\}
3630 {
do ix^db=ixomin^db,ixomax^db\}
3634 ^
c&f(ix^d,
m^
c_)=wc(ix^d,
mom(idim))*w(ix^d,
m^
c_)-w(ix^d,mag(idim))*w(ix^d,
b^
c_)\
3635 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3637 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3640 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3641 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3643 ^
c&f(ix^d,
b^
c_)=w(ix^d,
mom(idim))*w(ix^d,
b^
c_)-w(ix^d,mag(idim))*w(ix^d,
m^
c_)\
3647 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3648 {
do ix^db=ixomin^db,ixomax^db\}
3649 if(total_energy)
then
3651 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3652 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3655 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*w(ix^d,
b^
c_)-vhall(ix^d,^
c)*w(ix^d,mag(idim))\
3660 {
do ix^db=ixomin^db,ixomax^db\}
3661 f(ix^d,mag(idim))=w(ix^d,
psi_)
3663 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3668 {
do ix^db=ixomin^db,ixomax^db\}
3669 f(ix^d,
r_e)=w(ix^d,
mom(idim))*wc(ix^d,
r_e)
3674 f(ixo^s,
fip_) = w(ixo^s,
mom(idim)) * wc(ixo^s,
fip_)
3678 {
do ix^db=ixomin^db,ixomax^db\}
3684 if(use_perp_flux)
then
3686 call eos%get_Rfactor(w,x,ixi^l,ixi^l,r)
3687 te(ixi^s)=w(ixi^s,
p_)/(r(ixi^s)*rho_loc(ixi^s))
3688 {
do ix^db=ixomin^db,ixomax^db\}
3690 bvec(ix^d,idir)=w(ix^d,mag(idir))
3693 call mhd_get_hyperbolic_tc_geometry(ixi^l,ixo^l,te,bvec,bgradt,gradtperp_mag,nperp)
3697 {
do ix^db=ixomin^db,ixomax^db\}
3698 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qpar_)*w(ix^d,mag(idim))/(dsqrt(^
c&w(ix^d,
b^
c_)**2+)+smalldouble)
3700 if(use_perp_flux)
then
3701 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qperp_)*nperp(ix^d,idim)
3706 end subroutine mhd_get_flux
3710 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3715 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3717 double precision,
intent(in) :: wc(ixi^s,nw)
3719 double precision,
intent(in) :: w(ixi^s,nw)
3720 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3721 double precision,
intent(out) :: f(ixi^s,nwflux)
3723 double precision :: vhall(ixi^s,1:
ndir)
3724 double precision :: adiabs(ixi^s), gammas(ixi^s)
3737 {
do ix^db=ixomin^db,ixomax^db\}
3743 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3748 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3749 {
do ix^db=ixomin^db,ixomax^db\}
3751 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*w(ix^d,
b^
c_)-vhall(ix^d,^
c)*w(ix^d,mag(idim))\
3755 {
do ix^db=ixomin^db,ixomax^db\}
3756 f(ix^d,mag(idim))=w(ix^d,
psi_)
3758 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3762 f(ixo^s,
fip_) = w(ixo^s,
mom(idim)) * wc(ixo^s,
fip_)
3766 {
do ix^db=ixomin^db,ixomax^db\}
3770 end subroutine mhd_get_flux_noe
3773 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3777 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3779 double precision,
intent(in) :: wc(ixi^s,nw)
3781 double precision,
intent(in) :: w(ixi^s,nw)
3782 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3783 double precision,
intent(out) :: f(ixi^s,nwflux)
3785 double precision :: vhall(ixi^s,1:
ndir)
3786 double precision :: r(ixi^s), te(ixi^s), rho_loc(ixi^s)
3787 double precision :: bvec(ixi^s,1:
ndir)
3788 double precision :: bgradt(ixi^s), gradtperp_mag(ixi^s)
3789 double precision :: nperp(ixi^s,1:
ndir)
3790 logical :: use_perp_flux
3791 integer :: iw, ix^
d, idir
3793 {
do ix^db=ixomin^db,ixomax^db\}
3806 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3807 {
do ix^db=ixomin^db,ixomax^db\}
3809 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*w(ix^d,
b^
c_)-vhall(ix^d,^
c)*w(ix^d,mag(idim))\
3813 {
do ix^db=ixomin^db,ixomax^db\}
3814 f(ix^d,mag(idim))=w(ix^d,
psi_)
3816 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3820 f(ixo^s,
fip_) = w(ixo^s,
mom(idim)) * wc(ixo^s,
fip_)
3824 {
do ix^db=ixomin^db,ixomax^db\}
3829 if(use_perp_flux)
then
3831 call eos%get_Rfactor(w,x,ixi^l,ixi^l,r)
3832 te(ixi^s)=w(ixi^s,
p_)/(r(ixi^s)*rho_loc(ixi^s))
3833 {
do ix^db=ixomin^db,ixomax^db\}
3835 bvec(ix^d,idir)=w(ix^d,mag(idir))
3838 call mhd_get_hyperbolic_tc_geometry(ixi^l,ixo^l,te,bvec,bgradt,gradtperp_mag,nperp)
3841 {
do ix^db=ixomin^db,ixomax^db\}
3842 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qpar_)*w(ix^d,mag(idim))/(dsqrt(^
c&w(ix^d,
b^
c_)**2+)+smalldouble)
3844 if(use_perp_flux)
then
3845 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qperp_)*nperp(ix^d,idim)
3850 end subroutine mhd_get_flux_hde
3857 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
3861 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3863 double precision,
intent(in) :: wc(ixi^s,nw)
3865 double precision,
intent(in) :: w(ixi^s,nw)
3866 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3867 double precision,
intent(out) :: f(ixi^s,nwflux)
3869 double precision :: vhall(ixi^s,1:
ndir)
3870 double precision :: ptotal, btotal(ixo^s,1:
ndir)
3871 double precision :: r(ixi^s), te(ixi^s), rho_loc(ixi^s)
3872 double precision :: bvec(ixi^s,1:
ndir)
3873 double precision :: bgradt(ixi^s), gradtperp_mag(ixi^s)
3874 double precision :: nperp(ixi^s,1:
ndir)
3875 logical :: use_perp_flux
3876 integer :: iw, ix^
d, idir
3878 {
do ix^db=ixomin^db,ixomax^db\}
3886 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3890 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
3894 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
3895 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3897 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
3901 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3904 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
3911 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
3912 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
3917 {
do ix^db=ixomin^db,ixomax^db\}
3918 f(ix^d,mag(idim))=w(ix^d,
psi_)
3920 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3925 {
do ix^db=ixomin^db,ixomax^db\}
3926 f(ix^d,
r_e)=w(ix^d,
mom(idim))*wc(ix^d,
r_e)
3931 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3932 {
do ix^db=ixomin^db,ixomax^db\}
3934 ^
c&f(ix^d,
b^
c_)=f(ix^d,
b^
c_)+vhall(ix^d,idim)*btotal(ix^d,^
c)-btotal(ix^d,idim)*vhall(ix^d,^
c)\
3935 if(total_energy)
then
3937 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
3938 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3943 f(ixo^s,
fip_) = w(ixo^s,
mom(idim)) * wc(ixo^s,
fip_)
3947 {
do ix^db=ixomin^db,ixomax^db\}
3952 if(use_perp_flux)
then
3954 call eos%get_Rfactor(w,x,ixi^l,ixi^l,r)
3955 te(ixi^s)=w(ixi^s,
p_)/(r(ixi^s)*rho_loc(ixi^s))
3956 {
do ix^db=ixomin^db,ixomax^db\}
3958 bvec(ix^d,idir)=btotal(ix^d,idir)
3961 call mhd_get_hyperbolic_tc_geometry(ixi^l,ixo^l,te,bvec,bgradt,gradtperp_mag,nperp)
3964 {
do ix^db=ixomin^db,ixomax^db\}
3965 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qpar_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
3967 if(use_perp_flux)
then
3968 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qperp_)*nperp(ix^d,idim)
3973 end subroutine mhd_get_flux_split
3976 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
3980 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3982 double precision,
intent(in) :: wc(ixi^s,nw)
3984 double precision,
intent(in) :: w(ixi^s,nw)
3985 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3986 double precision,
intent(out) :: f(ixi^s,nwflux)
3987 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
3988 double precision :: r(ixi^s), te(ixi^s), rho_loc(ixi^s)
3989 double precision :: bvec(ixi^s,1:
ndir)
3990 double precision :: bgradt(ixi^s), gradtperp_mag(ixi^s)
3991 double precision :: nperp(ixi^s,1:
ndir)
3992 logical :: use_perp_flux
3993 integer :: iw, ix^
d, idir
3995 {
do ix^db=ixomin^db,ixomax^db\}
4000 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4001 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4002 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4007 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4012 e2=(^
c&e(ix^
d,^
c)**2+)
4019 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4020 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4021 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4024 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4025 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4034 eos%gamma*w(ix^
d,
p_)*eos%inv_gamma_minus_1)+sa(ix^
d,idim)
4038 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4040 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+w(ix^
d,
p_)+half*((^
c&w(ix^
d,
b^
c_)**2+)+e2*inv_squared_c)
4047 {
do ix^db=ixomin^db,ixomax^db\}
4048 f(ix^d,mag(idim))=w(ix^d,
psi_)
4050 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4054 f(ixo^s,
fip_) = w(ixo^s,
mom(idim)) * wc(ixo^s,
fip_)
4058 {
do ix^db=ixomin^db,ixomax^db\}
4063 if(use_perp_flux)
then
4065 call eos%get_Rfactor(w,x,ixi^l,ixi^l,r)
4066 te(ixi^s)=w(ixi^s,
p_)/(r(ixi^s)*rho_loc(ixi^s))
4067 {
do ix^db=ixomin^db,ixomax^db\}
4069 bvec(ix^d,idir)=w(ix^d,mag(idir))
4072 call mhd_get_hyperbolic_tc_geometry(ixi^l,ixo^l,te,bvec,bgradt,gradtperp_mag,nperp)
4075 {
do ix^db=ixomin^db,ixomax^db\}
4076 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qpar_)*w(ix^d,mag(idim))/(dsqrt(^
c&w(ix^d,
b^
c_)**2+)+smalldouble)
4078 if(use_perp_flux)
then
4079 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
qperp_)*nperp(ix^d,idim)
4084 end subroutine mhd_get_flux_semirelati
4086 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4091 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4093 double precision,
intent(in) :: wc(ixi^s,nw)
4095 double precision,
intent(in) :: w(ixi^s,nw)
4096 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4097 double precision,
intent(out) :: f(ixi^s,nwflux)
4099 double precision :: adiabs(ixi^s), gammas(ixi^s)
4100 double precision :: e(ixo^s,1:
ndir),e2
4113 {
do ix^db=ixomin^db,ixomax^db\}
4118 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4119 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4120 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4121 e2=(^
c&e(ix^
d,^
c)**2+)
4126 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4136 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4138 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)+half*((^
c&w(ix^
d,
b^
c_)**2+)+e2*inv_squared_c)
4145 {
do ix^db=ixomin^db,ixomax^db\}
4146 f(ix^d,mag(idim))=w(ix^d,
psi_)
4148 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4152 f(ixo^s,
fip_) = w(ixo^s,
mom(idim)) * wc(ixo^s,
fip_)
4156 {
do ix^db=ixomin^db,ixomax^db\}
4160 end subroutine mhd_get_flux_semirelati_noe
4168 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x)
4170 integer,
intent(in) :: ixi^
l, ixo^
l
4171 double precision,
intent(in) :: qdt
4172 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4173 double precision,
intent(inout) :: w(ixi^s,1:nw)
4175 double precision :: tmp(ixi^s),btot2(ixi^s)
4176 double precision :: jxbxb(ixi^s,1:3)
4178 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4181 where (btot2(ixo^s)>smalldouble )
4182 tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=
ndim+1) / btot2(ixo^s)
4189 w(ixo^s,
e_)=w(ixo^s,
e_)- qdt*tmp(ixo^s)
4191 end subroutine add_source_ambipolar_internal_energy
4194 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4197 integer,
intent(in) :: ixi^
l, ixo^
l
4198 double precision,
intent(in) :: w(ixi^s,nw)
4199 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4200 double precision,
intent(out) :: res(ixi^s,1:3)
4202 double precision :: btot(ixi^s,1:3)
4203 double precision :: current(ixi^s,7-2*
ndir:3)
4204 double precision :: tmp(ixi^s),b2(ixi^s)
4205 integer :: idir, idirmin
4215 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4219 btot(ixo^s, idir) = w(ixo^s,mag(idir))
4223 tmp(ixo^s)= sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4224 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4226 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4229 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4234 where (b2(ixo^s)<smalldouble )
4235 res(ixo^s,idir) = zero
4238 end subroutine mhd_get_jxbxb
4244 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4248 integer,
intent(in) :: ixi^
l,ixo^
l,igrid,nflux
4249 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4250 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4251 double precision,
intent(in) :: my_dt
4252 logical,
intent(in) :: fix_conserve_at_step
4254 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4255 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4256 double precision :: fe(ixi^s,
sdim:3)
4257 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4258 integer :: i, ixa^
l, ie_
4265 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4282 btot(ixa^s,1:3) = 0.d0
4284 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4288 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4289 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4291 wres(ixo^s,
e_)=-tmp2(ixo^s)
4298 ff(ixa^s,1) = tmp(ixa^s,2)
4299 ff(ixa^s,2) = -tmp(ixa^s,1)
4301 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4302 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4303 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4306 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4308 ixamin^
d=ixomin^
d-1;
4309 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4319 ff(ixa^s,2) = tmp(ixa^s,3)
4320 ff(ixa^s,3) = -tmp(ixa^s,2)
4321 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4322 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4324 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4327 ff(ixa^s,1) = -tmp(ixa^s,3)
4329 ff(ixa^s,3) = tmp(ixa^s,1)
4330 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4331 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4332 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4338 ff(ixa^s,2) = tmp(ixa^s,3)
4339 ff(ixa^s,3) = -tmp(ixa^s,2)
4340 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4341 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4343 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4345 ff(ixa^s,1) = -tmp(ixa^s,3)
4347 ff(ixa^s,3) = tmp(ixa^s,1)
4348 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4349 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4350 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4355 ff(ixa^s,1) = tmp(ixa^s,2)
4356 ff(ixa^s,2) = -tmp(ixa^s,1)
4358 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4359 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4360 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4365 if(fix_conserve_at_step)
then
4366 fluxall=my_dt*fluxall
4373 end subroutine sts_set_source_ambipolar
4376 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4379 integer,
intent(in) :: ixi^
l, ixo^
l
4380 double precision,
intent(in) :: w(ixi^s,1:nw)
4381 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4383 double precision,
intent(in) :: ecc(ixi^s,1:3)
4384 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4385 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4387 integer :: hxc^
l,ixc^
l,ixa^
l
4388 integer :: idim1,idim2,idir,ix^
d
4394 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4396 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4397 ixamin^
d=ixcmin^
d+ix^
d;
4398 ixamax^
d=ixcmax^
d+ix^
d;
4399 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4401 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4407 ixcmin^d=ixomin^d-1;
4414 hxc^l=ixc^l-kr(idim2,^d);
4416 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4417 +lvc(idim1,idim2,idir)&
4422 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4425 end subroutine update_faces_ambipolar
4431 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4434 integer,
intent(in) :: ixi^
l, ixo^
l
4435 double precision,
dimension(ixI^S,1:3),
intent(inout) :: ff
4436 double precision,
intent(out) :: src(ixi^s)
4438 double precision :: ffc(ixi^s,1:
ndim)
4439 double precision :: dxinv(
ndim)
4440 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4448 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4450 ixbmin^
d=ixcmin^
d+ix^
d;
4451 ixbmax^
d=ixcmax^
d+ix^
d;
4454 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4456 call mpistop(
"to generalize using volume averaging")
4459 ff(ixi^s,1:ndim)=0.d0
4461 ixb^l=ixo^l-kr(idims,^d);
4462 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4464 if({ ix^d==0 .and. ^d==idims | .or.})
then
4465 ixbmin^d=ixcmin^d-ix^d;
4466 ixbmax^d=ixcmax^d-ix^d;
4467 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4470 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4473 if(slab_uniform)
then
4475 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4476 ixb^l=ixo^l-kr(idims,^d);
4477 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4481 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4482 ixb^l=ixo^l-kr(idims,^d);
4483 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4485 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4487 end subroutine get_flux_on_cell_face
4491 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4494 integer,
intent(in) :: ixi^
l, ixo^
l
4495 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4496 double precision,
intent(in) :: w(ixi^s,1:nw)
4497 double precision :: dtnew
4499 double precision :: coef
4500 double precision :: dxarr(
ndim)
4501 double precision :: tmp(ixi^s)
4507 coef = maxval(dabs(tmp(ixo^s)))
4514 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4516 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4519 end function get_ambipolar_dt
4527 integer,
intent(in) :: ixi^
l, ixo^
l
4528 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4529 double precision,
intent(inout) :: res(ixi^s)
4530 double precision :: tmp(ixi^s)
4531 double precision :: rho(ixi^s)
4538 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4543 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4551 integer,
intent(in) :: ixi^
l, ixo^
l
4552 double precision,
intent(in) :: qdt,dtfactor
4553 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4554 double precision,
intent(inout) :: w(ixi^s,1:nw)
4555 logical,
intent(in) :: qsourcesplit
4556 logical,
intent(inout) :: active
4563 if (.not. qsourcesplit)
then
4567 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4571 call add_equi_terms(qdt,dtfactor,ixi^
l,ixo^
l,wct,w,x,wctprim)
4577 call add_hyperbolic_tc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4586 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wct,w,x,wctprim)
4590 if (abs(
mhd_eta)>smalldouble)
then
4592 call add_source_res_exp(qdt,ixi^
l,ixo^
l,wct,w,x)
4597 call add_source_ambi_exp(qdt,ixi^
l,ixo^
l,wct,w,x)
4602 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4608 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4612 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4619 select case (type_divb)
4624 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4627 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4630 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4631 case (divb_janhunen)
4633 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4634 case (divb_lindejanhunen)
4636 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4637 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4638 case (divb_lindepowel)
4640 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4641 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4642 case (divb_lindeglm)
4644 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4645 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4646 case (divb_multigrid)
4651 call mpistop(
'Unknown divB fix')
4658 w,x,qsourcesplit,active,
rc_fl)
4668 w,x,gravity_energy,qsourcesplit,active)
4677 call mhd_add_radiation_source(qdt,ixi^
l,ixo^
l,wct,wctprim,w,x,qsourcesplit,active)
4681 if(eos%eos_type ==
'PI')
then
4682 if(.not.qsourcesplit)
then
4684 call eos%update_eos(ixi^
l,ixo^
l,w,x)
4688 end subroutine mhd_add_source
4690 subroutine mhd_add_radiation_source(qdt,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4696 integer,
intent(in) :: ixi^
l, ixo^
l
4697 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
4698 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw)
4699 double precision,
intent(inout) :: w(ixi^s,1:nw)
4700 logical,
intent(in) :: qsourcesplit
4701 logical,
intent(inout) :: active
4707 end subroutine mhd_add_radiation_source
4710 subroutine add_equi_terms(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x,wCTprim)
4715 integer,
intent(in) :: ixi^
l, ixo^
l
4716 double precision,
intent(in) :: qdt,dtfactor
4717 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4718 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
4719 double precision,
intent(inout) :: w(ixi^s,1:nw)
4721 double precision :: divv(ixi^s)
4722 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4723 double precision :: gravity_field(ixi^s,1:
ndim)
4735 divv(ixo^s)=divv(ixo^s)*eos%gamma*eos%inv_gamma_minus_1
4746 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*
block%equi_vars(ixo^s,
equi_rho0_,0)*gravity_field(ixo^s,idir)*eos%inv_gamma_minus_1
4755 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4760 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*axb(ixo^s,idir)*eos%inv_gamma_minus_1
4766 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*
block%equi_vars(ixo^s,
equi_rho0_,0)*gravity_field(ixo^s,idir)*eos%inv_gamma_minus_1
4775 w(ixo^s,
e_)=w(ixo^s,
e_)-qdt*wctprim(ixo^s,
mom(idir))*
block%equi_vars(ixo^s,
equi_rho0_,0)*gravity_field(ixo^s,idir)*eos%inv_gamma_minus_1
4779 end subroutine add_equi_terms
4781 subroutine mhd_get_hyperbolic_tc_geometry(ixI^L,ixO^L,Te,Bvec,bgradT,gradTperp_mag,nperp)
4784 integer,
intent(in) :: ixi^
l,ixo^
l
4785 double precision,
intent(in) :: te(ixi^s)
4786 double precision,
intent(in) :: bvec(ixi^s,1:
ndir)
4787 double precision,
intent(out) :: bgradt(ixi^s), gradtperp_mag(ixi^s)
4788 double precision,
intent(out) :: nperp(ixi^s,1:
ndir)
4790 double precision :: bmag, bunitvec(
ndir), gradt(
ndir), gradt_perp(
ndir)
4791 double precision :: gradt_cell(ixi^s,1:
ndir)
4792 integer :: ix^
d, idir
4797 call gradient(te,ixi^
l,ixo^
l,idir,gradt_cell(ixi^s,idir))
4802 do ix2=ixomin2,ixomax2
4803 do ix1=ixomin1,ixomax1
4806 bmag=bmag+bvec(ix^
d,idir)**2
4810 if(bmag>smalldouble)
then
4812 bunitvec(idir)=bvec(ix^
d,idir)/bmag
4820 gradt(1)=((8.d0*(te(ix1+1,ix2)-te(ix1-1,ix2))-te(ix1+2,ix2)+te(ix1-2,ix2))/12.d0)/
block%ds(ix^
d,1)
4821 gradt(2)=((8.d0*(te(ix1,ix2+1)-te(ix1,ix2-1))-te(ix1,ix2+2)+te(ix1,ix2-2))/12.d0)/
block%ds(ix^
d,2)
4822 if(
ndir>2) gradt(3)=zero
4825 gradt(idir)=gradt_cell(ix^
d,idir)
4831 bgradt(ix^
d)=bgradt(ix^
d)+bunitvec(idir)*gradt(idir)
4835 gradt_perp(idir)=gradt(idir)-bgradt(ix^
d)*bunitvec(idir)
4838 gradtperp_mag(ix^
d)=zero
4840 gradtperp_mag(ix^
d)=gradtperp_mag(ix^
d)+gradt_perp(idir)**2
4842 gradtperp_mag(ix^
d)=dsqrt(gradtperp_mag(ix^
d))
4844 if(gradtperp_mag(ix^
d)>smalldouble)
then
4846 nperp(ix^
d,idir)=gradt_perp(idir)/gradtperp_mag(ix^
d)
4849 gradtperp_mag(ix^
d)=zero
4851 nperp(ix^
d,idir)=zero
4858 do ix3=ixomin3,ixomax3
4859 do ix2=ixomin2,ixomax2
4860 do ix1=ixomin1,ixomax1
4861 bmag=dsqrt(bvec(ix^
d,1)**2+bvec(ix^
d,2)**2+bvec(ix^
d,3)**2)
4862 if(bmag>smalldouble)
then
4863 bunitvec(1)=bvec(ix^
d,1)/bmag
4864 bunitvec(2)=bvec(ix^
d,2)/bmag
4865 bunitvec(3)=bvec(ix^
d,3)/bmag
4873 gradt(1)=((8.d0*(te(ix1+1,ix2,ix3)-te(ix1-1,ix2,ix3))-te(ix1+2,ix2,ix3)+te(ix1-2,ix2,ix3))/12.d0)/
block%ds(ix^
d,1)
4874 gradt(2)=((8.d0*(te(ix1,ix2+1,ix3)-te(ix1,ix2-1,ix3))-te(ix1,ix2+2,ix3)+te(ix1,ix2-2,ix3))/12.d0)/
block%ds(ix^
d,2)
4875 gradt(3)=((8.d0*(te(ix1,ix2,ix3+1)-te(ix1,ix2,ix3-1))-te(ix1,ix2,ix3+2)+te(ix1,ix2,ix3-2))/12.d0)/
block%ds(ix^
d,3)
4878 gradt(idir)=gradt_cell(ix^
d,idir)
4884 bgradt(ix^
d)=bgradt(ix^
d)+bunitvec(idir)*gradt(idir)
4888 gradt_perp(idir)=gradt(idir)-bgradt(ix^
d)*bunitvec(idir)
4891 gradtperp_mag(ix^
d)=dsqrt(gradt_perp(1)**2+gradt_perp(2)**2+gradt_perp(3)**2)
4892 if(gradtperp_mag(ix^
d)>smalldouble)
then
4894 nperp(ix^
d,idir)=gradt_perp(idir)/gradtperp_mag(ix^
d)
4897 gradtperp_mag(ix^
d)=zero
4899 nperp(ix^
d,idir)=zero
4906 end subroutine mhd_get_hyperbolic_tc_geometry
4908 subroutine add_hyperbolic_tc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4911 integer,
intent(in) :: ixi^
l,ixo^
l
4912 double precision,
intent(in) :: qdt
4913 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4914 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4915 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4917 double precision,
dimension(ixI^S) :: r,te,rho_loc,pth_loc
4918 double precision,
dimension(ixI^S,1:ndir) :: bvec
4919 double precision,
dimension(ixI^S) :: bgradt, gradtperp_mag
4920 double precision,
dimension(ixI^S,1:ndir) :: nperp
4921 double precision,
dimension(ixI^S) :: gradt_geom
4922 double precision,
parameter :: lnlambda_perp = 20.d0
4923 double precision,
parameter :: xe_prefac_cgs = 4.753567596681522d6
4924 double precision :: kappa_t5,kappa_t5_perp,kappa_t5_perp_eff
4925 double precision :: kappa_t7,f_sat,kappat5_bgradt,kappat5_gradtperp,tau,b2,fb,gradt1
4926 double precision :: bmag_loc,tloc,tcond,nloc_code,cchi,chi
4927 double precision :: cmax(
ndim),c2,cfast2,avmincs2(
ndim),inv_rho
4928 logical :: use_perp_source
4929 integer :: ix^
d,idir
4931 cchi = 0.823d0*(xe_prefac_cgs/lnlambda_perp) * &
4933 call eos%get_Rfactor(wct,x,ixi^
l,ixi^
l,r)
4934 {
do ix^db=iximin^db,iximax^db\}
4939 rho_loc(ix^
d)=wctprim(ix^
d,
rho_)
4940 pth_loc(ix^
d)=wctprim(ix^
d,
p_)
4942 te(ix^
d)=pth_loc(ix^
d)/(r(ix^
d)*rho_loc(ix^
d))
4946 {
do ix^db=ixomin^db,ixomax^db\}
4948 bvec(ix^d,idir)=wct(ix^d,mag(idir))+block%B0(ix^d,idir,0)
4952 {
do ix^db=ixomin^db,ixomax^db\}
4954 bvec(ix^d,idir)=wct(ix^d,mag(idir))
4959 call mhd_get_hyperbolic_tc_geometry(ixi^l,ixo^l,te,bvec,bgradt,gradtperp_mag,nperp)
4963 if(.not.slab_uniform)
then
4964 call gradient(te,ixi^l,ixo^l,1,gradt_geom)
4966 do ix1=ixomin1,ixomax1
4969 kappa_t7=kappa_t5*te(ix1)
4973 tcond = max(tcond, block%wextra(ix1,
tcoff_))
4976 kappa_t7=kappa_t5*tcond
4978 if(slab_uniform)
then
4979 gradt1=((8.d0*(te(ix1+1)-te(ix1-1))-te(ix1+2)+te(ix1-2))/12.d0)/block%ds(ix1,1)
4981 gradt1=gradt_geom(ix1)
4985 b2=b2+bvec(ix1,idir)**2
4987 if(b2>smalldouble**2)
then
4988 bgradt(ix1)=bvec(ix1,1)*gradt1/dsqrt(b2)
4992 kappat5_bgradt=kappa_t5*bgradt(ix1)
4993 inv_rho=1.d0/rho_loc(ix1)
4994 c2=eos%gamma*pth_loc(ix1)*inv_rho
4995 cfast2 = b2*inv_rho + c2
4996 avmincs2(1) = cfast2**2 - 4.0d0*c2*bvec(ix1,1)**2*inv_rho
4997 cmax(1) = sqrt(half*(cfast2 + sqrt(dabs(avmincs2(1)))))
4999 f_sat=one/(one+dabs(kappat5_bgradt)/(1.5d0*rho_loc(ix^d)*(pth_loc(ix^d)/rho_loc(ix^d))**1.5d0))
5000 tau=max(4.d0*dt, f_sat*kappa_t7*courantpar**2/(pth_loc(ix^d)*eos%inv_gamma_minus_1*cmax(1)**2))
5001 w(ix^d,
qpar_)=w(ix^d,
qpar_)-qdt*(f_sat*kappat5_bgradt+wct(ix^d,
qpar_))/tau
5004 max(4.d0*dt, kappa_t7*courantpar**2/(pth_loc(ix^d)*eos%inv_gamma_minus_1*cmax(1)**2))
5009 do ix2=ixomin2,ixomax2
5010 do ix1=ixomin1,ixomax1
5013 kappa_t7=kappa_t5*te(ix^d)
5017 tcond=max(tcond, block%wextra(ix^d,
tcoff_))
5020 kappa_t7 = kappa_t5*tcond
5022 kappat5_bgradt=kappa_t5*bgradt(ix^d)
5025 b2 = b2 + bvec(ix^d,idir)**2
5027 if(use_perp_source)
then
5037 kappa_t5_perp_eff=(one-fb)*kappa_t5
5038 kappa_t5_perp=kappa_t5_perp_eff
5040 bmag_loc = dsqrt(b2)
5041 tloc = max(te(ix^d), smalldouble)
5042 nloc_code = max(rho_loc(ix^d), smalldouble)
5043 chi = cchi*bmag_loc*tloc**1.5d0/nloc_code
5044 kappa_t5_perp_eff = kappa_t5/(one+chi**2)
5045 kappa_t5_perp = kappa_t5_perp_eff
5049 kappat5_gradtperp=kappa_t5_perp*gradtperp_mag(ix^d)
5051 inv_rho=1.d0/rho_loc(ix^d)
5052 c2=eos%gamma*pth_loc(ix^d)*inv_rho
5053 cfast2 = b2*inv_rho + c2
5055 avmincs2(idir)=cfast2**2-4.0d0*c2*bvec(ix^d,idir)**2*inv_rho
5056 cmax(idir)=sqrt(half*(cfast2+sqrt(dabs(avmincs2(idir)))))\
5059 f_sat=one/(one+dabs(kappat5_bgradt)/(1.5d0*rho_loc(ix^d)*(pth_loc(ix^d)/rho_loc(ix^d))**1.5d0))
5060 tau=max(4.d0*dt, f_sat*kappa_t7*courantpar**2/(pth_loc(ix^d)*eos%inv_gamma_minus_1*maxval(cmax(:))**2))
5061 w(ix^d,
qpar_)=w(ix^d,
qpar_)-qdt*(f_sat*kappat5_bgradt+wct(ix^d,
qpar_))/tau
5062 if(use_perp_source)
then
5066 tau=max(4.d0*dt, kappa_t7*courantpar**2/(pth_loc(ix^d)*eos%inv_gamma_minus_1*maxval(cmax(:))**2))
5068 if(use_perp_source)
then
5076 do ix3=ixomin3,ixomax3
5077 do ix2=ixomin2,ixomax2
5078 do ix1=ixomin1,ixomax1
5081 kappa_t7=kappa_t5*te(ix^d)
5085 tcond=max(tcond, block%wextra(ix^d,
tcoff_))
5088 kappa_t7 = kappa_t5*tcond
5090 kappat5_bgradt=kappa_t5*bgradt(ix^d)
5093 b2 = b2 + bvec(ix^d,idir)**2
5095 if(use_perp_source)
then
5105 kappa_t5_perp_eff=(one-fb)*kappa_t5
5106 kappa_t5_perp=kappa_t5_perp_eff
5108 bmag_loc = dsqrt(b2)
5109 tloc = max(te(ix^d), smalldouble)
5110 nloc_code = max(rho_loc(ix^d), smalldouble)
5111 chi = cchi*bmag_loc*tloc**1.5d0/nloc_code
5112 kappa_t5_perp_eff = kappa_t5/(one+chi**2)
5113 kappa_t5_perp = kappa_t5_perp_eff
5117 kappat5_gradtperp=kappa_t5_perp*gradtperp_mag(ix^d)
5119 inv_rho=1.d0/rho_loc(ix^d)
5120 c2=eos%gamma*pth_loc(ix^d)*inv_rho
5121 cfast2 = b2*inv_rho + c2
5123 avmincs2(idir)=cfast2**2-4.0d0*c2*bvec(ix^d,idir)**2*inv_rho
5124 cmax(idir)=sqrt(half*(cfast2+sqrt(dabs(avmincs2(idir)))))\
5127 f_sat=one/(one+dabs(kappat5_bgradt)/(1.5d0*rho_loc(ix^d)*(pth_loc(ix^d)/rho_loc(ix^d))**1.5d0))
5128 tau=max(4.d0*dt, f_sat*kappa_t7*courantpar**2/(pth_loc(ix^d)*eos%inv_gamma_minus_1*maxval(cmax(:))**2))
5129 w(ix^d,
qpar_)=w(ix^d,
qpar_)-qdt*(f_sat*kappat5_bgradt+wct(ix^d,
qpar_))/tau
5130 if(use_perp_source)
then
5134 tau=max(4.d0*dt, kappa_t7*courantpar**2/(pth_loc(ix^d)*eos%inv_gamma_minus_1*maxval(cmax(:))**2))
5136 if(use_perp_source)
then
5144 end subroutine add_hyperbolic_tc_source
5148 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
5150 integer,
intent(in) :: ixi^
l, ixo^
l
5151 double precision,
intent(in) :: w(ixi^s,1:nw)
5152 double precision,
intent(inout) :: jxb(ixi^s,3)
5153 double precision :: a(ixi^s,3),
b(ixi^s,3)
5155 double precision :: current(ixi^s,7-2*
ndir:3)
5156 integer :: idir, idirmin
5161 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
5165 b(ixo^s, idir) = w(ixo^s,mag(idir))
5174 a(ixo^s,idir)=current(ixo^s,idir)
5178 end subroutine get_lorentz_force
5182 integer,
intent(in) :: ixi^
l, ixo^
l
5183 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
5184 double precision,
intent(out) :: rho(ixi^s)
5189 rho(ixo^s) = w(ixo^s,
rho_)
5195 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
5198 integer,
intent(in) :: ixi^
l,ixo^
l, ie
5199 double precision,
intent(inout) :: w(ixi^s,1:nw)
5200 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5201 character(len=*),
intent(in) :: subname
5203 double precision :: rho(ixi^s)
5205 logical :: flag(ixi^s,1:nw)
5210 flag(ixo^s,ie)=.true.
5212 where(w(ixo^s,ie)<
small_e) flag(ixo^s,ie)=.true.
5214 if(any(flag(ixo^s,ie)))
then
5218 where(flag(ixo^s,ie)) w(ixo^s,ie)=
small_e - &
5221 where(flag(ixo^s,ie)) w(ixo^s,ie)=
small_e
5227 w(ixo^s,
e_)=w(ixo^s,
e_)*eos%gamma_minus_1
5230 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
5236 end subroutine mhd_handle_small_ei
5240 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x,wCTprim)
5243 integer,
intent(in) :: ixi^
l, ixo^
l
5244 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5245 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5246 double precision,
intent(inout) :: w(ixi^s,1:nw)
5248 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
5260 a(ixo^s,idir)=
block%J0(ixo^s,idir)
5265 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5268 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5274 if(total_energy)
then
5277 b(ixo^s,:)=wctprim(ixo^s,mag(:))
5286 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5289 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5294 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5298 call mhd_getv_hall(wct,x,ixi^
l,ixo^
l,a,.true.)
5303 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5306 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5310 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5318 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5323 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5329 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5331 end subroutine add_source_b0split
5334 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5338 integer,
intent(in) :: ixi^
l, ixo^
l
5339 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5340 double precision,
intent(inout) :: w(ixi^s,1:nw)
5341 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5343 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5344 integer :: idir, idirmin, ix^
d
5348 {
do ix^db=iximin^db,iximax^db\}
5350 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5351 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5352 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5354 call divvector(e,ixi^l,ixo^l,dive)
5356 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5359 {
do ix^db=ixomin^db,ixomax^db\}
5360 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5361 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5362 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5363 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5364 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5365 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5369 end subroutine add_source_semirelativistic
5372 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5376 integer,
intent(in) :: ixi^
l, ixo^
l
5377 double precision,
intent(in) :: qdt
5378 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5379 double precision,
intent(inout) :: w(ixi^s,1:nw)
5380 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5382 double precision :: divv(ixi^s), tmp
5394 {
do ix^db=ixomin^db,ixomax^db\}
5396 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5402 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x)
5405 if(fix_small_values)
then
5406 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5408 end subroutine add_source_internal_e
5411 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5416 integer,
intent(in) :: ixi^
l, ixo^
l
5417 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5418 double precision,
intent(inout) :: w(ixi^s,1:nw)
5419 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5421 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5422 double precision :: current(ixi^s,7-2*
ndir:3)
5423 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5424 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5425 integer :: idir, idirmin, idims, ix^
d
5430 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5442 j(ixo^s,idir)=current(ixo^s,idir)
5521 call add_source_ambipolar_internal_energy(qdt,ixi^
l,ixo^
l,wct,w,x)
5524 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_hydrodynamic_e')
5526 end subroutine add_source_hydrodynamic_e
5532 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5537 integer,
intent(in) :: ixi^
l, ixo^
l
5538 double precision,
intent(in) :: qdt
5539 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5540 double precision,
intent(inout) :: w(ixi^s,1:nw)
5542 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim
5543 double precision :: tmp(ixi^s),tmp2(ixi^s)
5546 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5547 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5548 double precision :: lapl_vec(ixi^s,1:
ndir)
5554 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5555 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5562 gradeta(ixo^s,1:
ndim)=zero
5567 gradeta(ixo^s,idim)=tmp(ixo^s)
5574 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5581 tmp(ixo^s)=lapl_vec(ixo^s,idir)*eta(ixo^s)
5585 do jdir=1,
ndim;
do kdir=idirmin,3
5586 if (
lvc(idir,jdir,kdir)/=0)
then
5587 if (
lvc(idir,jdir,kdir)==1)
then
5588 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5590 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5597 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5598 if(total_energy)
then
5599 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5605 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5608 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5610 end subroutine add_source_res1
5614 subroutine add_source_res_exp(qdt,ixI^L,ixO^L,wCT,w,x)
5619 integer,
intent(in) :: ixi^
l, ixo^
l
5620 double precision,
intent(in) :: qdt
5621 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5622 double precision,
intent(inout) :: w(ixi^s,1:nw)
5625 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5626 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5627 integer :: ixa^
l,idir,idirmin,idirmin1
5631 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5632 call mpistop(
"Error in add_source_res_exp: Non-conforming input limits")
5642 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5647 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5656 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5659 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5664 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5666 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5668 if(total_energy)
then
5671 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5672 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5675 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5679 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res_exp')
5680 end subroutine add_source_res_exp
5685 subroutine add_source_ambi_exp(qdt,ixI^L,ixO^L,wCT,w,x)
5690 integer,
intent(in) :: ixi^
l, ixo^
l
5691 double precision,
intent(in) :: qdt
5692 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5693 double precision,
intent(inout) :: w(ixi^s,1:nw)
5695 double precision :: current(ixi^s,1:3),curlj(ixi^s,1:3)
5696 double precision :: tmpvec(ixi^s,1:3),tmp(ixi^s),btot2(ixi^s)
5697 integer :: ixa^
l,idir,idirmin1
5701 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5702 call mpistop(
"Error in add_source_ambi_exp: Non-conforming input limits")
5706 call mhd_get_jxbxb(wct,x,ixi^
l,ixa^
l,current)
5720 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5723 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5730 where (btot2(ixa^s)>smalldouble )
5731 tmp(ixa^s) = sum(current(ixa^s,1:3)**2,dim=
ndim+1) / btot2(ixa^s)
5738 tmp(ixo^s)=-qdt*tmp(ixo^s)
5739 if(total_energy)
then
5742 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5743 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5746 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5750 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_ambi_exp')
5751 end subroutine add_source_ambi_exp
5755 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5759 integer,
intent(in) :: ixi^
l, ixo^
l
5760 double precision,
intent(in) :: qdt
5761 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5762 double precision,
intent(inout) :: w(ixi^s,1:nw)
5764 double precision :: current(ixi^s,7-2*
ndir:3)
5765 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5766 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5769 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5770 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5773 tmpvec(ixa^s,1:
ndir)=zero
5775 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5779 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5782 tmpvec(ixa^s,1:
ndir)=zero
5783 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5787 tmpvec2(ixa^s,1:
ndir)=zero
5788 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5791 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5794 if(total_energy)
then
5797 tmpvec2(ixa^s,1:
ndir)=zero
5798 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5799 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5800 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5801 end do;
end do;
end do
5803 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5804 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5807 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5809 end subroutine add_source_hyperres
5811 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5818 integer,
intent(in) :: ixi^
l, ixo^
l
5819 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5820 double precision,
intent(inout) :: w(ixi^s,1:nw)
5822 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5843 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5846 if(total_energy)
then
5855 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5864 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5868 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5870 end subroutine add_source_glm
5873 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5876 integer,
intent(in) :: ixi^
l, ixo^
l
5877 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5878 double precision,
intent(inout) :: w(ixi^s,1:nw)
5880 double precision :: divb(ixi^s), ba(1:
ndir)
5881 integer :: idir, ix^
d
5887 {
do ix^db=ixomin^db,ixomax^db\}
5892 if (total_energy)
then
5898 {
do ix^db=ixomin^db,ixomax^db\}
5900 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5902 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5903 if (total_energy)
then
5905 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5910 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5912 end subroutine add_source_powel
5914 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5919 integer,
intent(in) :: ixi^
l, ixo^
l
5920 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5921 double precision,
intent(inout) :: w(ixi^s,1:nw)
5923 double precision :: divb(ixi^s)
5924 integer :: idir, ix^
d
5929 {
do ix^db=ixomin^db,ixomax^db\}
5934 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5936 end subroutine add_source_janhunen
5938 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5943 integer,
intent(in) :: ixi^
l, ixo^
l
5944 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5945 double precision,
intent(inout) :: w(ixi^s,1:nw)
5947 double precision :: divb(ixi^s),graddivb(ixi^s)
5948 integer :: idim, idir, ixp^
l, i^
d, iside
5949 logical,
dimension(-1:1^D&) :: leveljump
5957 if(i^
d==0|.and.) cycle
5958 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5959 leveljump(i^
d)=.true.
5961 leveljump(i^
d)=.false.
5970 i^dd=kr(^dd,^d)*(2*iside-3);
5971 if (leveljump(i^dd))
then
5973 ixpmin^d=ixomin^d-i^d
5975 ixpmax^d=ixomax^d-i^d
5986 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5988 {
do i^db=ixpmin^db,ixpmax^db\}
5990 graddivb(i^d)=graddivb(i^d)*
divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5992 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5994 if (typedivbdiff==
'all' .and. total_energy)
then
5996 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
6001 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
6003 end subroutine add_source_linde
6010 integer,
intent(in) :: ixi^
l, ixo^
l
6011 double precision,
intent(in) :: w(ixi^s,1:nw)
6012 double precision :: divb(ixi^s), dsurface(ixi^s)
6014 double precision :: invb(ixo^s)
6015 integer :: ixa^
l,idims
6017 call get_divb(w,ixi^
l,ixo^
l,divb)
6019 where(invb(ixo^s)/=0.d0)
6020 invb(ixo^s)=1.d0/invb(ixo^s)
6023 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
6025 ixamin^
d=ixomin^
d-1;
6026 ixamax^
d=ixomax^
d-1;
6027 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
6029 ixa^
l=ixo^
l-
kr(idims,^
d);
6030 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
6032 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
6033 block%dvolume(ixo^s)/dsurface(ixo^s)
6044 integer,
intent(in) :: ixo^
l, ixi^
l
6045 double precision,
intent(in) :: w(ixi^s,1:nw)
6046 integer,
intent(out) :: idirmin
6049 double precision :: current(ixi^s,7-2*
ndir:3)
6050 integer :: idir, idirmin0
6056 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
6057 block%J0(ixo^s,idirmin0:3)
6061 subroutine mhd_get_dt(wprim,ixI^L,ixO^L,dtnew,dx^D,x)
6069 integer,
intent(in) :: ixi^
l, ixo^
l
6070 double precision,
intent(inout) :: dtnew
6071 double precision,
intent(in) ::
dx^
d
6072 double precision,
intent(in) :: wprim(ixi^s,1:nw)
6073 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6075 double precision :: dxarr(
ndim)
6076 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
6077 integer :: idirmin,idim
6095 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
6098 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
6120 dtnew=min(
dtdiffpar*get_ambipolar_dt(wprim,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
6131 end subroutine mhd_get_dt
6136 subroutine mhd_fld_implicit_update(dtfactor,qdt,qtC,psa,psb)
6141 double precision,
intent(in) :: qdt
6142 double precision,
intent(in) :: qtc
6143 double precision,
intent(in) :: dtfactor
6146 end subroutine mhd_fld_implicit_update
6148 subroutine mhd_fld_evaluate_implicit(qtC,psa)
6152 double precision,
intent(in) :: qtc
6155 end subroutine mhd_fld_evaluate_implicit
6162 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6168 integer,
intent(in) :: ixi^
l, ixo^
l
6169 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6170 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6172 double precision :: adiabs(ixi^s), gammas(ixi^s)
6173 double precision :: tmp,tmp1,invr,cot
6175 integer :: mr_,mphi_
6176 integer :: br_,bphi_
6179 br_=mag(1); bphi_=mag(1)-1+
phi_
6196 {
do ix^db=ixomin^db,ixomax^db\}
6199 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6204 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6206 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6209 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6210 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6211 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6212 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6213 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6215 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6216 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6217 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6220 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6225 {
do ix^db=ixomin^db,ixomax^db\}
6227 if(local_timestep)
then
6228 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6233 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6235 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6239 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6242 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6243 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6247 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6253 cot=1.d0/tan(x(ix^d,2))
6257 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6258 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6260 if(.not.stagger_grid)
then
6261 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6263 tmp=tmp+wprim(ix^d,
psi_)*cot
6265 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6270 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6271 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6272 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6274 if(.not.stagger_grid)
then
6275 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6277 tmp=tmp+wprim(ix^d,
psi_)*cot
6279 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6282 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6283 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6284 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6285 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6286 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6288 if(.not.stagger_grid)
then
6289 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6290 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6291 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6292 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6293 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6300 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6303 end subroutine mhd_add_source_geom
6310 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6316 integer,
intent(in) :: ixi^
l, ixo^
l
6317 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6318 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6320 double precision :: adiabs(ixi^s), gammas(ixi^s)
6321 double precision :: tmp,tmp1,tmp2,invr,cot,ef(ixo^s,1:
ndir)
6323 integer :: mr_,mphi_
6324 integer :: br_,bphi_
6327 br_=mag(1); bphi_=mag(1)-1+
phi_
6344 {
do ix^db=ixomin^db,ixomax^db\}
6347 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6354 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)
6358 ef(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
6359 ef(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
6360 ef(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6365 ef(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6371 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
6372 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&ef(ix^
d,^
c)**2+)*inv_squared_c) -&
6373 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
6374 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6375 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6376 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+ef(ix^
d,
phi_)*ef(ix^
d,1)*inv_squared_c)
6378 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6379 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6380 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6383 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
6384 (^
c&ef(ix^
d,^
c)**2+)*inv_squared_c))
6389 {
do ix^db=ixomin^db,ixomax^db\}
6391 if(local_timestep)
then
6392 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
6398 ef(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
6399 ef(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
6400 ef(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6404 ef(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6411 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&ef(ix^d,^
c)**2+)*inv_squared_c)
6413 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&ef(ix^d,^
c)**2+)*inv_squared_c)
6417 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
6420 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
6421 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
6422 wprim(ix^d,
b^ce_)**2-ef(ix^d,^ce)**2*inv_squared_c+))
6426 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
6432 cot=1.d0/tan(x(ix^d,2))
6436 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
6437 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+ef(ix^d,1)*ef(ix^d,2)*inv_squared_c)
6439 if(.not.stagger_grid)
then
6440 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6442 tmp=tmp+wprim(ix^d,
psi_)*cot
6444 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6450 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6451 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+ef(ix^d,1)*ef(ix^d,2)*inv_squared_c&
6452 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6453 -wprim(ix^d,b3_)**2-ef(ix^d,3)**2*inv_squared_c)*cot)
6455 if(.not.stagger_grid)
then
6456 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6458 tmp=tmp+wprim(ix^d,
psi_)*cot
6460 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6463 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6464 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6465 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6466 +ef(ix^d,3)*ef(ix^d,1)*inv_squared_c&
6467 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6468 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6469 +ef(ix^d,2)*ef(ix^d,3)*inv_squared_c)*cot)
6471 if(.not.stagger_grid)
then
6472 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6473 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6474 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6475 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6476 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6483 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6486 end subroutine mhd_add_source_geom_semirelati
6495 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6501 integer,
intent(in) :: ixi^
l, ixo^
l
6502 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6503 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6505 double precision :: tmp,tmp1,tmp2,invr,cot
6506 double precision :: adiabs(ixi^s), gammas(ixi^s)
6508 integer :: mr_,mphi_
6509 integer :: br_,bphi_
6525 br_=mag(1); bphi_=mag(1)-1+
phi_
6530 {
do ix^db=ixomin^db,ixomax^db\}
6533 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6538 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6540 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6544 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6545 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6549 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6550 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6551 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6553 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(
block%B0(ix^
d,
phi_,0)*wprim(ix^
d,br_)+wprim(ix^
d,bphi_)*
block%B0(ix^
d,
r_,0))
6556 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6557 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6558 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6560 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6566 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6571 {
do ix^db=ixomin^db,ixomax^db\}
6573 if(local_timestep)
then
6574 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6578 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6579 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6582 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6583 if(b0field) w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp2*invr
6587 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6588 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6589 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6591 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6592 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6597 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6603 cot=1.d0/tan(x(ix^d,2))
6608 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6609 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6610 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6612 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6613 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6616 if(.not.stagger_grid)
then
6618 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6619 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6621 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6624 tmp=tmp+wprim(ix^d,
psi_)*cot
6626 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6632 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6633 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6634 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6635 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2-two*block%B0(ix^d,3,0)*wprim(ix^d,b3_))*cot)
6637 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6638 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6639 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6642 if(.not.stagger_grid)
then
6644 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6645 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6647 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6650 tmp=tmp+wprim(ix^d,
psi_)*cot
6652 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6656 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6657 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6658 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6659 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6660 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6661 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6662 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6663 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6664 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6666 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6667 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6668 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6669 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6670 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6673 if(.not.stagger_grid)
then
6675 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6676 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6677 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6678 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6679 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6680 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6681 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6682 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6683 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6685 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6686 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6687 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6688 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6689 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6697 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6700 end subroutine mhd_add_source_geom_split
6705 integer,
intent(in) :: ixi^
l, ixo^
l
6706 double precision,
intent(in) :: w(ixi^s, nw)
6707 double precision :: mge(ixo^s)
6710 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6712 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6716 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall,partial)
6720 integer,
intent(in) :: ixi^
l, ixo^
l
6721 double precision,
intent(in) :: w(ixi^s,nw)
6722 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6723 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6724 logical,
intent(in),
optional :: partial
6726 double precision :: current(ixi^s,7-2*
ndir:3)
6727 double precision :: rho(ixi^s)
6728 integer :: idir, idirmin, ix^
d
6729 logical :: use_partial
6732 if(
present(partial)) use_partial=partial
6734 if(.not.use_partial)
then
6745 do idir = idirmin,
ndir
6746 {
do ix^db=ixomin^db,ixomax^db\}
6747 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6751 end subroutine mhd_getv_hall
6753 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6756 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6757 double precision,
intent(in) :: qt
6758 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6759 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6762 double precision :: db(ixo^s), dpsi(ixo^s)
6766 {
do ix^db=ixomin^db,ixomax^db\}
6767 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6768 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6769 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6770 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6779 {
do ix^db=ixomin^db,ixomax^db\}
6780 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6781 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6782 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6783 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6784 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6786 if(total_energy)
then
6787 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6788 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6790 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6792 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6795 if(total_energy)
then
6796 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6797 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6802 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6804 end subroutine mhd_modify_wlr
6806 subroutine mhd_boundary_adjust(igrid,psb)
6808 integer,
intent(in) :: igrid
6811 integer :: ib, idims, iside, ixo^
l, i^
d
6820 i^
d=
kr(^
d,idims)*(2*iside-3);
6821 if (neighbor_type(i^
d,igrid)/=1) cycle
6822 ib=(idims-1)*2+iside
6840 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6845 end subroutine mhd_boundary_adjust
6847 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6850 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6851 double precision,
intent(inout) :: w(ixg^s,1:nw)
6852 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6854 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6855 integer :: ix^
d,ixf^
l
6861 if(total_energy)
call eos%to_primitive(ixg^
l,ixo^
l,w,x)
6869 do ix1=ixfmax1,ixfmin1,-1
6870 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6871 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6872 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6875 do ix1=ixfmax1,ixfmin1,-1
6876 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6877 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6878 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6879 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6880 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6881 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6882 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6896 do ix1=ixfmax1,ixfmin1,-1
6897 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6898 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6899 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6900 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6901 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6902 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6905 do ix1=ixfmax1,ixfmin1,-1
6906 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6907 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6908 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6909 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6910 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6911 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6912 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6913 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6914 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6915 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6916 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6917 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6918 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6919 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6920 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6921 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6922 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6923 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6928 if(total_energy)
call eos%to_conserved(ixg^
l,ixo^
l,w,x)
6931 if(total_energy)
call eos%to_primitive(ixg^
l,ixo^
l,w,x)
6939 do ix1=ixfmin1,ixfmax1
6940 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6941 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6942 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6945 do ix1=ixfmin1,ixfmax1
6946 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6947 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6948 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6949 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6950 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6951 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6952 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6966 do ix1=ixfmin1,ixfmax1
6967 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6968 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6969 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6970 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6971 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6972 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6975 do ix1=ixfmin1,ixfmax1
6976 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6977 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6978 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6979 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6980 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6981 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6982 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6983 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6984 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6985 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6986 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6987 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6988 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6989 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6990 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6991 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6992 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6993 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6998 if(total_energy)
call eos%to_conserved(ixg^
l,ixo^
l,w,x)
7001 if(total_energy)
call eos%to_primitive(ixg^
l,ixo^
l,w,x)
7009 do ix2=ixfmax2,ixfmin2,-1
7010 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
7011 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
7012 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
7015 do ix2=ixfmax2,ixfmin2,-1
7016 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
7017 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
7018 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
7019 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
7020 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
7021 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
7022 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
7036 do ix2=ixfmax2,ixfmin2,-1
7037 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
7038 ix2+1,ixfmin3:ixfmax3,mag(2)) &
7039 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
7040 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
7041 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
7042 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
7045 do ix2=ixfmax2,ixfmin2,-1
7046 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
7047 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
7048 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
7049 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
7050 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
7051 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
7052 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
7053 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
7054 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
7055 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
7056 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
7057 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
7058 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
7059 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
7060 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
7061 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
7062 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
7063 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
7068 if(total_energy)
call eos%to_conserved(ixg^
l,ixo^
l,w,x)
7071 if(total_energy)
call eos%to_primitive(ixg^
l,ixo^
l,w,x)
7079 do ix2=ixfmin2,ixfmax2
7080 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
7081 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
7082 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
7085 do ix2=ixfmin2,ixfmax2
7086 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
7087 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
7088 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
7089 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
7090 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
7091 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
7092 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
7106 do ix2=ixfmin2,ixfmax2
7107 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
7108 ix2-1,ixfmin3:ixfmax3,mag(2)) &
7109 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
7110 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
7111 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
7112 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
7115 do ix2=ixfmin2,ixfmax2
7116 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
7117 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
7118 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
7119 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
7120 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
7121 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
7122 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
7123 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
7124 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
7125 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
7126 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
7127 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
7128 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
7129 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
7130 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
7131 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
7132 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
7133 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
7138 if(total_energy)
call eos%to_conserved(ixg^
l,ixo^
l,w,x)
7142 if(total_energy)
call eos%to_primitive(ixg^
l,ixo^
l,w,x)
7152 do ix3=ixfmax3,ixfmin3,-1
7153 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
7154 ixfmin2:ixfmax2,ix3+1,mag(3)) &
7155 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
7156 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
7157 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
7158 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
7161 do ix3=ixfmax3,ixfmin3,-1
7162 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
7163 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
7164 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
7165 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
7166 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
7167 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7168 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
7169 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
7170 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7171 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
7172 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
7173 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
7174 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
7175 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
7176 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
7177 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
7178 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
7179 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
7183 if(total_energy)
call eos%to_conserved(ixg^
l,ixo^
l,w,x)
7186 if(total_energy)
call eos%to_primitive(ixg^
l,ixo^
l,w,x)
7196 do ix3=ixfmin3,ixfmax3
7197 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
7198 ixfmin2:ixfmax2,ix3-1,mag(3)) &
7199 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
7200 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
7201 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
7202 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
7205 do ix3=ixfmin3,ixfmax3
7206 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
7207 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
7208 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
7209 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
7210 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
7211 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7212 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
7213 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
7214 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
7215 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
7216 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
7217 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
7218 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
7219 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
7220 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
7221 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
7222 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
7223 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
7227 if(total_energy)
call eos%to_conserved(ixg^
l,ixo^
l,w,x)
7230 call mpistop(
"Special boundary is not defined for this region")
7233 end subroutine fixdivb_boundary
7242 double precision,
intent(in) :: qdt
7243 double precision,
intent(in) :: qt
7244 logical,
intent(inout) :: active
7247 integer,
parameter :: max_its = 50
7248 double precision :: residual_it(max_its), max_divb
7249 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
7250 double precision :: res
7251 double precision,
parameter :: max_residual = 1
d-3
7252 double precision,
parameter :: residual_reduction = 1
d-10
7253 integer :: iigrid, igrid
7254 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
7257 mg%operator_type = mg_laplacian
7265 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7266 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7269 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
7270 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7272 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7273 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7276 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7277 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7281 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
7282 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7283 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7291 do iigrid = 1, igridstail
7292 igrid = igrids(iigrid);
7295 lvl =
mg%boxes(id)%lvl
7296 nc =
mg%box_size_lvl(lvl)
7302 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
7304 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
7305 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
7310 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
7313 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
7314 if (
mype == 0) print *,
"iteration vs residual"
7317 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
7318 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
7319 if (residual_it(n) < residual_reduction * max_divb)
exit
7321 if (
mype == 0 .and. n > max_its)
then
7322 print *,
"divb_multigrid warning: not fully converged"
7323 print *,
"current amplitude of divb: ", residual_it(max_its)
7324 print *,
"multigrid smallest grid: ", &
7325 mg%domain_size_lvl(:,
mg%lowest_lvl)
7326 print *,
"note: smallest grid ideally has <= 8 cells"
7327 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
7328 print *,
"note: dx/dy/dz should be similar"
7332 call mg_fas_vcycle(
mg, max_res=res)
7333 if (res < max_residual)
exit
7335 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
7340 do iigrid = 1, igridstail
7341 igrid = igrids(iigrid);
7350 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
7354 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
7356 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
7358 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
7371 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
7372 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
7375 if(total_energy)
then
7377 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
7380 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
7390 subroutine mhd_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7394 integer,
intent(in) :: ixi^
l, ixo^
l
7395 double precision,
intent(in) :: qt,qdt
7397 double precision,
intent(in) :: wp(ixi^s,1:nw)
7398 type(state) :: sct, s
7399 type(ct_velocity) :: vcts
7400 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7401 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7403 double precision :: circ(ixi^s,1:
ndim)
7405 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7406 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
7407 integer :: idim1,idim2,idir,iwdim1,iwdim2
7409 associate(bfaces=>s%ws,x=>s%x)
7416 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7423 i1kr^
d=
kr(idim1,^
d);
7426 i2kr^
d=
kr(idim2,^
d);
7429 if (
lvc(idim1,idim2,idir)==1)
then
7431 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7433 {
do ix^db=ixcmin^db,ixcmax^db\}
7434 fe(ix^
d,idir)=quarter*&
7435 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
7436 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
7438 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
7443 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
7451 if(
associated(usr_set_electric_field)) &
7452 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7454 circ(ixi^s,1:ndim)=zero
7459 ixcmin^d=ixomin^d-kr(idim1,^d);
7461 ixa^l=ixc^l-kr(idim2,^d);
7464 if(lvc(idim1,idim2,idir)==1)
then
7466 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7469 else if(lvc(idim1,idim2,idir)==-1)
then
7471 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7477 {
do ix^db=ixcmin^db,ixcmax^db\}
7479 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7481 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7488 end subroutine mhd_update_faces_average
7491 subroutine mhd_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7496 integer,
intent(in) :: ixi^
l, ixo^
l
7497 double precision,
intent(in) :: qt, qdt
7499 double precision,
intent(in) :: wp(ixi^s,1:nw)
7500 type(state) :: sct, s
7501 type(ct_velocity) :: vcts
7502 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7503 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7505 double precision :: circ(ixi^s,1:
ndim)
7507 double precision :: ecc(ixi^s,
sdim:3)
7508 double precision :: ein(ixi^s,
sdim:3)
7510 double precision :: el(ixi^s),er(ixi^s)
7512 double precision :: elc,erc
7514 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7516 double precision :: jce(ixi^s,
sdim:3)
7518 double precision :: xs(ixgs^t,1:
ndim)
7519 double precision :: gradi(ixgs^t)
7520 integer :: ixc^
l,ixa^
l
7521 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7523 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7526 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7532 {
do ix^db=iximin^db,iximax^db\}
7535 ecc(ix^
d,1)=(wp(ix^
d,b2_)+
block%B0(ix^
d,2,0))*wp(ix^
d,m3_)-(wp(ix^
d,b3_)+
block%B0(ix^
d,3,0))*wp(ix^
d,m2_)
7536 ecc(ix^
d,2)=(wp(ix^
d,b3_)+
block%B0(ix^
d,3,0))*wp(ix^
d,m1_)-(wp(ix^
d,b1_)+
block%B0(ix^
d,1,0))*wp(ix^
d,m3_)
7537 ecc(ix^
d,3)=(wp(ix^
d,b1_)+
block%B0(ix^
d,1,0))*wp(ix^
d,m2_)-(wp(ix^
d,b2_)+
block%B0(ix^
d,2,0))*wp(ix^
d,m1_)
7540 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7547 {
do ix^db=iximin^db,iximax^db\}
7550 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7551 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7552 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7555 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7569 i1kr^d=kr(idim1,^d);
7572 i2kr^d=kr(idim2,^d);
7575 if (lvc(idim1,idim2,idir)==1)
then
7577 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7580 {
do ix^db=ixcmin^db,ixcmax^db\}
7581 fe(ix^d,idir)=quarter*&
7582 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7583 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7588 ixamax^d=ixcmax^d+i1kr^d;
7589 {
do ix^db=ixamin^db,ixamax^db\}
7590 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7591 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7594 do ix^db=ixcmin^db,ixcmax^db\}
7595 if(vnorm(ix^d,idim1)>0.d0)
then
7597 else if(vnorm(ix^d,idim1)<0.d0)
then
7598 elc=el({ix^d+i1kr^d})
7600 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7602 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7604 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7605 erc=er({ix^d+i1kr^d})
7607 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7609 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7614 ixamax^d=ixcmax^d+i2kr^d;
7615 {
do ix^db=ixamin^db,ixamax^db\}
7616 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7617 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7620 do ix^db=ixcmin^db,ixcmax^db\}
7621 if(vnorm(ix^d,idim2)>0.d0)
then
7623 else if(vnorm(ix^d,idim2)<0.d0)
then
7624 elc=el({ix^d+i2kr^d})
7626 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7628 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7630 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7631 erc=er({ix^d+i2kr^d})
7633 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7635 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7639 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7644 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7658 if (lvc(idim1,idim2,idir)==0) cycle
7660 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7661 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7664 xs(ixa^s,:)=x(ixa^s,:)
7665 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7666 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7667 if (lvc(idim1,idim2,idir)==1)
then
7668 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7670 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7677 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7679 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7683 {
do ix^db=ixomin^db,ixomax^db\}
7684 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7685 +ein(ix1,ix2-1,ix3-1,idir))
7686 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7687 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7689 else if(idir==2)
then
7690 {
do ix^db=ixomin^db,ixomax^db\}
7691 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7692 +ein(ix1-1,ix2,ix3-1,idir))
7693 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7694 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7697 {
do ix^db=ixomin^db,ixomax^db\}
7698 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7699 +ein(ix1-1,ix2-1,ix3,idir))
7700 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7701 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7707 {
do ix^db=ixomin^db,ixomax^db\}
7708 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7709 +ein(ix1-1,ix2-1,idir))
7710 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7711 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7722 if(
associated(usr_set_electric_field)) &
7723 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7725 circ(ixi^s,1:ndim)=zero
7730 ixcmin^d=ixomin^d-kr(idim1,^d);
7732 ixa^l=ixc^l-kr(idim2,^d);
7735 if(lvc(idim1,idim2,idir)==1)
then
7737 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7740 else if(lvc(idim1,idim2,idir)==-1)
then
7742 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7748 {
do ix^db=ixcmin^db,ixcmax^db\}
7750 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7752 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7759 end subroutine mhd_update_faces_contact
7762 subroutine mhd_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7767 integer,
intent(in) :: ixi^
l, ixo^
l
7768 double precision,
intent(in) :: qt, qdt
7770 double precision,
intent(in) :: wp(ixi^s,1:nw)
7771 type(state) :: sct, s
7772 type(ct_velocity) :: vcts
7773 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7774 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7776 double precision :: vtill(ixi^s,2)
7777 double precision :: vtilr(ixi^s,2)
7778 double precision :: bfacetot(ixi^s,
ndim)
7779 double precision :: btill(ixi^s,
ndim)
7780 double precision :: btilr(ixi^s,
ndim)
7781 double precision :: cp(ixi^s,2)
7782 double precision :: cm(ixi^s,2)
7783 double precision :: circ(ixi^s,1:
ndim)
7785 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7786 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7787 integer :: idim1,idim2,idir,ix^
d
7789 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7790 cbarmax=>vcts%cbarmax)
7803 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7819 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7823 idim2=mod(idir+1,3)+1
7825 jxc^
l=ixc^
l+
kr(idim1,^
d);
7826 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7830 vtill(ixi^s,2),vtilr(ixi^s,2))
7833 vtill(ixi^s,1),vtilr(ixi^s,1))
7839 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7840 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7842 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7843 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7846 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7849 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7853 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7854 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7856 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7857 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7861 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7862 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7863 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7864 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7865 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7866 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7867 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7868 /(cp(ixc^s,2)+cm(ixc^s,2))
7871 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7875 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7889 circ(ixi^s,1:
ndim)=zero
7894 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7898 if(
lvc(idim1,idim2,idir)/=0)
then
7899 hxc^
l=ixc^
l-
kr(idim2,^
d);
7901 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7902 +
lvc(idim1,idim2,idir)&
7908 {
do ix^db=ixcmin^db,ixcmax^db\}
7910 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7912 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7918 end subroutine mhd_update_faces_hll
7921 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
7926 integer,
intent(in) :: ixi^
l, ixo^
l
7928 double precision,
intent(in) :: wp(ixi^s,1:nw)
7929 type(state),
intent(in) :: sct, s
7931 double precision :: jce(ixi^s,
sdim:3)
7934 double precision :: jcc(ixi^s,7-2*
ndir:3)
7936 double precision :: xs(ixgs^t,1:
ndim)
7938 double precision :: eta(ixi^s)
7939 double precision :: gradi(ixgs^t)
7940 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7942 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7948 if (
lvc(idim1,idim2,idir)==0) cycle
7950 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7951 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7954 xs(ixb^s,:)=x(ixb^s,:)
7955 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7956 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7957 if (
lvc(idim1,idim2,idir)==1)
then
7958 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7960 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7967 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7975 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7976 jcc(ixc^s,idir)=0.d0
7978 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7979 ixamin^
d=ixcmin^
d+ix^
d;
7980 ixamax^
d=ixcmax^
d+ix^
d;
7981 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7983 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7984 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7989 end subroutine get_resistive_electric_field
7992 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7995 integer,
intent(in) :: ixi^
l, ixo^
l
7996 double precision,
intent(in) :: w(ixi^s,1:nw)
7997 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7998 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
8000 double precision :: jxbxb(ixi^s,1:3)
8001 integer :: idir,ixa^
l,ixc^
l,ix^
d
8004 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
8011 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
8014 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
8015 ixamin^
d=ixcmin^
d+ix^
d;
8016 ixamax^
d=ixcmax^
d+ix^
d;
8017 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
8019 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
8022 end subroutine get_ambipolar_electric_field
8028 integer,
intent(in) :: ixo^
l
8038 do ix^db=ixomin^db,ixomax^db\}
8040 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
8041 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
8042 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
8043 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
8044 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
8045 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
8048 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
8049 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
8050 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
8051 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
8094 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
8095 double precision,
intent(inout) :: ws(ixis^s,1:nws)
8096 double precision,
intent(in) :: x(ixi^s,1:
ndim)
8098 double precision :: adummy(ixis^s,1:3)
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
subroutine cak_init(phys_gamma)
Initialize the module.
subroutine cak_get_dt(wprim, ixil, ixol, dtnew, dxd, x)
Check time step for total radiation contribution.
subroutine cak_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
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
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...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
PI (partial-ionisation) ionisation-degree backend for the eos% family.
Equation of state for AMRVAC, handled through a single eos_container object.
Module for flux conservation near refinement boundaries.
subroutine, public store_flux(igrid, fc, idimlim, nwfluxin)
subroutine, public store_edge(igrid, ixil, fe, idimlim)
Module for flux limited diffusion (FLD)-approximation in Radiation-(Magneto)hydrodynamics simulations...
double precision, public fld_bisect_tol
Tolerance for bisection method for Energy sourceterms This is a percentage of the minimum of gas- and...
subroutine, public fld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
get dt limit for radiation force and FLD explicit source additions NOTE: w is primitive on entry
double precision, public fld_diff_tol
Tolerance for radiative Energy diffusion.
character(len=40) fld_fluxlimiter
flux limiter choice
character(len=40) fld_opal_table
double precision, public fld_kappa0
Opacity value when using constant opacity.
subroutine, public add_fld_rad_force(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO This subroutine handles th...
character(len=40) fld_opacity_law
switches for opacity
character(len=40) fld_interaction_method
Which method to find the root for the energy interaction polynomial.
subroutine, public fld_get_radpress(w, x, ixil, ixol, rad_pressure, fl)
Returns Radiation Pressure as tensor NOTE: w is primitive on entry.
logical fld_radforce_split
source split for energy interact and radforce:
subroutine, public fld_implicit_update(dtfactor, qdt, qtc, psa, psb, fl)
Calling all subroutines to perform the multigrid method Communicates rad_e and diff_coeff to multigri...
subroutine, public fld_evaluate_implicit(qtc, psa, fl)
inplace update of psa==>F_im(psa)
subroutine, public fld_init()
Initialising FLD-module Read opacities Initialise Multigrid and adimensionalise kappa.
integer nth_for_diff_mg
diffusion coefficient stencil control
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, nth_in)
subroutine laplacian_of_vector(qvec, ixil, ixol, lapl_qvec)
integer, parameter cylindrical
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...
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
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.
double precision const_kappae
double precision arad_norm
Normalised radiation constant.
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.
pure subroutine cross_product(ixil, ixol, a, b, axb)
Cross product of two vectors.
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.
double precision unit_opacity
Physical scaling factor for Opacity.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision unit_mass
Physical scaling factor for mass.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
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.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision const_rad_a
Physical factors useful for radiation fld.
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.
logical local_timestep
each cell has its own timestep or not
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
double precision, dimension(:), allocatable, parameter d
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 const_sigmasb
integer nwauxio
Number of auxiliary variables that are only included in output.
double precision unit_velocity
Physical scaling factor for velocity.
double precision small_r_e
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.
double precision unit_radflux
Physical scaling factor for radiation flux.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
spatial steps for all dimensions at all levels
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
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 max_blocks
The maximum number of grid blocks in a processor.
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
integer phys_trac_finegrid
integer, parameter unitconvert
double precision unit_erad
Physical scaling factor for radiation energy density.
integer number_equi_vars
number of equilibrium set variables, besides the mag field
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for including gravity in (magneto)hydrodynamics simulations.
subroutine gravity_get_dt(wprim, ixil, ixol, dtnew, dxd, x)
subroutine gravity_init()
Initialize the module.
subroutine gravity_add_source(qdt, ixil, ixol, wct, wctprim, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
subroutine, public mhd_get_trad(w, x, ixil, ixol, trad)
Calculates radiation temperature.
integer, public, protected c_
logical, public, protected mhd_gravity
Whether gravity is added.
integer, public, protected fip_
Index of the FIP passive scalar rho*fip in conserved form, fip in primitive form.
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
logical, public, protected mhd_glm_extended
Whether extended GLM-MHD is used with additional sources.
logical, public mhd_hyperbolic_tc_constant
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
integer, public, protected qpar_
Index of the field-aligned heat flux q_parallel.
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
subroutine, public mhd_e_to_ei(ixil, ixol, w, x)
Transform total energy to internal energy.
double precision, public mhd_adiab
The adiabatic constant.
double precision, public divbdiff
Coefficient of diffusive divB cleaning.
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
double precision, public, protected mhd_hyperbolic_tc_bmin
Field-strength transition scale for perpendicular closure.
double precision, public, protected rr
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
integer, public, protected mhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
logical, public numerical_resistive_heating
Whether numerical resistive heating is included when solving partial energy equation.
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
logical, public, protected mhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected mhd_semirelativistic
Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
integer, public, protected mhd_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public, protected mhd_n_tracer
Number of tracer species.
integer, public, protected te_
Indices of temperature.
integer, public, protected m
integer, public equi_rho0_
equi vars indices in the stateequi_vars array
integer, public, protected mhd_trac_type
Which TRAC method is used.
logical, public, protected mhd_cak_force
Whether plasma is partially ionized Whether CAK radiation line force is activated.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
integer, public, protected ne_
Index of the electron number density for LTE module.
logical, public, protected mhd_hall
Whether Hall-MHD is used.
type(te_fluid), allocatable, public te_fl_mhd
type of fluid for thermal emission synthesis
logical, public, protected mhd_ambipolar
Whether Ambipolar term is used.
logical, public, protected mhd_hyperbolic_tc
Whether thermal conduction is used.
logical, public, protected mhd_hyperbolic_tc_sat
Whether saturation is considered for hyperbolic TC.
double precision, public, protected mhd_hyperbolic_tc_kappa_perp_factor
Relative perpendicular hyperbolic-TC coefficient in fixed/strong-field limit: kappa_perp0 = mhd_hyper...
logical, public has_equi_rho_and_p
whether split off equilibrium density and pressure
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixil, ixol)
Compute 2 times total magnetic energy.
logical, public, protected mhd_radiation_fld
Whether radiation-gas interaction is handled using flux limited diffusion.
subroutine, public multiplyambicoef(ixil, ixol, res, w, x)
multiply res by the ambipolar coefficient The ambipolar coefficient is calculated as -mhd_eta_ambi/rh...
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
logical, public, protected mhd_viscosity
Whether viscosity is added.
subroutine, public mhd_get_pradiation_from_prim(w, x, ixil, ixol, prad)
Calculate radiation pressure within ixO^L.
double precision, public, protected mhd_reduced_c
Reduced speed of light for semirelativistic MHD: 2% of light speed.
logical, public, protected mhd_energy
Whether an energy equation is used.
logical, public, protected mhd_ambipolar_exp
Whether Ambipolar term is implemented explicitly.
double precision, public mhd_hyperbolic_tc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
logical, public, protected mhd_glm
Whether GLM-MHD is used to control div B.
type(fld_fluid), allocatable, public fld_fl
Radiation fluid object (gas-EoS callbacks for FLD), wired in mhd_link_eos.
logical, public clean_initial_divb
clean initial divB
double precision, public mhd_eta
The MHD resistivity.
logical, public divbwave
Add divB wave in Roe solver.
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
double precision, public, protected mhd_trac_mask
Height of the mask used in the TRAC method.
procedure(mask_subroutine), pointer, public usr_mask_ambipolar
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
integer, public equi_pe0_
subroutine, public mhd_get_csrad2_prim(w, x, ixil, ixol, csound)
Calculate modified squared fast wave speed for FLD NOTE: w is primitive on entry here!...
integer, public, protected qperp_
Index of the perpendicular heat flux q_perp.
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
logical, public, protected mhd_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
logical, public, protected mhd_particles
Whether particles module is added.
integer, public, protected b
subroutine, public mhd_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
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)...
double precision, public mhd_etah
Hall resistivity.
integer, public, protected mhd_hyperbolic_tc_perp_mode
Perpendicular hyperbolic-TC closure mode: 0 = off, 1 = fixed anisotropy, 2 = field-strength-dependent...
subroutine, public mhd_get_v(w, x, ixil, ixol, v)
Calculate v vector.
double precision, public mhd_eta_ambi
The MHD ambipolar coefficient.
logical, public, protected mhd_fip
Whether FIP passive scalar is enabled.
logical, public, protected mhd_hydrodynamic_e
Whether hydrodynamic energy is solved instead of total energy.
integer, public, protected r_e
Index of the radiation energy.
subroutine, public mhd_phys_init()
logical, public, protected mhd_trac
Whether TRAC method is used.
subroutine, public mhd_get_csrad2(w, x, ixil, ixol, csound)
Calculate modified squared sound speed for FLD NOTE: only for diagnostic purposes,...
subroutine, public mhd_get_pthermal_plus_pradiation(w, x, ixil, ixol, pth_plus_prad)
Calculates the sum of the gas pressure and the max Prad tensor element.
type(rc_fluid), allocatable, public rc_fl
type of fluid for radiative cooling
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
integer, public, protected rho_
Index of the density (in the w array)
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
integer, public, protected tweight_
logical, public, protected mhd_ambipolar_sts
Whether Ambipolar term is implemented using supertimestepping.
logical, public, protected mhd_hyperbolic_tc_use_perp
Whether the perpendicular hyperbolic-TC channel is enabled.
subroutine, public mhd_ei_to_e(ixil, ixol, w, x)
Transform internal energy to total energy.
integer, public, protected e_
Index of the energy density (-1 if not present)
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine, public mhd_get_rho(w, x, ixil, ixol, rho)
integer, public, protected psi_
Indices of the GLM psi.
logical, public mhd_equi_thermal
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.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_add_source(qdt, dtfactor, ixil, ixol, wct, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rotating_frame_init()
Initialize the module.
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method which can be used for multiple source terms in the governing equatio...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (mhd implementation) Note: for multi-D MHD (1D MHD will use HD f...
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation) Note: also used in 1D MHD (or for neutrals i...
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
double precision time_htc_total
double precision time_htc0
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
procedure(special_resistivity), pointer usr_special_resistivity
procedure(set_adiab), pointer usr_set_adiab
procedure(set_adiab), pointer usr_set_gamma
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
The module add viscous source terms and check time step.
subroutine, public viscosity_get_dt(wprim, ixil, ixol, dtnew, dxd, x)
procedure(sub_add_source), pointer, public viscosity_add_source
subroutine, public viscosity_init(phys_wider_stencil)
Initialize the module.
Radiation fluid object: gas-EoS callbacks the FLD module needs, wired by the physics module at link t...
The data structure that contains information about a tree node/grid block.