22 double precision,
public ::
mhd_eta = 0.0d0
30 double precision,
protected :: small_e
41 double precision :: divbdiff = 0.8d0
46 double precision,
public,
protected ::
h_ion_fr=1d0
49 double precision,
public,
protected ::
he_ion_fr=1d0
56 double precision,
public,
protected ::
rr=1d0
58 double precision :: gamma_1, inv_gamma_1
60 double precision :: inv_squared_c0, inv_squared_c
67 integer,
public,
protected ::
rho_
69 integer,
allocatable,
public,
protected ::
mom(:)
71 integer,
public,
protected :: ^
c&m^C_
73 integer,
public,
protected ::
e_
75 integer,
public,
protected :: ^
c&b^C_
77 integer,
public,
protected ::
p_
79 integer,
public,
protected ::
q_
81 integer,
public,
protected ::
psi_
83 integer,
public,
protected ::
te_
88 integer,
allocatable,
public,
protected ::
tracer(:)
96 integer,
parameter :: divb_none = 0
97 integer,
parameter :: divb_multigrid = -1
98 integer,
parameter :: divb_glm = 1
99 integer,
parameter :: divb_powel = 2
100 integer,
parameter :: divb_janhunen = 3
101 integer,
parameter :: divb_linde = 4
102 integer,
parameter :: divb_lindejanhunen = 5
103 integer,
parameter :: divb_lindepowel = 6
104 integer,
parameter :: divb_lindeglm = 7
105 integer,
parameter :: divb_ct = 8
135 logical,
public,
protected ::
mhd_glm = .false.
170 logical :: compactres = .false.
184 logical :: total_energy = .true.
188 logical :: gravity_energy
190 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
192 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
194 character(len=std_len) :: typedivbdiff =
'all'
205 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
207 integer,
intent(in) :: ixi^
l, ixo^
l
208 double precision,
intent(in) :: x(ixi^s,1:
ndim)
209 double precision,
intent(in) :: w(ixi^s,1:nw)
210 double precision,
intent(inout) :: res(ixi^s)
211 end subroutine mask_subroutine
249 subroutine mhd_read_params(files)
252 character(len=*),
intent(in) :: files(:)
268 do n = 1,
size(files)
269 open(
unitpar, file=trim(files(n)), status=
"old")
270 read(
unitpar, mhd_list,
end=111)
274 end subroutine mhd_read_params
277 subroutine mhd_write_info(fh)
279 integer,
intent(in) :: fh
282 integer,
parameter :: n_par = 1
283 double precision :: values(n_par)
284 integer,
dimension(MPI_STATUS_SIZE) :: st
285 character(len=name_len) :: names(n_par)
287 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
291 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
292 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
293 end subroutine mhd_write_info
320 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
331 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
335 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
339 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
343 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F'
347 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
351 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
355 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
359 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
363 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho0=F when mhd_energy=F'
367 if(
mype==0)
write(*,*)
'WARNING: set has_equi_pe0=F when mhd_energy=F'
373 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
379 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
403 phys_total_energy=total_energy
406 gravity_energy=.false.
408 gravity_energy=.true.
411 gravity_energy=.false.
417 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
422 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
431 type_divb = divb_none
434 type_divb = divb_multigrid
436 mg%operator_type = mg_laplacian
443 case (
'powel',
'powell')
444 type_divb = divb_powel
446 type_divb = divb_janhunen
448 type_divb = divb_linde
449 case (
'lindejanhunen')
450 type_divb = divb_lindejanhunen
452 type_divb = divb_lindepowel
456 type_divb = divb_lindeglm
461 call mpistop(
'Unknown divB fix')
466 allocate(start_indices(number_species),stop_indices(number_species))
473 mom(:) = var_set_momentum(
ndir)
479 e_ = var_set_energy()
488 mag(:) = var_set_bfield(
ndir)
492 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
508 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
513 te_ = var_set_auxvar(
'Te',
'Te')
522 stop_indices(1)=nwflux
553 allocate(iw_vector(nvector))
554 iw_vector(1) =
mom(1) - 1
555 iw_vector(2) = mag(1) - 1
558 if (.not.
allocated(flux_type))
then
559 allocate(flux_type(
ndir, nwflux))
560 flux_type = flux_default
561 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
562 call mpistop(
"phys_check error: flux_type has wrong shape")
565 if(nwflux>mag(
ndir))
then
567 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
572 flux_type(:,
psi_)=flux_special
574 flux_type(idir,mag(idir))=flux_special
578 flux_type(idir,mag(idir))=flux_tvdlf
584 phys_get_dt => mhd_get_dt
587 phys_get_cmax => mhd_get_cmax_semirelati
589 phys_get_cmax => mhd_get_cmax_semirelati_noe
593 phys_get_cmax => mhd_get_cmax_origin
595 phys_get_cmax => mhd_get_cmax_origin_noe
598 phys_get_a2max => mhd_get_a2max
599 phys_get_tcutoff => mhd_get_tcutoff
600 phys_get_h_speed => mhd_get_h_speed
602 phys_get_cbounds => mhd_get_cbounds_split_rho
604 phys_get_cbounds => mhd_get_cbounds_semirelati
606 phys_get_cbounds => mhd_get_cbounds
609 phys_to_primitive => mhd_to_primitive_hde
611 phys_to_conserved => mhd_to_conserved_hde
615 phys_to_primitive => mhd_to_primitive_semirelati
617 phys_to_conserved => mhd_to_conserved_semirelati
620 phys_to_primitive => mhd_to_primitive_semirelati_noe
622 phys_to_conserved => mhd_to_conserved_semirelati_noe
627 phys_to_primitive => mhd_to_primitive_split_rho
629 phys_to_conserved => mhd_to_conserved_split_rho
632 phys_to_primitive => mhd_to_primitive_inte
634 phys_to_conserved => mhd_to_conserved_inte
637 phys_to_primitive => mhd_to_primitive_origin
639 phys_to_conserved => mhd_to_conserved_origin
642 phys_to_primitive => mhd_to_primitive_origin_noe
644 phys_to_conserved => mhd_to_conserved_origin_noe
649 phys_get_flux => mhd_get_flux_hde
652 phys_get_flux => mhd_get_flux_semirelati
654 phys_get_flux => mhd_get_flux_semirelati_noe
658 phys_get_flux => mhd_get_flux_split
660 phys_get_flux => mhd_get_flux
662 phys_get_flux => mhd_get_flux_noe
667 phys_add_source_geom => mhd_add_source_geom_semirelati
669 phys_add_source_geom => mhd_add_source_geom_split
671 phys_add_source_geom => mhd_add_source_geom
673 phys_add_source => mhd_add_source
674 phys_check_params => mhd_check_params
675 phys_write_info => mhd_write_info
678 phys_handle_small_values => mhd_handle_small_values_inte
679 mhd_handle_small_values => mhd_handle_small_values_inte
680 phys_check_w => mhd_check_w_inte
682 phys_handle_small_values => mhd_handle_small_values_hde
683 mhd_handle_small_values => mhd_handle_small_values_hde
684 phys_check_w => mhd_check_w_hde
686 phys_handle_small_values => mhd_handle_small_values_semirelati
687 mhd_handle_small_values => mhd_handle_small_values_semirelati
688 phys_check_w => mhd_check_w_semirelati
690 phys_handle_small_values => mhd_handle_small_values_split
691 mhd_handle_small_values => mhd_handle_small_values_split
692 phys_check_w => mhd_check_w_split
694 phys_handle_small_values => mhd_handle_small_values_origin
695 mhd_handle_small_values => mhd_handle_small_values_origin
696 phys_check_w => mhd_check_w_origin
698 phys_handle_small_values => mhd_handle_small_values_noe
699 mhd_handle_small_values => mhd_handle_small_values_noe
700 phys_check_w => mhd_check_w_noe
704 phys_get_pthermal => mhd_get_pthermal_inte
707 phys_get_pthermal => mhd_get_pthermal_hde
710 phys_get_pthermal => mhd_get_pthermal_semirelati
713 phys_get_pthermal => mhd_get_pthermal_origin
716 phys_get_pthermal => mhd_get_pthermal_noe
721 phys_set_equi_vars => set_equi_vars_grid
724 if(type_divb==divb_glm)
then
725 phys_modify_wlr => mhd_modify_wlr
731 phys_update_temperature => mhd_update_temperature
756 transverse_ghost_cells = 1
757 phys_get_ct_velocity => mhd_get_ct_velocity_average
758 phys_update_faces => mhd_update_faces_average
760 transverse_ghost_cells = 1
761 phys_get_ct_velocity => mhd_get_ct_velocity_contact
762 phys_update_faces => mhd_update_faces_contact
764 transverse_ghost_cells = 2
765 phys_get_ct_velocity => mhd_get_ct_velocity_hll
766 phys_update_faces => mhd_update_faces_hll
768 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
771 phys_modify_wlr => mhd_modify_wlr
773 phys_boundary_adjust => mhd_boundary_adjust
782 call mhd_physical_units()
788 call mpistop(
"thermal conduction needs mhd_energy=T")
791 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
794 call mpistop(
"radiative cooling needs mhd_energy=T")
809 if(phys_internal_e)
then
811 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
813 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
816 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
819 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
821 tc_fl%has_equi = .true.
822 tc_fl%get_temperature_equi => mhd_get_temperature_equi
823 tc_fl%get_rho_equi => mhd_get_rho_equi
825 tc_fl%has_equi = .false.
828 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
856 rc_fl%has_equi = .true.
857 rc_fl%get_rho_equi => mhd_get_rho_equi
858 rc_fl%get_pthermal_equi => mhd_get_pe_equi
860 rc_fl%has_equi = .false.
868 phys_te_images => mhd_te_images
884 if (particles_eta < zero) particles_eta =
mhd_eta
885 if (particles_etah < zero) particles_eta =
mhd_etah
887 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
888 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
902 phys_wider_stencil = 2
904 phys_wider_stencil = 1
912 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
924 phys_wider_stencil = 2
926 phys_wider_stencil = 1
940 subroutine mhd_te_images
945 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
947 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
949 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
951 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
954 call mpistop(
"Error in synthesize emission: Unknown convert_type")
956 end subroutine mhd_te_images
962 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
966 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
967 double precision,
intent(in) :: x(ixi^s,1:
ndim)
968 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
969 double precision,
intent(in) :: my_dt
970 logical,
intent(in) :: fix_conserve_at_step
972 end subroutine mhd_sts_set_source_tc_mhd
974 subroutine mhd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
978 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
979 double precision,
intent(in) :: x(ixi^s,1:
ndim)
980 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
981 double precision,
intent(in) :: my_dt
982 logical,
intent(in) :: fix_conserve_at_step
984 end subroutine mhd_sts_set_source_tc_hd
986 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
993 integer,
intent(in) :: ixi^
l, ixo^
l
994 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
995 double precision,
intent(in) :: w(ixi^s,1:nw)
996 double precision :: dtnew
999 end function mhd_get_tc_dt_mhd
1001 function mhd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1008 integer,
intent(in) :: ixi^
l, ixo^
l
1009 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1010 double precision,
intent(in) :: w(ixi^s,1:nw)
1011 double precision :: dtnew
1014 end function mhd_get_tc_dt_hd
1016 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
1019 integer,
intent(in) :: ixi^
l,ixo^
l
1020 double precision,
intent(inout) :: w(ixi^s,1:nw)
1021 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1022 integer,
intent(in) :: step
1023 character(len=140) :: error_msg
1025 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1026 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1027 end subroutine mhd_tc_handle_small_e
1030 subroutine tc_params_read_mhd(fl)
1032 type(tc_fluid),
intent(inout) :: fl
1034 double precision :: tc_k_para=0d0
1035 double precision :: tc_k_perp=0d0
1038 logical :: tc_perpendicular=.false.
1039 logical :: tc_saturate=.false.
1040 character(len=std_len) :: tc_slope_limiter=
"MC"
1042 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1046 read(
unitpar, tc_list,
end=111)
1050 fl%tc_perpendicular = tc_perpendicular
1051 fl%tc_saturate = tc_saturate
1052 fl%tc_k_para = tc_k_para
1053 fl%tc_k_perp = tc_k_perp
1054 select case(tc_slope_limiter)
1056 fl%tc_slope_limiter = 0
1059 fl%tc_slope_limiter = 1
1062 fl%tc_slope_limiter = 2
1065 fl%tc_slope_limiter = 3
1068 fl%tc_slope_limiter = 4
1070 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1072 end subroutine tc_params_read_mhd
1076 subroutine rc_params_read(fl)
1079 type(rc_fluid),
intent(inout) :: fl
1081 double precision :: cfrac=0.1d0
1084 double precision :: rad_cut_hgt=0.5d0
1085 double precision :: rad_cut_dey=0.15d0
1088 integer :: ncool = 4000
1090 logical :: tfix=.false.
1092 logical :: rc_split=.false.
1093 logical :: rad_cut=.false.
1095 character(len=std_len) :: coolcurve=
'JCcorona'
1097 character(len=std_len) :: coolmethod=
'exact'
1099 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1103 read(
unitpar, rc_list,
end=111)
1108 fl%coolcurve=coolcurve
1109 fl%coolmethod=coolmethod
1112 fl%rc_split=rc_split
1115 fl%rad_cut_hgt=rad_cut_hgt
1116 fl%rad_cut_dey=rad_cut_dey
1117 end subroutine rc_params_read
1121 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1124 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1125 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1127 double precision :: delx(ixi^s,1:
ndim)
1128 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1129 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1135 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1139 hxo^
l=ixo^
l-
kr(idims,^
d);
1145 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1148 xshift^
d=half*(one-
kr(^
d,idims));
1155 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1159 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1162 end subroutine set_equi_vars_grid_faces
1165 subroutine set_equi_vars_grid(igrid)
1169 integer,
intent(in) :: igrid
1175 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1177 end subroutine set_equi_vars_grid
1180 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1182 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1183 double precision,
intent(in) :: w(ixi^s, 1:nw)
1184 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1185 double precision :: wnew(ixo^s, 1:nwc)
1192 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1198 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1202 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1206 if(
b0field .and. total_energy)
then
1207 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1208 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1212 end function convert_vars_splitting
1214 subroutine mhd_check_params
1222 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1223 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1227 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1228 inv_gamma_1=1.d0/gamma_1
1233 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1238 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1243 end subroutine mhd_check_params
1245 subroutine mhd_physical_units()
1247 double precision :: mp,kb,miu0,c_lightspeed
1248 double precision :: a,
b
1259 c_lightspeed=const_c
1400 end subroutine mhd_physical_units
1402 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1405 logical,
intent(in) :: primitive
1406 logical,
intent(inout) :: flag(ixi^s,1:nw)
1407 integer,
intent(in) :: ixi^
l, ixo^
l
1408 double precision,
intent(in) :: w(ixi^s,nw)
1410 double precision :: tmp,b2,
b(ixo^s,1:
ndir)
1411 double precision :: v(ixo^s,1:
ndir),gamma2,inv_rho
1422 {
do ix^db=ixomin^db,ixomax^db \}
1423 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1426 {
do ix^db=ixomin^db,ixomax^db \}
1427 b2=(^
c&w(ix^d,
b^
c_)**2+)
1428 if(b2>smalldouble)
then
1433 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1434 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1435 inv_rho = 1d0/w(ix^d,
rho_)
1437 b2=b2*inv_rho*inv_squared_c
1439 gamma2=1.d0/(1.d0+b2)
1441 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1444 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1445 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1446 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1451 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1457 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1458 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1459 if(tmp<small_e) flag(ix^d,
e_)=.true.
1465 end subroutine mhd_check_w_semirelati
1467 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1470 logical,
intent(in) :: primitive
1471 integer,
intent(in) :: ixi^
l, ixo^
l
1472 double precision,
intent(in) :: w(ixi^s,nw)
1473 logical,
intent(inout) :: flag(ixi^s,1:nw)
1478 {
do ix^db=ixomin^db,ixomax^db\}
1484 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1488 end subroutine mhd_check_w_origin
1490 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1493 logical,
intent(in) :: primitive
1494 integer,
intent(in) :: ixi^
l, ixo^
l
1495 double precision,
intent(in) :: w(ixi^s,nw)
1496 logical,
intent(inout) :: flag(ixi^s,1:nw)
1498 double precision :: tmp
1502 {
do ix^db=ixomin^db,ixomax^db\}
1508 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1513 end subroutine mhd_check_w_split
1515 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1518 logical,
intent(in) :: primitive
1519 integer,
intent(in) :: ixi^
l, ixo^
l
1520 double precision,
intent(in) :: w(ixi^s,nw)
1521 logical,
intent(inout) :: flag(ixi^s,1:nw)
1526 {
do ix^db=ixomin^db,ixomax^db\}
1530 end subroutine mhd_check_w_noe
1532 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1535 logical,
intent(in) :: primitive
1536 integer,
intent(in) :: ixi^
l, ixo^
l
1537 double precision,
intent(in) :: w(ixi^s,nw)
1538 logical,
intent(inout) :: flag(ixi^s,1:nw)
1543 {
do ix^db=ixomin^db,ixomax^db\}
1548 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1552 end subroutine mhd_check_w_inte
1554 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1557 logical,
intent(in) :: primitive
1558 integer,
intent(in) :: ixi^
l, ixo^
l
1559 double precision,
intent(in) :: w(ixi^s,nw)
1560 logical,
intent(inout) :: flag(ixi^s,1:nw)
1565 {
do ix^db=ixomin^db,ixomax^db\}
1570 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1574 end subroutine mhd_check_w_hde
1577 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1579 integer,
intent(in) :: ixi^
l, ixo^
l
1580 double precision,
intent(inout) :: w(ixi^s, nw)
1581 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1585 {
do ix^db=ixomin^db,ixomax^db\}
1587 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1589 +(^
c&w(ix^
d,
b^
c_)**2+))
1594 end subroutine mhd_to_conserved_origin
1597 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1599 integer,
intent(in) :: ixi^
l, ixo^
l
1600 double precision,
intent(inout) :: w(ixi^s, nw)
1601 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1605 {
do ix^db=ixomin^db,ixomax^db\}
1610 end subroutine mhd_to_conserved_origin_noe
1613 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1615 integer,
intent(in) :: ixi^
l, ixo^
l
1616 double precision,
intent(inout) :: w(ixi^s, nw)
1617 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1621 {
do ix^db=ixomin^db,ixomax^db\}
1623 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1629 end subroutine mhd_to_conserved_hde
1632 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1634 integer,
intent(in) :: ixi^
l, ixo^
l
1635 double precision,
intent(inout) :: w(ixi^s, nw)
1636 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1640 {
do ix^db=ixomin^db,ixomax^db\}
1642 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1647 end subroutine mhd_to_conserved_inte
1650 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1652 integer,
intent(in) :: ixi^
l, ixo^
l
1653 double precision,
intent(inout) :: w(ixi^s, nw)
1654 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1656 double precision :: rho
1659 {
do ix^db=ixomin^db,ixomax^db\}
1662 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1663 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1664 +(^
c&w(ix^
d,
b^
c_)**2+))
1669 end subroutine mhd_to_conserved_split_rho
1672 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1674 integer,
intent(in) :: ixi^
l, ixo^
l
1675 double precision,
intent(inout) :: w(ixi^s, nw)
1676 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1678 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1681 {
do ix^db=ixomin^db,ixomax^db\}
1683 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1684 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1685 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1686 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1687 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1688 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1693 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1694 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1695 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1703 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1707 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1709 +(^
c&w(ix^
d,
b^
c_)**2+)&
1710 +(^
c&e(ix^
d,^
c)**2+)*inv_squared_c)
1718 end subroutine mhd_to_conserved_semirelati
1720 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1722 integer,
intent(in) :: ixi^
l, ixo^
l
1723 double precision,
intent(inout) :: w(ixi^s, nw)
1724 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1726 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1729 {
do ix^db=ixomin^db,ixomax^db\}
1731 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1732 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1733 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1734 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1735 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1736 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1741 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1742 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1743 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1753 end subroutine mhd_to_conserved_semirelati_noe
1756 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1758 integer,
intent(in) :: ixi^
l, ixo^
l
1759 double precision,
intent(inout) :: w(ixi^s, nw)
1760 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1762 double precision :: inv_rho
1767 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1770 {
do ix^db=ixomin^db,ixomax^db\}
1771 inv_rho = 1.d0/w(ix^
d,
rho_)
1775 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1777 +(^
c&w(ix^
d,
b^
c_)**2+)))
1780 end subroutine mhd_to_primitive_origin
1783 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1785 integer,
intent(in) :: ixi^
l, ixo^
l
1786 double precision,
intent(inout) :: w(ixi^s, nw)
1787 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1789 double precision :: inv_rho
1794 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1797 {
do ix^db=ixomin^db,ixomax^db\}
1798 inv_rho = 1.d0/w(ix^
d,
rho_)
1803 end subroutine mhd_to_primitive_origin_noe
1806 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1808 integer,
intent(in) :: ixi^
l, ixo^
l
1809 double precision,
intent(inout) :: w(ixi^s, nw)
1810 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1812 double precision :: inv_rho
1817 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1820 {
do ix^db=ixomin^db,ixomax^db\}
1821 inv_rho = 1d0/w(ix^
d,
rho_)
1828 end subroutine mhd_to_primitive_hde
1831 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1833 integer,
intent(in) :: ixi^
l, ixo^
l
1834 double precision,
intent(inout) :: w(ixi^s, nw)
1835 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1837 double precision :: inv_rho
1842 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1845 {
do ix^db=ixomin^db,ixomax^db\}
1847 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1849 inv_rho = 1.d0/w(ix^
d,
rho_)
1853 end subroutine mhd_to_primitive_inte
1856 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1858 integer,
intent(in) :: ixi^
l, ixo^
l
1859 double precision,
intent(inout) :: w(ixi^s, nw)
1860 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1862 double precision :: inv_rho
1867 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1870 {
do ix^db=ixomin^db,ixomax^db\}
1875 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1877 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1880 end subroutine mhd_to_primitive_split_rho
1883 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1885 integer,
intent(in) :: ixi^
l, ixo^
l
1886 double precision,
intent(inout) :: w(ixi^s, nw)
1887 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1889 double precision ::
b(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
1894 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1897 {
do ix^db=ixomin^db,ixomax^db\}
1898 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1899 if(b2>smalldouble)
then
1907 inv_rho=1.d0/w(ix^
d,
rho_)
1909 b2=b2*inv_rho*inv_squared_c
1911 gamma2=1.d0/(1.d0+b2)
1913 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1917 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1921 b(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1922 b(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1923 b(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1927 b(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1933 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1935 +(^
c&w(ix^
d,
b^
c_)**2+)&
1936 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
1940 end subroutine mhd_to_primitive_semirelati
1943 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1945 integer,
intent(in) :: ixi^
l, ixo^
l
1946 double precision,
intent(inout) :: w(ixi^s, nw)
1947 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1949 double precision ::
b(ixo^s,1:
ndir),tmp,b2,gamma2,inv_rho
1950 integer :: ix^
d, idir
1954 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1957 {
do ix^db=ixomin^db,ixomax^db\}
1958 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1959 if(b2>smalldouble)
then
1967 inv_rho=1.d0/w(ix^
d,
rho_)
1969 b2=b2*inv_rho*inv_squared_c
1971 gamma2=1.d0/(1.d0+b2)
1973 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1976 end subroutine mhd_to_primitive_semirelati_noe
1981 integer,
intent(in) :: ixi^
l, ixo^
l
1982 double precision,
intent(inout) :: w(ixi^s, nw)
1983 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1988 {
do ix^db=ixomin^db,ixomax^db\}
1991 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1993 +(^
c&w(ix^
d,
b^
c_)**2+))
1996 {
do ix^db=ixomin^db,ixomax^db\}
1998 w(ix^d,
e_)=w(ix^d,
e_)&
1999 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2000 +(^
c&w(ix^d,
b^
c_)**2+))
2007 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
2009 integer,
intent(in) :: ixi^
l, ixo^
l
2010 double precision,
intent(inout) :: w(ixi^s, nw)
2011 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2015 {
do ix^db=ixomin^db,ixomax^db\}
2021 end subroutine mhd_ei_to_e_hde
2024 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
2026 integer,
intent(in) :: ixi^
l, ixo^
l
2027 double precision,
intent(inout) :: w(ixi^s, nw)
2028 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2030 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2031 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2033 end subroutine mhd_ei_to_e_semirelati
2038 integer,
intent(in) :: ixi^
l, ixo^
l
2039 double precision,
intent(inout) :: w(ixi^s, nw)
2040 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2045 {
do ix^db=ixomin^db,ixomax^db\}
2048 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2050 +(^
c&w(ix^
d,
b^
c_)**2+))
2053 {
do ix^db=ixomin^db,ixomax^db\}
2055 w(ix^d,
e_)=w(ix^d,
e_)&
2056 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2057 +(^
c&w(ix^d,
b^
c_)**2+))
2061 if(fix_small_values)
then
2062 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2068 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2070 integer,
intent(in) :: ixi^
l, ixo^
l
2071 double precision,
intent(inout) :: w(ixi^s, nw)
2072 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2076 {
do ix^db=ixomin^db,ixomax^db\}
2082 if(fix_small_values)
then
2083 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2086 end subroutine mhd_e_to_ei_hde
2089 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2091 integer,
intent(in) :: ixi^
l, ixo^
l
2092 double precision,
intent(inout) :: w(ixi^s, nw)
2093 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2095 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2096 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2098 end subroutine mhd_e_to_ei_semirelati
2100 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2103 logical,
intent(in) :: primitive
2104 integer,
intent(in) :: ixi^
l,ixo^
l
2105 double precision,
intent(inout) :: w(ixi^s,1:nw)
2106 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2107 character(len=*),
intent(in) :: subname
2109 double precision ::
b(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2110 double precision :: tmp, b2, gamma2, inv_rho
2112 logical :: flag(ixi^s,1:nw)
2121 {
do ix^db=ixomin^db,ixomax^db\}
2122 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2123 if(b2>smalldouble)
then
2130 inv_rho=1.d0/w(ix^
d,
rho_)
2132 b2=b2*inv_rho*inv_squared_c
2134 gamma2=1.d0/(1.d0+b2)
2136 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
2139 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2140 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2141 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2145 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2151 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2152 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2153 +(^
c&w(ix^
d,
b^
c_)**2+)&
2154 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
2161 select case (small_values_method)
2163 {
do ix^db=ixomin^db,ixomax^db\}
2164 if(flag(ix^d,
rho_))
then
2165 w(ix^d,
rho_) = small_density
2166 ^
c&w(ix^d,
m^
c_)=0.d0\
2170 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2172 if(flag(ix^d,
e_))
then
2173 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2174 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2181 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2184 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2186 w(ixo^s,
e_)=pressure(ixo^s)
2187 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2188 {
do ix^db=ixomin^db,ixomax^db\}
2189 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2190 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2195 if(.not.primitive)
then
2197 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2198 w(ixo^s,
e_)=pressure(ixo^s)
2200 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2204 end subroutine mhd_handle_small_values_semirelati
2206 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2209 logical,
intent(in) :: primitive
2210 integer,
intent(in) :: ixi^
l,ixo^
l
2211 double precision,
intent(inout) :: w(ixi^s,1:nw)
2212 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2213 character(len=*),
intent(in) :: subname
2216 logical :: flag(ixi^s,1:nw)
2218 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2223 {
do ix^db=ixomin^db,ixomax^db\}
2227 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2239 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2241 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2244 {
do ix^db=iximin^db,iximax^db\}
2245 w(ix^d,
e_)=w(ix^d,
e_)&
2246 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2248 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2250 {
do ix^db=iximin^db,iximax^db\}
2251 w(ix^d,
e_)=w(ix^d,
e_)&
2252 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2256 if(.not.primitive)
then
2258 {
do ix^db=ixomin^db,ixomax^db\}
2260 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2261 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2264 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2268 end subroutine mhd_handle_small_values_origin
2270 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2273 logical,
intent(in) :: primitive
2274 integer,
intent(in) :: ixi^
l,ixo^
l
2275 double precision,
intent(inout) :: w(ixi^s,1:nw)
2276 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2277 character(len=*),
intent(in) :: subname
2279 double precision :: rho
2281 logical :: flag(ixi^s,1:nw)
2283 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2288 {
do ix^db=ixomin^db,ixomax^db\}
2293 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2300 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2306 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2308 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2311 {
do ix^db=iximin^db,iximax^db\}
2313 w(ix^d,
e_)=w(ix^d,
e_)&
2314 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2316 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2318 {
do ix^db=iximin^db,iximax^db\}
2320 w(ix^d,
e_)=w(ix^d,
e_)&
2321 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2325 if(.not.primitive)
then
2327 {
do ix^db=ixomin^db,ixomax^db\}
2329 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2330 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2331 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2334 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2338 end subroutine mhd_handle_small_values_split
2340 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2343 logical,
intent(in) :: primitive
2344 integer,
intent(in) :: ixi^
l,ixo^
l
2345 double precision,
intent(inout) :: w(ixi^s,1:nw)
2346 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2347 character(len=*),
intent(in) :: subname
2350 logical :: flag(ixi^s,1:nw)
2352 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2357 {
do ix^db=ixomin^db,ixomax^db\}
2358 if(flag(ix^
d,
rho_))
then
2360 ^
c&w(ix^
d,
m^
c_)=0.d0\
2365 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2370 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2372 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2374 if(.not.primitive)
then
2376 {
do ix^db=ixomin^db,ixomax^db\}
2378 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2381 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2385 end subroutine mhd_handle_small_values_inte
2387 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2390 logical,
intent(in) :: primitive
2391 integer,
intent(in) :: ixi^
l,ixo^
l
2392 double precision,
intent(inout) :: w(ixi^s,1:nw)
2393 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2394 character(len=*),
intent(in) :: subname
2397 logical :: flag(ixi^s,1:nw)
2399 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2404 {
do ix^db=ixomin^db,ixomax^db\}
2408 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2414 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2416 if(.not.primitive)
then
2418 {
do ix^db=ixomin^db,ixomax^db\}
2422 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2426 end subroutine mhd_handle_small_values_noe
2428 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2431 logical,
intent(in) :: primitive
2432 integer,
intent(in) :: ixi^
l,ixo^
l
2433 double precision,
intent(inout) :: w(ixi^s,1:nw)
2434 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2435 character(len=*),
intent(in) :: subname
2438 logical :: flag(ixi^s,1:nw)
2440 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2445 {
do ix^db=ixomin^db,ixomax^db\}
2446 if(flag(ix^
d,
rho_))
then
2448 ^
c&w(ix^
d,
m^
c_)=0.d0\
2453 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2458 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2460 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2462 if(.not.primitive)
then
2464 {
do ix^db=ixomin^db,ixomax^db\}
2466 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2469 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2473 end subroutine mhd_handle_small_values_hde
2479 integer,
intent(in) :: ixi^
l, ixo^
l
2480 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2481 double precision,
intent(out) :: v(ixi^s,
ndir)
2483 double precision :: rho(ixi^s)
2488 rho(ixo^s)=1.d0/rho(ixo^s)
2491 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2497 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2500 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2501 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2502 double precision,
intent(inout) :: cmax(ixi^s)
2504 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2510 {
do ix^db=ixomin^db,ixomax^db \}
2525 cfast2=b2*inv_rho+cmax(ix^
d)
2526 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2527 if(avmincs2<zero) avmincs2=zero
2528 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2532 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2534 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2537 {
do ix^db=ixomin^db,ixomax^db \}
2551 b2=(^
c&w(ix^d,
b^
c_)**2+)
2552 cfast2=b2*inv_rho+cmax(ix^d)
2553 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2554 if(avmincs2<zero) avmincs2=zero
2555 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2559 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2561 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2565 end subroutine mhd_get_cmax_origin
2568 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2572 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2573 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2574 double precision,
intent(inout) :: cmax(ixi^s)
2576 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2577 double precision :: adiabs(ixo^s), gammas(ixo^s)
2593 {
do ix^db=ixomin^db,ixomax^db \}
2601 cmax(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*rho**(gammas(ix^
d)-1.d0)
2604 cfast2=b2*inv_rho+cmax(ix^
d)
2605 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2606 if(avmincs2<zero) avmincs2=zero
2607 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2611 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2613 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2616 {
do ix^db=ixomin^db,ixomax^db \}
2624 cmax(ix^d)=gammas(ix^d)*adiabs(ix^d)*rho**(gammas(ix^d)-1.d0)
2626 b2=(^
c&w(ix^d,
b^
c_)**2+)
2627 cfast2=b2*inv_rho+cmax(ix^d)
2628 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2629 if(avmincs2<zero) avmincs2=zero
2630 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2634 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2636 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2640 end subroutine mhd_get_cmax_origin_noe
2643 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2646 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2647 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2648 double precision,
intent(inout):: cmax(ixi^s)
2650 double precision :: csound, avmincs2, idim_alfven_speed2
2651 double precision :: inv_rho, alfven_speed2, gamma2
2654 {
do ix^db=ixomin^db,ixomax^db \}
2655 inv_rho=1.d0/w(ix^
d,
rho_)
2656 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2657 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2658 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2661 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2664 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2665 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2666 if(avmincs2<zero) avmincs2=zero
2668 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2669 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2672 end subroutine mhd_get_cmax_semirelati
2675 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2679 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2680 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2681 double precision,
intent(inout):: cmax(ixi^s)
2683 double precision :: adiabs(ixo^s), gammas(ixo^s)
2684 double precision :: csound, avmincs2, idim_alfven_speed2
2685 double precision :: inv_rho, alfven_speed2, gamma2
2699 {
do ix^db=ixomin^db,ixomax^db \}
2700 inv_rho=1.d0/w(ix^
d,
rho_)
2701 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2702 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2703 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2704 csound=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
2705 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2708 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2709 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2710 if(avmincs2<zero) avmincs2=zero
2712 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2713 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2716 end subroutine mhd_get_cmax_semirelati_noe
2718 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2721 integer,
intent(in) :: ixi^
l, ixo^
l
2722 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2723 double precision,
intent(inout) :: a2max(
ndim)
2724 double precision :: a2(ixi^s,
ndim,nw)
2725 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2730 hxo^
l=ixo^
l-
kr(i,^
d);
2731 gxo^
l=hxo^
l-
kr(i,^
d);
2732 jxo^
l=ixo^
l+
kr(i,^
d);
2733 kxo^
l=jxo^
l+
kr(i,^
d);
2734 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2735 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2736 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2738 end subroutine mhd_get_a2max
2741 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2744 integer,
intent(in) :: ixi^
l,ixo^
l
2745 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2747 double precision,
intent(inout) :: w(ixi^s,1:nw)
2748 double precision,
intent(out) :: tco_local,tmax_local
2750 double precision,
parameter :: trac_delta=0.25d0
2751 double precision :: te(ixi^s),lts(ixi^s)
2752 double precision,
dimension(1:ndim) :: bdir, bunitvec
2753 double precision,
dimension(ixI^S,1:ndim) :: gradt
2754 double precision :: ltrc,ltrp,altr
2755 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2756 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2759 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2762 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2765 tmax_local=maxval(te(ixo^s))
2773 do ix1=ixomin1,ixomax1
2774 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2775 if(lts(ix1)>trac_delta)
then
2776 tco_local=max(tco_local,te(ix1))
2788 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2789 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2790 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2791 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2793 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2804 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2811 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2816 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2817 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2821 if(bdir(1)/=0.d0)
then
2822 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2824 block%special_values(3)=0.d0
2826 if(bdir(2)/=0.d0)
then
2827 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2829 block%special_values(4)=0.d0
2833 if(bdir(1)/=0.d0)
then
2834 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2835 (bdir(3)/bdir(1))**2)
2837 block%special_values(3)=0.d0
2839 if(bdir(2)/=0.d0)
then
2840 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2841 (bdir(3)/bdir(2))**2)
2843 block%special_values(4)=0.d0
2845 if(bdir(3)/=0.d0)
then
2846 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2847 (bdir(2)/bdir(3))**2)
2849 block%special_values(5)=0.d0
2854 block%special_values(1)=zero
2855 {
do ix^db=ixomin^db,ixomax^db\}
2857 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2859 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2862 if(bdir(1)/=0.d0)
then
2863 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2867 if(bdir(2)/=0.d0)
then
2868 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2873 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2874 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2877 if(bdir(1)/=0.d0)
then
2878 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2882 if(bdir(2)/=0.d0)
then
2883 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2887 if(bdir(3)/=0.d0)
then
2888 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2893 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2894 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2896 if(lts(ix^d)>trac_delta)
then
2897 block%special_values(1)=max(block%special_values(1),te(ix^d))
2900 block%special_values(2)=tmax_local
2919 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2920 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2921 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2925 {
do ix^db=ixpmin^db,ixpmax^db\}
2926 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2928 if(bdir(1)/=0.d0)
then
2929 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2933 if(bdir(2)/=0.d0)
then
2934 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2940 if(bdir(1)/=0.d0)
then
2941 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2945 if(bdir(2)/=0.d0)
then
2946 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2950 if(bdir(3)/=0.d0)
then
2951 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2957 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2959 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2960 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2963 {
do ix^db=ixpmin^db,ixpmax^db\}
2965 if(w(ix^d,iw_mag(1))/=0.d0)
then
2966 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)
2970 if(w(ix^d,iw_mag(2))/=0.d0)
then
2971 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)
2977 if(w(ix^d,iw_mag(1))/=0.d0)
then
2978 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+&
2979 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
2983 if(w(ix^d,iw_mag(2))/=0.d0)
then
2984 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+&
2985 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
2989 if(w(ix^d,iw_mag(3))/=0.d0)
then
2990 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+&
2991 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
2997 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2999 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
3000 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
3006 {
do ix^db=ixpmin^db,ixpmax^db\}
3008 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
3009 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
3010 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
3013 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
3014 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
3015 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
3016 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
3022 call mpistop(
"unknown mhd_trac_type")
3025 end subroutine mhd_get_tcutoff
3028 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
3031 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3032 double precision,
intent(in) :: wprim(ixi^s, nw)
3033 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3034 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
3036 double precision :: csound(ixi^s,
ndim)
3037 double precision,
allocatable :: tmp(:^
d&)
3038 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
3042 allocate(tmp(ixa^s))
3044 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
3045 csound(ixa^s,id)=tmp(ixa^s)
3048 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
3049 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
3050 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
3051 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))
3055 ixamax^
d=ixcmax^
d+
kr(id,^
d);
3056 ixamin^
d=ixcmin^
d+
kr(id,^
d);
3057 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)))
3058 ixamax^
d=ixcmax^
d-
kr(id,^
d);
3059 ixamin^
d=ixcmin^
d-
kr(id,^
d);
3060 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)))
3065 ixamax^
d=jxcmax^
d+
kr(id,^
d);
3066 ixamin^
d=jxcmin^
d+
kr(id,^
d);
3067 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)))
3068 ixamax^
d=jxcmax^
d-
kr(id,^
d);
3069 ixamin^
d=jxcmin^
d-
kr(id,^
d);
3070 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)))
3074 end subroutine mhd_get_h_speed
3077 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3080 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3081 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3082 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3083 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3084 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3085 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3086 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3088 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3089 double precision :: umean, dmean, tmp1, tmp2, tmp3
3096 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3097 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3098 if(
present(cmin))
then
3099 {
do ix^db=ixomin^db,ixomax^db\}
3100 tmp1=sqrt(wlp(ix^
d,
rho_))
3101 tmp2=sqrt(wrp(ix^
d,
rho_))
3102 tmp3=1.d0/(tmp1+tmp2)
3103 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3104 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3105 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3106 cmin(ix^
d,1)=umean-dmean
3107 cmax(ix^
d,1)=umean+dmean
3109 if(h_correction)
then
3110 {
do ix^db=ixomin^db,ixomax^db\}
3111 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3112 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3116 {
do ix^db=ixomin^db,ixomax^db\}
3117 tmp1=sqrt(wlp(ix^d,
rho_))
3118 tmp2=sqrt(wrp(ix^d,
rho_))
3119 tmp3=1.d0/(tmp1+tmp2)
3120 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3121 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3122 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3123 cmax(ix^d,1)=abs(umean)+dmean
3127 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3128 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3129 if(
present(cmin))
then
3130 {
do ix^db=ixomin^db,ixomax^db\}
3131 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3132 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3134 if(h_correction)
then
3135 {
do ix^db=ixomin^db,ixomax^db\}
3136 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3137 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3141 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3145 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3146 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3147 if(
present(cmin))
then
3148 {
do ix^db=ixomin^db,ixomax^db\}
3149 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3150 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3151 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3153 if(h_correction)
then
3154 {
do ix^db=ixomin^db,ixomax^db\}
3155 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3156 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3160 {
do ix^db=ixomin^db,ixomax^db\}
3161 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3162 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3167 end subroutine mhd_get_cbounds
3170 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3173 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3174 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3175 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3176 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3177 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3178 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3179 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3181 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3186 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3187 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3189 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3190 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
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(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3196 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3199 {
do ix^db=ixomin^db,ixomax^db\}
3200 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3201 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3205 end subroutine mhd_get_cbounds_semirelati
3208 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3211 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3212 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3213 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3214 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3215 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3216 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3217 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3219 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3220 double precision :: umean, dmean, tmp1, tmp2, tmp3
3227 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3228 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3229 if(
present(cmin))
then
3230 {
do ix^db=ixomin^db,ixomax^db\}
3233 tmp3=1.d0/(tmp1+tmp2)
3234 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3235 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3236 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3237 cmin(ix^
d,1)=umean-dmean
3238 cmax(ix^
d,1)=umean+dmean
3240 if(h_correction)
then
3241 {
do ix^db=ixomin^db,ixomax^db\}
3242 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3243 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3247 {
do ix^db=ixomin^db,ixomax^db\}
3250 tmp3=1.d0/(tmp1+tmp2)
3251 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3252 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3253 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3254 cmax(ix^d,1)=abs(umean)+dmean
3258 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3259 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3260 if(
present(cmin))
then
3261 {
do ix^db=ixomin^db,ixomax^db\}
3262 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3263 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3265 if(h_correction)
then
3266 {
do ix^db=ixomin^db,ixomax^db\}
3267 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3268 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3272 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3276 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3277 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3278 if(
present(cmin))
then
3279 {
do ix^db=ixomin^db,ixomax^db\}
3280 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3281 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3282 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3284 if(h_correction)
then
3285 {
do ix^db=ixomin^db,ixomax^db\}
3286 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3287 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3291 {
do ix^db=ixomin^db,ixomax^db\}
3292 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3293 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3298 end subroutine mhd_get_cbounds_split_rho
3301 subroutine mhd_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3304 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3305 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3306 double precision,
intent(in) :: cmax(ixi^s)
3307 double precision,
intent(in),
optional :: cmin(ixi^s)
3308 type(ct_velocity),
intent(inout):: vcts
3310 end subroutine mhd_get_ct_velocity_average
3312 subroutine mhd_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3315 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3316 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3317 double precision,
intent(in) :: cmax(ixi^s)
3318 double precision,
intent(in),
optional :: cmin(ixi^s)
3319 type(ct_velocity),
intent(inout):: vcts
3321 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3323 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3325 end subroutine mhd_get_ct_velocity_contact
3327 subroutine mhd_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3330 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3331 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3332 double precision,
intent(in) :: cmax(ixi^s)
3333 double precision,
intent(in),
optional :: cmin(ixi^s)
3334 type(ct_velocity),
intent(inout):: vcts
3336 integer :: idime,idimn
3338 if(.not.
allocated(vcts%vbarC))
then
3339 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3340 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3343 if(
present(cmin))
then
3344 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3345 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3347 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3348 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3351 idimn=mod(idim,
ndir)+1
3352 idime=mod(idim+1,
ndir)+1
3354 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3355 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3356 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3357 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3358 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3360 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3361 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3362 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3363 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3364 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3366 end subroutine mhd_get_ct_velocity_hll
3369 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3373 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3374 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3375 double precision,
intent(out):: csound(ixo^s)
3377 double precision :: adiabs(ixo^s), gammas(ixo^s)
3378 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3398 {
do ix^db=ixomin^db,ixomax^db \}
3399 inv_rho=1.d0/w(ix^
d,
rho_)
3403 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3406 cfast2=b2*inv_rho+csound(ix^
d)
3407 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3409 if(avmincs2<zero) avmincs2=zero
3410 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3412 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3416 {
do ix^db=ixomin^db,ixomax^db \}
3417 inv_rho=1.d0/w(ix^d,
rho_)
3421 csound(ix^d)=gammas(ix^d)*adiabs(ix^d)*w(ix^d,
rho_)**(gammas(ix^d)-1.d0)
3423 b2=(^
c&w(ix^d,
b^
c_)**2+)
3424 cfast2=b2*inv_rho+csound(ix^d)
3425 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3426 if(avmincs2<zero) avmincs2=zero
3427 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3429 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3434 end subroutine mhd_get_csound_prim
3437 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3440 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3441 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3442 double precision,
intent(out):: csound(ixo^s)
3444 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3451 {
do ix^db=ixomin^db,ixomax^db \}
3458 cfast2=b2*inv_rho+csound(ix^
d)
3459 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3461 if(avmincs2<zero) avmincs2=zero
3462 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3464 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3468 {
do ix^db=ixomin^db,ixomax^db \}
3474 b2=(^
c&w(ix^d,
b^
c_)**2+)
3475 cfast2=b2*inv_rho+csound(ix^d)
3476 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3477 if(avmincs2<zero) avmincs2=zero
3478 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3480 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3485 end subroutine mhd_get_csound_prim_split
3488 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3491 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3493 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3494 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3496 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3499 {
do ix^db=ixomin^db,ixomax^db\}
3500 inv_rho = 1.d0/w(ix^
d,
rho_)
3503 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3504 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3505 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3506 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3509 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3510 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3511 if(avmincs2<zero) avmincs2=zero
3513 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3516 end subroutine mhd_get_csound_semirelati
3519 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3523 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3525 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3526 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3528 double precision :: adiabs(ixo^s), gammas(ixo^s)
3529 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3542 {
do ix^db=ixomin^db,ixomax^db\}
3543 inv_rho = 1.d0/w(ix^
d,
rho_)
3545 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3546 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3547 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3548 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3549 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3552 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3553 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3554 if(avmincs2<zero) avmincs2=zero
3556 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3559 end subroutine mhd_get_csound_semirelati_noe
3562 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3566 integer,
intent(in) :: ixi^
l, ixo^
l
3567 double precision,
intent(in) :: w(ixi^s,nw)
3568 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3569 double precision,
intent(out):: pth(ixi^s)
3571 double precision :: adiabs(ixo^s), gammas(ixo^s)
3584 {
do ix^db=ixomin^db,ixomax^db\}
3588 pth(ix^
d)=adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)
3592 end subroutine mhd_get_pthermal_noe
3595 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3599 integer,
intent(in) :: ixi^
l, ixo^
l
3600 double precision,
intent(in) :: w(ixi^s,nw)
3601 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3602 double precision,
intent(out):: pth(ixi^s)
3606 {
do ix^db= ixomin^db,ixomax^db\}
3610 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3615 if(check_small_values.and..not.fix_small_values)
then
3616 {
do ix^db= ixomin^db,ixomax^db\}
3617 if(pth(ix^d)<small_pressure)
then
3618 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3619 " encountered when call mhd_get_pthermal_inte"
3620 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3621 write(*,*)
"Location: ", x(ix^d,:)
3622 write(*,*)
"Cell number: ", ix^d
3624 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3627 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3628 write(*,*)
"Saving status at the previous time step"
3634 end subroutine mhd_get_pthermal_inte
3637 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3641 integer,
intent(in) :: ixi^
l, ixo^
l
3642 double precision,
intent(in) :: w(ixi^s,nw)
3643 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3644 double precision,
intent(out):: pth(ixi^s)
3648 {
do ix^db=ixomin^db,ixomax^db\}
3651 +(^
c&w(ix^
d,
b^
c_)**2+)))
3653 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3654 +(^
c&w(ix^
d,
b^
c_)**2+)))
3662 if(check_small_values.and..not.fix_small_values)
then
3663 {
do ix^db=ixomin^db,ixomax^db\}
3664 if(pth(ix^d)<small_pressure)
then
3665 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3666 " encountered when call mhd_get_pthermal"
3667 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3668 write(*,*)
"Location: ", x(ix^d,:)
3669 write(*,*)
"Cell number: ", ix^d
3671 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3674 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3675 write(*,*)
"Saving status at the previous time step"
3681 end subroutine mhd_get_pthermal_origin
3684 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3688 integer,
intent(in) :: ixi^
l, ixo^
l
3689 double precision,
intent(in) :: w(ixi^s,nw)
3690 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3691 double precision,
intent(out):: pth(ixi^s)
3693 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3696 {
do ix^db=ixomin^db,ixomax^db\}
3697 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3698 if(b2>smalldouble)
then
3706 inv_rho=1.d0/w(ix^
d,
rho_)
3708 b2=b2*inv_rho*inv_squared_c
3710 gamma2=1.d0/(1.d0+b2)
3712 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3716 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3717 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3718 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3722 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3728 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3729 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3730 +(^
c&w(ix^
d,
b^
c_)**2+)&
3731 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3735 if(check_small_values.and..not.fix_small_values)
then
3736 {
do ix^db=ixomin^db,ixomax^db\}
3737 if(pth(ix^d)<small_pressure)
then
3738 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3739 " encountered when call mhd_get_pthermal_semirelati"
3740 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3741 write(*,*)
"Location: ", x(ix^d,:)
3742 write(*,*)
"Cell number: ", ix^d
3744 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3747 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3748 write(*,*)
"Saving status at the previous time step"
3754 end subroutine mhd_get_pthermal_semirelati
3757 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3761 integer,
intent(in) :: ixi^
l, ixo^
l
3762 double precision,
intent(in) :: w(ixi^s,nw)
3763 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3764 double precision,
intent(out):: pth(ixi^s)
3768 {
do ix^db= ixomin^db,ixomax^db\}
3769 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3772 if(check_small_values.and..not.fix_small_values)
then
3773 {
do ix^db= ixomin^db,ixomax^db\}
3774 if(pth(ix^d)<small_pressure)
then
3775 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3776 " encountered when call mhd_get_pthermal_hde"
3777 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3778 write(*,*)
"Location: ", x(ix^d,:)
3779 write(*,*)
"Cell number: ", ix^d
3781 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3784 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3785 write(*,*)
"Saving status at the previous time step"
3791 end subroutine mhd_get_pthermal_hde
3794 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3796 integer,
intent(in) :: ixi^
l, ixo^
l
3797 double precision,
intent(in) :: w(ixi^s, 1:nw)
3798 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3799 double precision,
intent(out):: res(ixi^s)
3800 res(ixo^s) = w(ixo^s,
te_)
3801 end subroutine mhd_get_temperature_from_te
3804 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3806 integer,
intent(in) :: ixi^
l, ixo^
l
3807 double precision,
intent(in) :: w(ixi^s, 1:nw)
3808 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3809 double precision,
intent(out):: res(ixi^s)
3811 double precision :: r(ixi^s)
3814 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3815 end subroutine mhd_get_temperature_from_eint
3818 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3820 integer,
intent(in) :: ixi^
l, ixo^
l
3821 double precision,
intent(in) :: w(ixi^s, 1:nw)
3822 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3823 double precision,
intent(out):: res(ixi^s)
3825 double precision :: r(ixi^s),rho(ixi^s)
3830 res(ixo^s)=res(ixo^s)/(r(ixo^s)*rho(ixo^s))
3832 end subroutine mhd_get_temperature_from_etot
3834 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3836 integer,
intent(in) :: ixi^
l, ixo^
l
3837 double precision,
intent(in) :: w(ixi^s, 1:nw)
3838 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3839 double precision,
intent(out):: res(ixi^s)
3841 double precision :: r(ixi^s)
3847 end subroutine mhd_get_temperature_from_eint_with_equi
3849 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3851 integer,
intent(in) :: ixi^
l, ixo^
l
3852 double precision,
intent(in) :: w(ixi^s, 1:nw)
3853 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3854 double precision,
intent(out):: res(ixi^s)
3856 double precision :: r(ixi^s)
3862 end subroutine mhd_get_temperature_equi
3864 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3866 integer,
intent(in) :: ixi^
l, ixo^
l
3867 double precision,
intent(in) :: w(ixi^s, 1:nw)
3868 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3869 double precision,
intent(out):: res(ixi^s)
3871 end subroutine mhd_get_rho_equi
3873 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3875 integer,
intent(in) :: ixi^
l, ixo^
l
3876 double precision,
intent(in) :: w(ixi^s, 1:nw)
3877 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3878 double precision,
intent(out):: res(ixi^s)
3880 end subroutine mhd_get_pe_equi
3883 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3887 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3889 double precision,
intent(in) :: wc(ixi^s,nw)
3891 double precision,
intent(in) :: w(ixi^s,nw)
3892 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3893 double precision,
intent(out) :: f(ixi^s,nwflux)
3895 double precision :: vhall(ixi^s,1:
ndir)
3896 double precision :: ptotal
3900 {
do ix^db=ixomin^db,ixomax^db\}
3913 {
do ix^db=ixomin^db,ixomax^db\}
3917 ^
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_)\
3918 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3920 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3923 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3924 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3926 ^
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_)\
3930 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3931 {
do ix^db=ixomin^db,ixomax^db\}
3932 if(total_energy)
then
3934 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3935 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3938 ^
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))\
3942 {
do ix^db=ixomin^db,ixomax^db\}
3943 f(ix^d,mag(idim))=w(ix^d,
psi_)
3945 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3950 {
do ix^db=ixomin^db,ixomax^db\}
3956 {
do ix^db=ixomin^db,ixomax^db\}
3957 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(^d&w({ix^d},
b^d_)**2+)+smalldouble)
3962 end subroutine mhd_get_flux
3965 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3970 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3972 double precision,
intent(in) :: wc(ixi^s,nw)
3974 double precision,
intent(in) :: w(ixi^s,nw)
3975 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3976 double precision,
intent(out) :: f(ixi^s,nwflux)
3978 double precision :: vhall(ixi^s,1:
ndir)
3979 double precision :: adiabs(ixo^s), gammas(ixo^s)
3992 {
do ix^db=ixomin^db,ixomax^db\}
3998 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+)
4003 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4004 {
do ix^db=ixomin^db,ixomax^db\}
4006 ^
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))\
4010 {
do ix^db=ixomin^db,ixomax^db\}
4011 f(ix^d,mag(idim))=w(ix^d,
psi_)
4013 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4018 {
do ix^db=ixomin^db,ixomax^db\}
4023 end subroutine mhd_get_flux_noe
4026 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
4030 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4032 double precision,
intent(in) :: wc(ixi^s,nw)
4034 double precision,
intent(in) :: w(ixi^s,nw)
4035 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4036 double precision,
intent(out) :: f(ixi^s,nwflux)
4038 double precision :: vhall(ixi^s,1:
ndir)
4041 {
do ix^db=ixomin^db,ixomax^db\}
4054 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4055 {
do ix^db=ixomin^db,ixomax^db\}
4057 ^
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))\
4061 {
do ix^db=ixomin^db,ixomax^db\}
4062 f(ix^d,mag(idim))=w(ix^d,
psi_)
4064 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4069 {
do ix^db=ixomin^db,ixomax^db\}
4075 {
do ix^db=ixomin^db,ixomax^db\}
4076 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(^d&w({ix^d},
b^d_)**2+)+smalldouble)
4081 end subroutine mhd_get_flux_hde
4084 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
4088 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4090 double precision,
intent(in) :: wc(ixi^s,nw)
4092 double precision,
intent(in) :: w(ixi^s,nw)
4093 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4094 double precision,
intent(out) :: f(ixi^s,nwflux)
4096 double precision :: vhall(ixi^s,1:
ndir)
4097 double precision :: ptotal, btotal(ixo^s,1:
ndir)
4100 {
do ix^db=ixomin^db,ixomax^db\}
4108 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
4112 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
4116 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
4117 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4119 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
4123 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4126 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
4133 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
4134 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
4139 {
do ix^db=ixomin^db,ixomax^db\}
4140 f(ix^d,mag(idim))=w(ix^d,
psi_)
4142 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4147 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4148 {
do ix^db=ixomin^db,ixomax^db\}
4150 ^
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))\
4151 if(total_energy)
then
4153 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
4154 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
4160 {
do ix^db=ixomin^db,ixomax^db\}
4165 {
do ix^db=ixomin^db,ixomax^db\}
4166 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
4171 end subroutine mhd_get_flux_split
4174 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
4178 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4180 double precision,
intent(in) :: wc(ixi^s,nw)
4182 double precision,
intent(in) :: w(ixi^s,nw)
4183 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4184 double precision,
intent(out) :: f(ixi^s,nwflux)
4186 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
4189 {
do ix^db=ixomin^db,ixomax^db\}
4194 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4195 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4196 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4201 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4206 e2=(^
c&e(ix^
d,^
c)**2+)
4213 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4214 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4215 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4218 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4219 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4232 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4234 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)
4241 {
do ix^db=ixomin^db,ixomax^db\}
4242 f(ix^d,mag(idim))=w(ix^d,
psi_)
4244 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4249 {
do ix^db=ixomin^db,ixomax^db\}
4254 {
do ix^db=ixomin^db,ixomax^db\}
4255 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*w(ix^d,mag(idim))/(dsqrt(^d&w({ix^d},
b^d_)**2+)+smalldouble)
4260 end subroutine mhd_get_flux_semirelati
4262 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4267 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4269 double precision,
intent(in) :: wc(ixi^s,nw)
4271 double precision,
intent(in) :: w(ixi^s,nw)
4272 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4273 double precision,
intent(out) :: f(ixi^s,nwflux)
4275 double precision :: adiabs(ixo^s), gammas(ixo^s)
4276 double precision :: e(ixo^s,1:
ndir),e2
4289 {
do ix^db=ixomin^db,ixomax^db\}
4294 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4295 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4296 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4297 e2=(^
c&e(ix^
d,^
c)**2+)
4302 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4311 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4313 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)
4320 {
do ix^db=ixomin^db,ixomax^db\}
4321 f(ix^d,mag(idim))=w(ix^d,
psi_)
4323 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4328 {
do ix^db=ixomin^db,ixomax^db\}
4333 end subroutine mhd_get_flux_semirelati_noe
4339 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4341 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4342 double precision,
intent(in) :: qdt
4343 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4344 double precision,
intent(inout) :: w(ixi^s,1:nw)
4345 double precision :: tmp(ixi^s)
4346 double precision :: jxbxb(ixi^s,1:3)
4348 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4351 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4353 end subroutine add_source_ambipolar_internal_energy
4355 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4358 integer,
intent(in) :: ixi^
l, ixo^
l
4359 double precision,
intent(in) :: w(ixi^s,nw)
4360 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4361 double precision,
intent(out) :: res(:^
d&,:)
4363 double precision :: btot(ixi^s,1:3)
4364 double precision :: current(ixi^s,7-2*
ndir:3)
4365 double precision :: tmp(ixi^s),b2(ixi^s)
4366 integer :: idir, idirmin
4375 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4378 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4381 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4382 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4384 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4387 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4389 end subroutine mhd_get_jxbxb
4396 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4400 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4401 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4402 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4403 double precision,
intent(in) :: my_dt
4404 logical,
intent(in) :: fix_conserve_at_step
4406 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4407 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4408 double precision :: fe(ixi^s,
sdim:3)
4409 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4410 integer :: i, ixa^
l, ie_
4416 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4426 btot(ixa^s,1:3)=0.d0
4432 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4435 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4436 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4438 wres(ixo^s,
e_)=-tmp2(ixo^s)
4444 ff(ixa^s,1) = tmp(ixa^s,2)
4445 ff(ixa^s,2) = -tmp(ixa^s,1)
4447 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4448 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4449 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4452 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4454 ixamin^
d=ixomin^
d-1;
4455 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4464 ff(ixa^s,2) = tmp(ixa^s,3)
4465 ff(ixa^s,3) = -tmp(ixa^s,2)
4466 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4467 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4469 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4471 ff(ixa^s,1) = -tmp(ixa^s,3)
4473 ff(ixa^s,3) = tmp(ixa^s,1)
4474 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4475 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4476 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4480 ff(ixa^s,1) = tmp(ixa^s,2)
4481 ff(ixa^s,2) = -tmp(ixa^s,1)
4483 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4484 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4485 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4490 if(fix_conserve_at_step)
then
4491 fluxall=my_dt*fluxall
4498 end subroutine sts_set_source_ambipolar
4501 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4504 integer,
intent(in) :: ixi^
l, ixo^
l
4505 double precision,
intent(in) :: w(ixi^s,1:nw)
4506 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4508 double precision,
intent(in) :: ecc(ixi^s,1:3)
4509 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4510 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4512 integer :: hxc^
l,ixc^
l,ixa^
l
4513 integer :: idim1,idim2,idir,ix^
d
4519 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4521 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4522 ixamin^
d=ixcmin^
d+ix^
d;
4523 ixamax^
d=ixcmax^
d+ix^
d;
4524 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4526 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4532 ixcmin^d=ixomin^d-1;
4540 hxc^l=ixc^l-kr(idim2,^d);
4542 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4543 +lvc(idim1,idim2,idir)&
4548 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4551 end subroutine update_faces_ambipolar
4555 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4558 integer,
intent(in) :: ixi^
l, ixo^
l
4559 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4560 double precision,
intent(out) :: src(ixi^s)
4562 double precision :: ffc(ixi^s,1:
ndim)
4563 double precision :: dxinv(
ndim)
4564 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4570 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4572 ixbmin^
d=ixcmin^
d+ix^
d;
4573 ixbmax^
d=ixcmax^
d+ix^
d;
4576 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4578 ff(ixi^s,1:ndim)=0.d0
4580 ixb^l=ixo^l-kr(idims,^d);
4581 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4583 if({ ix^d==0 .and. ^d==idims | .or.})
then
4584 ixbmin^d=ixcmin^d-ix^d;
4585 ixbmax^d=ixcmax^d-ix^d;
4586 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4589 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4592 if(slab_uniform)
then
4594 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4595 ixb^l=ixo^l-kr(idims,^d);
4596 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4600 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4601 ixb^l=ixo^l-kr(idims,^d);
4602 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4604 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4606 end subroutine get_flux_on_cell_face
4610 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4613 integer,
intent(in) :: ixi^
l, ixo^
l
4614 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4615 double precision,
intent(in) :: w(ixi^s,1:nw)
4616 double precision :: dtnew
4618 double precision :: coef
4619 double precision :: dxarr(
ndim)
4620 double precision :: tmp(ixi^s)
4625 coef = maxval(abs(tmp(ixo^s)))
4632 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4634 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4637 end function get_ambipolar_dt
4645 integer,
intent(in) :: ixi^
l, ixo^
l
4646 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4647 double precision,
intent(inout) :: res(ixi^s)
4648 double precision :: tmp(ixi^s)
4649 double precision :: rho(ixi^s)
4658 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4662 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4669 integer,
intent(in) :: ixi^
l, ixo^
l
4670 double precision,
intent(in) :: qdt,dtfactor
4671 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4672 double precision,
intent(inout) :: w(ixi^s,1:nw)
4673 logical,
intent(in) :: qsourcesplit
4674 logical,
intent(inout) :: active
4681 if (.not. qsourcesplit)
then
4685 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4689 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4695 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4701 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4705 if (abs(
mhd_eta)>smalldouble)
then
4707 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4712 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4718 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4722 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4729 select case (type_divb)
4734 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4737 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4740 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4741 case (divb_janhunen)
4743 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4744 case (divb_lindejanhunen)
4746 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4747 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4748 case (divb_lindepowel)
4750 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4751 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4752 case (divb_lindeglm)
4754 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4755 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4756 case (divb_multigrid)
4761 call mpistop(
'Unknown divB fix')
4768 w,x,qsourcesplit,active,
rc_fl)
4778 w,x,gravity_energy,qsourcesplit,active)
4787 if(.not.qsourcesplit)
then
4789 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4793 end subroutine mhd_add_source
4795 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4799 integer,
intent(in) :: ixi^
l, ixo^
l
4800 double precision,
intent(in) :: qdt,dtfactor
4801 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4802 double precision,
intent(inout) :: w(ixi^s,1:nw)
4803 double precision :: divv(ixi^s)
4819 end subroutine add_pe0_divv
4821 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4823 integer,
intent(in) :: ixi^
l,ixo^
l
4824 double precision,
intent(in) :: qdt
4825 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4826 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4827 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4829 double precision :: r(ixi^s),te(ixi^s),rho_loc(ixi^s),pth_loc(ixi^s)
4830 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,bdir(
ndim),bunitvec(
ndim)
4834 {
do ix^db=iximin^db,iximax^db\}
4838 rho_loc(ix^
d)=wctprim(ix^
d,
rho_)
4843 pth_loc(ix^
d)=wctprim(ix^
d,
p_)
4845 te(ix^
d)=pth_loc(ix^
d)/(r(ix^
d)*rho_loc(ix^
d))
4851 do ix1=ixomin1,ixomax1
4853 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4855 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4858 sigma_t7=sigma_t5*te(ix^d)
4862 sigma_t7=sigma_t5*te(ix^d)
4864 sigmat5_bgradt=sigma_t5*(8.d0*(te(ix1+1)-te(ix1-1))-te(ix1+2)+te(ix1-2))/12.d0/block%ds(ix^d,1)
4866 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4867 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4868 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4870 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4871 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4876 do ix2=ixomin2,ixomax2
4877 do ix1=ixomin1,ixomax1
4879 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4881 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4884 sigma_t7=sigma_t5*te(ix^d)
4888 sigma_t7=sigma_t5*te(ix^d)
4891 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4893 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4895 if(bdir(1)/=0.d0)
then
4896 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
4900 if(bdir(2)/=0.d0)
then
4901 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
4905 sigmat5_bgradt=sigma_t5*(&
4906 bunitvec(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)&
4907 +bunitvec(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))
4909 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4910 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4911 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4913 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4914 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4920 do ix3=ixomin3,ixomax3
4921 do ix2=ixomin2,ixomax2
4922 do ix1=ixomin1,ixomax1
4924 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4926 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4929 sigma_t7=sigma_t5*te(ix^d)
4933 sigma_t7=sigma_t5*te(ix^d)
4936 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4938 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4940 if(bdir(1)/=0.d0)
then
4941 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
4945 if(bdir(2)/=0.d0)
then
4946 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
4950 if(bdir(3)/=0.d0)
then
4951 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
4955 sigmat5_bgradt=sigma_t5*(&
4956 bunitvec(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)&
4957 +bunitvec(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)&
4958 +bunitvec(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))
4960 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4961 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4962 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4964 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4965 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4971 end subroutine add_hypertc_source
4974 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4976 integer,
intent(in) :: ixi^
l, ixo^
l
4977 double precision,
intent(in) :: w(ixi^s,1:nw)
4978 double precision,
intent(inout) :: jxb(ixi^s,3)
4979 double precision :: a(ixi^s,3),
b(ixi^s,3)
4981 double precision :: current(ixi^s,7-2*
ndir:3)
4982 integer :: idir, idirmin
4987 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4991 b(ixo^s, idir) = w(ixo^s,mag(idir))
5000 a(ixo^s,idir)=current(ixo^s,idir)
5004 end subroutine get_lorentz_force
5008 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
5010 integer,
intent(in) :: ixi^
l, ixo^
l
5011 double precision,
intent(in) :: w(ixi^s, nw)
5012 double precision,
intent(out) :: gamma_a2(ixo^s)
5013 double precision :: rho(ixi^s)
5019 rho(ixo^s) = w(ixo^s,
rho_)
5022 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
5023 end subroutine mhd_gamma2_alfven
5027 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
5029 integer,
intent(in) :: ixi^
l, ixo^
l
5030 double precision,
intent(in) :: w(ixi^s, nw)
5031 double precision :: gamma_a(ixo^s)
5033 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
5034 gamma_a = sqrt(gamma_a)
5035 end function mhd_gamma_alfven
5039 integer,
intent(in) :: ixi^
l, ixo^
l
5040 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
5041 double precision,
intent(out) :: rho(ixi^s)
5046 rho(ixo^s) = w(ixo^s,
rho_)
5052 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
5055 integer,
intent(in) :: ixi^
l,ixo^
l, ie
5056 double precision,
intent(inout) :: w(ixi^s,1:nw)
5057 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5058 character(len=*),
intent(in) :: subname
5060 double precision :: rho(ixi^s)
5062 logical :: flag(ixi^s,1:nw)
5066 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
5067 flag(ixo^s,ie)=.true.
5069 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
5071 if(any(flag(ixo^s,ie)))
then
5075 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
5078 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
5084 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
5087 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
5093 end subroutine mhd_handle_small_ei
5095 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
5099 integer,
intent(in) :: ixi^
l, ixo^
l
5100 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5101 double precision,
intent(inout) :: w(ixi^s,1:nw)
5103 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
5112 end subroutine mhd_update_temperature
5115 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
5118 integer,
intent(in) :: ixi^
l, ixo^
l
5119 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5120 double precision,
intent(inout) :: w(ixi^s,1:nw)
5122 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
5134 a(ixo^s,idir)=
block%J0(ixo^s,idir)
5139 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5142 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5148 if(total_energy)
then
5151 b(ixo^s,:)=wct(ixo^s,mag(:))
5160 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5163 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5167 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5171 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5176 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5181 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5183 end subroutine add_source_b0split
5186 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5190 integer,
intent(in) :: ixi^
l, ixo^
l
5191 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5192 double precision,
intent(inout) :: w(ixi^s,1:nw)
5193 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5195 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5196 integer :: idir, idirmin, ix^
d
5200 {
do ix^db=iximin^db,iximax^db\}
5202 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5203 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5204 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5206 call divvector(e,ixi^l,ixo^l,dive)
5208 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5211 {
do ix^db=ixomin^db,ixomax^db\}
5212 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5213 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5214 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5215 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5216 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5217 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5221 end subroutine add_source_semirelativistic
5224 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5228 integer,
intent(in) :: ixi^
l, ixo^
l
5229 double precision,
intent(in) :: qdt
5230 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5231 double precision,
intent(inout) :: w(ixi^s,1:nw)
5232 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5234 double precision :: divv(ixi^s), tmp
5246 {
do ix^db=ixomin^db,ixomax^db\}
5248 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5249 if(w(ix^
d,
e_)<small_e)
then
5254 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
5257 if(fix_small_values)
then
5258 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5260 end subroutine add_source_internal_e
5263 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5268 integer,
intent(in) :: ixi^
l, ixo^
l
5269 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5270 double precision,
intent(inout) :: w(ixi^s,1:nw)
5271 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5273 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5274 double precision :: current(ixi^s,7-2*
ndir:3)
5275 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5276 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5277 integer :: idir, idirmin, idims, ix^
d
5282 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5290 j(ixo^s,idir)=current(ixo^s,idir)
5308 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
5314 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5321 b(ixo^s,idir)=wct(ixo^s,
rho_)*(
b(ixo^s,idir)-gravity_field(ixo^s,idir))+j(ixo^s,idir)-jxb(ixo^s,idir)
5325 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5329 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5330 tmp(ixo^s)=sqrt(b2(ixo^s))
5331 where(tmp(ixo^s)>smalldouble)
5332 tmp(ixo^s)=1.d0/tmp(ixo^s)
5338 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5343 {
do ix^db=ixomin^db,ixomax^db\}
5345 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5347 b2(ix^
d)=vaoc/(1.d0+vaoc)
5350 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5353 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5357 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5360 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5361 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5364 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5367 end subroutine add_source_hydrodynamic_e
5373 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5378 integer,
intent(in) :: ixi^
l, ixo^
l
5379 double precision,
intent(in) :: qdt
5380 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5381 double precision,
intent(inout) :: w(ixi^s,1:nw)
5382 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5383 integer :: lxo^
l, kxo^
l
5385 double precision :: tmp(ixi^s),tmp2(ixi^s)
5388 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5389 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5398 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5399 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5406 gradeta(ixo^s,1:
ndim)=zero
5412 gradeta(ixo^s,idim)=tmp(ixo^s)
5419 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5426 tmp2(ixi^s)=bf(ixi^s,idir)
5428 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5429 jxo^
l=ixo^
l+
kr(idim,^
d);
5430 hxo^
l=ixo^
l-
kr(idim,^
d);
5431 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5432 tmp(ixo^s)=tmp(ixo^s)+&
5433 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5438 tmp2(ixi^s)=bf(ixi^s,idir)
5440 jxo^
l=ixo^
l+
kr(idim,^
d);
5441 hxo^
l=ixo^
l-
kr(idim,^
d);
5442 tmp(ixo^s)=tmp(ixo^s)+&
5443 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5448 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5452 do jdir=1,
ndim;
do kdir=idirmin,3
5453 if (
lvc(idir,jdir,kdir)/=0)
then
5454 if (
lvc(idir,jdir,kdir)==1)
then
5455 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5457 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5464 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5465 if(total_energy)
then
5466 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5472 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5475 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5477 end subroutine add_source_res1
5481 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5486 integer,
intent(in) :: ixi^
l, ixo^
l
5487 double precision,
intent(in) :: qdt
5488 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5489 double precision,
intent(inout) :: w(ixi^s,1:nw)
5492 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5493 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5494 integer :: ixa^
l,idir,idirmin,idirmin1
5498 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5499 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5509 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5514 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5523 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5526 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5531 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5533 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5535 if(total_energy)
then
5538 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5539 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5542 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5546 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5547 end subroutine add_source_res2
5551 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5555 integer,
intent(in) :: ixi^
l, ixo^
l
5556 double precision,
intent(in) :: qdt
5557 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5558 double precision,
intent(inout) :: w(ixi^s,1:nw)
5560 double precision :: current(ixi^s,7-2*
ndir:3)
5561 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5562 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5565 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5566 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5569 tmpvec(ixa^s,1:
ndir)=zero
5571 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5575 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5578 tmpvec(ixa^s,1:
ndir)=zero
5579 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5583 tmpvec2(ixa^s,1:
ndir)=zero
5584 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5587 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5590 if(total_energy)
then
5593 tmpvec2(ixa^s,1:
ndir)=zero
5594 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5595 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5596 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5597 end do;
end do;
end do
5599 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5600 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5603 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5605 end subroutine add_source_hyperres
5607 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5614 integer,
intent(in) :: ixi^
l, ixo^
l
5615 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5616 double precision,
intent(inout) :: w(ixi^s,1:nw)
5618 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5639 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5642 if(total_energy)
then
5651 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5660 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5664 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5666 end subroutine add_source_glm
5669 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5672 integer,
intent(in) :: ixi^
l, ixo^
l
5673 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5674 double precision,
intent(inout) :: w(ixi^s,1:nw)
5676 double precision :: divb(ixi^s), ba(1:
ndir)
5677 integer :: idir, ix^
d
5683 {
do ix^db=ixomin^db,ixomax^db\}
5688 if (total_energy)
then
5694 {
do ix^db=ixomin^db,ixomax^db\}
5696 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5698 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5699 if (total_energy)
then
5701 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5706 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5708 end subroutine add_source_powel
5710 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5715 integer,
intent(in) :: ixi^
l, ixo^
l
5716 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5717 double precision,
intent(inout) :: w(ixi^s,1:nw)
5719 double precision :: divb(ixi^s)
5720 integer :: idir, ix^
d
5725 {
do ix^db=ixomin^db,ixomax^db\}
5730 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5732 end subroutine add_source_janhunen
5734 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5739 integer,
intent(in) :: ixi^
l, ixo^
l
5740 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5741 double precision,
intent(inout) :: w(ixi^s,1:nw)
5743 double precision :: divb(ixi^s),graddivb(ixi^s)
5744 integer :: idim, idir, ixp^
l, i^
d, iside
5745 logical,
dimension(-1:1^D&) :: leveljump
5753 if(i^
d==0|.and.) cycle
5754 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5755 leveljump(i^
d)=.true.
5757 leveljump(i^
d)=.false.
5766 i^dd=kr(^dd,^d)*(2*iside-3);
5767 if (leveljump(i^dd))
then
5769 ixpmin^d=ixomin^d-i^d
5771 ixpmax^d=ixomax^d-i^d
5782 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5784 {
do i^db=ixpmin^db,ixpmax^db\}
5786 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5788 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5790 if (typedivbdiff==
'all' .and. total_energy)
then
5792 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5797 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5799 end subroutine add_source_linde
5806 integer,
intent(in) :: ixi^
l, ixo^
l
5807 double precision,
intent(in) :: w(ixi^s,1:nw)
5808 double precision :: divb(ixi^s), dsurface(ixi^s)
5810 double precision :: invb(ixo^s)
5811 integer :: ixa^
l,idims
5813 call get_divb(w,ixi^
l,ixo^
l,divb)
5815 where(invb(ixo^s)/=0.d0)
5816 invb(ixo^s)=1.d0/invb(ixo^s)
5819 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5821 ixamin^
d=ixomin^
d-1;
5822 ixamax^
d=ixomax^
d-1;
5823 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5825 ixa^
l=ixo^
l-
kr(idims,^
d);
5826 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5828 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5829 block%dvolume(ixo^s)/dsurface(ixo^s)
5840 integer,
intent(in) :: ixo^
l, ixi^
l
5841 double precision,
intent(in) :: w(ixi^s,1:nw)
5842 integer,
intent(out) :: idirmin
5845 double precision :: current(ixi^s,7-2*
ndir:3)
5846 integer :: idir, idirmin0
5852 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5853 block%J0(ixo^s,idirmin0:3)
5857 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5865 integer,
intent(in) :: ixi^
l, ixo^
l
5866 double precision,
intent(inout) :: dtnew
5867 double precision,
intent(in) ::
dx^
d
5868 double precision,
intent(in) :: w(ixi^s,1:nw)
5869 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5871 double precision :: dxarr(
ndim)
5872 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5873 integer :: idirmin,idim
5887 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5890 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5916 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5923 end subroutine mhd_get_dt
5926 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5932 integer,
intent(in) :: ixi^
l, ixo^
l
5933 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5934 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5936 double precision :: adiabs(ixo^s), gammas(ixo^s)
5937 double precision :: tmp,tmp1,invr,cot
5939 integer :: mr_,mphi_
5940 integer :: br_,bphi_
5943 br_=mag(1); bphi_=mag(1)-1+
phi_
5960 {
do ix^db=ixomin^db,ixomax^db\}
5963 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5968 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5970 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5973 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5974 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5975 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5976 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5977 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5979 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5980 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5981 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5984 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5989 {
do ix^db=ixomin^db,ixomax^db\}
5991 if(local_timestep)
then
5992 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5997 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5999 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6003 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6006 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6007 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6011 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6017 cot=1.d0/tan(x(ix^d,2))
6021 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6022 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6024 if(.not.stagger_grid)
then
6025 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6027 tmp=tmp+wprim(ix^d,
psi_)*cot
6029 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6034 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6035 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6036 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6038 if(.not.stagger_grid)
then
6039 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6041 tmp=tmp+wprim(ix^d,
psi_)*cot
6043 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6046 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6047 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6048 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6049 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6050 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6052 if(.not.stagger_grid)
then
6053 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6054 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6055 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6056 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6057 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6064 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6067 end subroutine mhd_add_source_geom
6070 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6076 integer,
intent(in) :: ixi^
l, ixo^
l
6077 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6078 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6080 double precision :: adiabs(ixo^s), gammas(ixo^s)
6081 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
6083 integer :: mr_,mphi_
6084 integer :: br_,bphi_
6087 br_=mag(1); bphi_=mag(1)-1+
phi_
6104 {
do ix^db=ixomin^db,ixomax^db\}
6107 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6114 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)
6118 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
6119 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
6120 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6125 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6131 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
6132 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
6133 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
6134 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6135 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6136 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
6138 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6139 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6140 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6143 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
6144 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
6149 {
do ix^db=ixomin^db,ixomax^db\}
6151 if(local_timestep)
then
6152 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
6158 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
6159 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
6160 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6164 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6171 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6173 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6177 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
6180 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
6181 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
6182 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
6186 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
6192 cot=1.d0/tan(x(ix^d,2))
6196 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
6197 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
6199 if(.not.stagger_grid)
then
6200 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6202 tmp=tmp+wprim(ix^d,
psi_)*cot
6204 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6210 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6211 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
6212 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6213 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
6215 if(.not.stagger_grid)
then
6216 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6218 tmp=tmp+wprim(ix^d,
psi_)*cot
6220 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6223 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6224 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6225 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6226 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
6227 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6228 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6229 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
6231 if(.not.stagger_grid)
then
6232 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6233 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6234 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6235 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6236 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6243 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6246 end subroutine mhd_add_source_geom_semirelati
6249 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6254 integer,
intent(in) :: ixi^
l, ixo^
l
6255 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6256 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6258 double precision :: tmp,tmp1,tmp2,invr,cot
6260 integer :: mr_,mphi_
6261 integer :: br_,bphi_
6264 br_=mag(1); bphi_=mag(1)-1+
phi_
6269 {
do ix^db=ixomin^db,ixomax^db\}
6272 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6276 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6278 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6279 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6280 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6281 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6282 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6284 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6285 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6286 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6289 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6294 {
do ix^db=ixomin^db,ixomax^db\}
6296 if(local_timestep)
then
6297 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6301 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6302 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6305 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6309 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6310 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6311 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6313 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6314 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6319 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6325 cot=1.d0/tan(x(ix^d,2))
6330 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6331 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6332 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6334 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6335 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6338 if(.not.stagger_grid)
then
6340 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6341 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6343 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6346 tmp=tmp+wprim(ix^d,
psi_)*cot
6348 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6354 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6355 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6356 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6357 +(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)
6359 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6360 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6361 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6364 if(.not.stagger_grid)
then
6366 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6367 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6369 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6372 tmp=tmp+wprim(ix^d,
psi_)*cot
6374 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6378 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6379 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6380 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6381 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6382 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6383 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6384 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6385 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6386 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6388 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6389 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6390 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6391 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6392 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6395 if(.not.stagger_grid)
then
6397 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6398 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6399 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6400 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6401 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6402 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6403 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6404 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6405 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6407 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6408 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6409 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6410 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6411 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6419 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6422 end subroutine mhd_add_source_geom_split
6427 integer,
intent(in) :: ixi^
l, ixo^
l
6428 double precision,
intent(in) :: w(ixi^s, nw)
6429 double precision :: mge(ixo^s)
6432 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6434 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6438 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6441 integer,
intent(in) :: ixi^
l, ixo^
l
6442 double precision,
intent(in) :: w(ixi^s,nw)
6443 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6444 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6446 double precision :: current(ixi^s,7-2*
ndir:3)
6447 double precision :: rho(ixi^s)
6448 integer :: idir, idirmin, ix^
d
6453 do idir = idirmin,
ndir
6454 {
do ix^db=ixomin^db,ixomax^db\}
6455 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6459 end subroutine mhd_getv_hall
6461 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6464 integer,
intent(in) :: ixi^
l, ixo^
l
6465 double precision,
intent(in) :: w(ixi^s,nw)
6466 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6467 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6470 double precision :: current(ixi^s,7-2*
ndir:3)
6471 integer :: idir, idirmin
6478 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6479 do idir = idirmin, 3
6483 end subroutine mhd_get_jambi
6485 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6488 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6489 double precision,
intent(in) :: qt
6490 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6491 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6494 double precision :: db(ixo^s), dpsi(ixo^s)
6498 {
do ix^db=ixomin^db,ixomax^db\}
6499 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6500 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6501 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6502 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6511 {
do ix^db=ixomin^db,ixomax^db\}
6512 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6513 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6514 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6515 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6516 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6518 if(total_energy)
then
6519 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6520 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6522 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6524 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6527 if(total_energy)
then
6528 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6529 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6534 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6536 end subroutine mhd_modify_wlr
6538 subroutine mhd_boundary_adjust(igrid,psb)
6540 integer,
intent(in) :: igrid
6543 integer :: ib, idims, iside, ixo^
l, i^
d
6552 i^
d=
kr(^
d,idims)*(2*iside-3);
6553 if (neighbor_type(i^
d,igrid)/=1) cycle
6554 ib=(idims-1)*2+iside
6572 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6577 end subroutine mhd_boundary_adjust
6579 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6582 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6583 double precision,
intent(inout) :: w(ixg^s,1:nw)
6584 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6586 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6587 integer :: ix^
d,ixf^
l
6600 do ix1=ixfmax1,ixfmin1,-1
6601 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6602 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6603 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6606 do ix1=ixfmax1,ixfmin1,-1
6607 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6608 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6609 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6610 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6611 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6612 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6613 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6627 do ix1=ixfmax1,ixfmin1,-1
6628 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6629 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6630 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6631 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6632 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6633 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6636 do ix1=ixfmax1,ixfmin1,-1
6637 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6638 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6639 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6640 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6641 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6642 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6643 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6644 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6645 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6646 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6647 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6648 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6649 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6650 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6651 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6652 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6653 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6654 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6668 do ix1=ixfmin1,ixfmax1
6669 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6670 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6671 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6674 do ix1=ixfmin1,ixfmax1
6675 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6676 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6677 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6678 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6679 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6680 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6681 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6695 do ix1=ixfmin1,ixfmax1
6696 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6697 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6698 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6699 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6700 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6701 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6704 do ix1=ixfmin1,ixfmax1
6705 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6706 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6707 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6708 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6709 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6710 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6711 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6712 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6713 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6714 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6715 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6716 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6717 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6718 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6719 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6720 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6721 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6722 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6736 do ix2=ixfmax2,ixfmin2,-1
6737 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6738 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6739 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6742 do ix2=ixfmax2,ixfmin2,-1
6743 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6744 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6745 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6746 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6747 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6748 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6749 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6763 do ix2=ixfmax2,ixfmin2,-1
6764 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6765 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6766 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6767 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6768 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6769 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6772 do ix2=ixfmax2,ixfmin2,-1
6773 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6774 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6775 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6776 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6777 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6778 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6779 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6780 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6781 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6782 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6783 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6784 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6785 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6786 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6787 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6788 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6789 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6790 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6804 do ix2=ixfmin2,ixfmax2
6805 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6806 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6807 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6810 do ix2=ixfmin2,ixfmax2
6811 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6812 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6813 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6814 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6815 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6816 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6817 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6831 do ix2=ixfmin2,ixfmax2
6832 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6833 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6834 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6835 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6836 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6837 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6840 do ix2=ixfmin2,ixfmax2
6841 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6842 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6843 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6844 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6845 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6846 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6847 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6848 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6849 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6850 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6851 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6852 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6853 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6854 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6855 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6856 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6857 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6858 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6875 do ix3=ixfmax3,ixfmin3,-1
6876 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6877 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6878 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6879 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6880 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6881 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6884 do ix3=ixfmax3,ixfmin3,-1
6885 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6886 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6887 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6888 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6889 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6890 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6891 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6892 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6893 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6894 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6895 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6896 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6897 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6898 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6899 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6900 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6901 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6902 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6917 do ix3=ixfmin3,ixfmax3
6918 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6919 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6920 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6921 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6922 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6923 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6926 do ix3=ixfmin3,ixfmax3
6927 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6928 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6929 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6930 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6931 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6932 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6933 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6934 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6935 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6936 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6937 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6938 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6939 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6940 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6941 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6942 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6943 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6944 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6950 call mpistop(
"Special boundary is not defined for this region")
6953 end subroutine fixdivb_boundary
6962 double precision,
intent(in) :: qdt
6963 double precision,
intent(in) :: qt
6964 logical,
intent(inout) :: active
6967 integer,
parameter :: max_its = 50
6968 double precision :: residual_it(max_its), max_divb
6969 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6970 double precision :: res
6971 double precision,
parameter :: max_residual = 1
d-3
6972 double precision,
parameter :: residual_reduction = 1
d-10
6973 integer :: iigrid, igrid
6974 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6977 mg%operator_type = mg_laplacian
6985 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6986 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6989 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6990 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6992 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6993 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6996 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6997 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7001 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
7002 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
7003 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
7011 do iigrid = 1, igridstail
7012 igrid = igrids(iigrid);
7015 lvl =
mg%boxes(id)%lvl
7016 nc =
mg%box_size_lvl(lvl)
7022 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
7024 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
7025 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
7030 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
7033 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
7034 if (
mype == 0) print *,
"iteration vs residual"
7037 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
7038 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
7039 if (residual_it(n) < residual_reduction * max_divb)
exit
7041 if (
mype == 0 .and. n > max_its)
then
7042 print *,
"divb_multigrid warning: not fully converged"
7043 print *,
"current amplitude of divb: ", residual_it(max_its)
7044 print *,
"multigrid smallest grid: ", &
7045 mg%domain_size_lvl(:,
mg%lowest_lvl)
7046 print *,
"note: smallest grid ideally has <= 8 cells"
7047 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
7048 print *,
"note: dx/dy/dz should be similar"
7052 call mg_fas_vcycle(
mg, max_res=res)
7053 if (res < max_residual)
exit
7055 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
7060 do iigrid = 1, igridstail
7061 igrid = igrids(iigrid);
7070 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
7074 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
7076 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
7078 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
7091 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
7092 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
7095 if(total_energy)
then
7097 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
7100 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
7110 subroutine mhd_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7114 integer,
intent(in) :: ixi^
l, ixo^
l
7115 double precision,
intent(in) :: qt,qdt
7117 double precision,
intent(in) :: wp(ixi^s,1:nw)
7118 type(state) :: sct, s
7119 type(ct_velocity) :: vcts
7120 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7121 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7123 double precision :: circ(ixi^s,1:
ndim)
7125 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7126 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
7127 integer :: idim1,idim2,idir,iwdim1,iwdim2
7129 associate(bfaces=>s%ws,x=>s%x)
7136 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7143 i1kr^
d=
kr(idim1,^
d);
7146 i2kr^
d=
kr(idim2,^
d);
7149 if (
lvc(idim1,idim2,idir)==1)
then
7151 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7153 {
do ix^db=ixcmin^db,ixcmax^db\}
7154 fe(ix^
d,idir)=quarter*&
7155 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
7156 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
7158 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
7163 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
7171 if(
associated(usr_set_electric_field)) &
7172 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7174 circ(ixi^s,1:ndim)=zero
7179 ixcmin^d=ixomin^d-kr(idim1,^d);
7181 ixa^l=ixc^l-kr(idim2,^d);
7184 if(lvc(idim1,idim2,idir)==1)
then
7186 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7189 else if(lvc(idim1,idim2,idir)==-1)
then
7191 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7197 {
do ix^db=ixcmin^db,ixcmax^db\}
7199 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7201 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7208 end subroutine mhd_update_faces_average
7211 subroutine mhd_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7216 integer,
intent(in) :: ixi^
l, ixo^
l
7217 double precision,
intent(in) :: qt, qdt
7219 double precision,
intent(in) :: wp(ixi^s,1:nw)
7220 type(state) :: sct, s
7221 type(ct_velocity) :: vcts
7222 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7223 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7225 double precision :: circ(ixi^s,1:
ndim)
7227 double precision :: ecc(ixi^s,
sdim:3)
7228 double precision :: ein(ixi^s,
sdim:3)
7230 double precision :: el(ixi^s),er(ixi^s)
7232 double precision :: elc,erc
7234 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7236 double precision :: jce(ixi^s,
sdim:3)
7238 double precision :: xs(ixgs^t,1:
ndim)
7239 double precision :: gradi(ixgs^t)
7240 integer :: ixc^
l,ixa^
l
7241 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7243 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7246 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7252 {
do ix^db=iximin^db,iximax^db\}
7255 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_)
7256 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_)
7257 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_)
7260 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7267 {
do ix^db=iximin^db,iximax^db\}
7270 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7271 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7272 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7275 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7289 i1kr^d=kr(idim1,^d);
7292 i2kr^d=kr(idim2,^d);
7295 if (lvc(idim1,idim2,idir)==1)
then
7297 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7300 {
do ix^db=ixcmin^db,ixcmax^db\}
7301 fe(ix^d,idir)=quarter*&
7302 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7303 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7308 ixamax^d=ixcmax^d+i1kr^d;
7309 {
do ix^db=ixamin^db,ixamax^db\}
7310 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7311 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7314 do ix^db=ixcmin^db,ixcmax^db\}
7315 if(vnorm(ix^d,idim1)>0.d0)
then
7317 else if(vnorm(ix^d,idim1)<0.d0)
then
7318 elc=el({ix^d+i1kr^d})
7320 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7322 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7324 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7325 erc=er({ix^d+i1kr^d})
7327 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7329 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7334 ixamax^d=ixcmax^d+i2kr^d;
7335 {
do ix^db=ixamin^db,ixamax^db\}
7336 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7337 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7340 do ix^db=ixcmin^db,ixcmax^db\}
7341 if(vnorm(ix^d,idim2)>0.d0)
then
7343 else if(vnorm(ix^d,idim2)<0.d0)
then
7344 elc=el({ix^d+i2kr^d})
7346 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7348 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7350 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7351 erc=er({ix^d+i2kr^d})
7353 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7355 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7359 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7364 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7378 if (lvc(idim1,idim2,idir)==0) cycle
7380 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7381 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7384 xs(ixa^s,:)=x(ixa^s,:)
7385 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7386 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7387 if (lvc(idim1,idim2,idir)==1)
then
7388 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7390 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7397 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7399 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7403 {
do ix^db=ixomin^db,ixomax^db\}
7404 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7405 +ein(ix1,ix2-1,ix3-1,idir))
7406 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7407 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7409 else if(idir==2)
then
7410 {
do ix^db=ixomin^db,ixomax^db\}
7411 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7412 +ein(ix1-1,ix2,ix3-1,idir))
7413 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7414 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7417 {
do ix^db=ixomin^db,ixomax^db\}
7418 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7419 +ein(ix1-1,ix2-1,ix3,idir))
7420 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7421 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7427 {
do ix^db=ixomin^db,ixomax^db\}
7428 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7429 +ein(ix1-1,ix2-1,idir))
7430 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7431 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7436 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7442 if(
associated(usr_set_electric_field)) &
7443 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7445 circ(ixi^s,1:ndim)=zero
7450 ixcmin^d=ixomin^d-kr(idim1,^d);
7452 ixa^l=ixc^l-kr(idim2,^d);
7455 if(lvc(idim1,idim2,idir)==1)
then
7457 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7460 else if(lvc(idim1,idim2,idir)==-1)
then
7462 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7468 {
do ix^db=ixcmin^db,ixcmax^db\}
7470 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7472 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7479 end subroutine mhd_update_faces_contact
7482 subroutine mhd_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7487 integer,
intent(in) :: ixi^
l, ixo^
l
7488 double precision,
intent(in) :: qt, qdt
7490 double precision,
intent(in) :: wp(ixi^s,1:nw)
7491 type(state) :: sct, s
7492 type(ct_velocity) :: vcts
7493 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7494 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7496 double precision :: vtill(ixi^s,2)
7497 double precision :: vtilr(ixi^s,2)
7498 double precision :: bfacetot(ixi^s,
ndim)
7499 double precision :: btill(ixi^s,
ndim)
7500 double precision :: btilr(ixi^s,
ndim)
7501 double precision :: cp(ixi^s,2)
7502 double precision :: cm(ixi^s,2)
7503 double precision :: circ(ixi^s,1:
ndim)
7505 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7506 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7507 integer :: idim1,idim2,idir,ix^
d
7509 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7510 cbarmax=>vcts%cbarmax)
7523 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7539 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7543 idim2=mod(idir+1,3)+1
7545 jxc^
l=ixc^
l+
kr(idim1,^
d);
7546 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7550 vtill(ixi^s,2),vtilr(ixi^s,2))
7553 vtill(ixi^s,1),vtilr(ixi^s,1))
7559 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7560 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7562 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7563 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7566 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7569 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7573 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7574 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7576 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7577 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7581 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7582 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7583 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7584 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7585 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7586 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7587 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7588 /(cp(ixc^s,2)+cm(ixc^s,2))
7591 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7595 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7609 circ(ixi^s,1:
ndim)=zero
7614 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7618 if(
lvc(idim1,idim2,idir)/=0)
then
7619 hxc^
l=ixc^
l-
kr(idim2,^
d);
7621 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7622 +
lvc(idim1,idim2,idir)&
7628 {
do ix^db=ixcmin^db,ixcmax^db\}
7630 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7632 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7638 end subroutine mhd_update_faces_hll
7641 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
7646 integer,
intent(in) :: ixi^
l, ixo^
l
7648 double precision,
intent(in) :: wp(ixi^s,1:nw)
7649 type(state),
intent(in) :: sct, s
7651 double precision :: jce(ixi^s,
sdim:3)
7654 double precision :: jcc(ixi^s,7-2*
ndir:3)
7656 double precision :: xs(ixgs^t,1:
ndim)
7658 double precision :: eta(ixi^s)
7659 double precision :: gradi(ixgs^t)
7660 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7662 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7668 if (
lvc(idim1,idim2,idir)==0) cycle
7670 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7671 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7674 xs(ixb^s,:)=x(ixb^s,:)
7675 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7676 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7677 if (
lvc(idim1,idim2,idir)==1)
then
7678 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7680 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7687 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7695 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7696 jcc(ixc^s,idir)=0.d0
7698 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7699 ixamin^
d=ixcmin^
d+ix^
d;
7700 ixamax^
d=ixcmax^
d+ix^
d;
7701 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7703 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7704 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7709 end subroutine get_resistive_electric_field
7712 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7715 integer,
intent(in) :: ixi^
l, ixo^
l
7716 double precision,
intent(in) :: w(ixi^s,1:nw)
7717 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7718 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7720 double precision :: jxbxb(ixi^s,1:3)
7721 integer :: idir,ixa^
l,ixc^
l,ix^
d
7724 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7731 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7734 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7735 ixamin^
d=ixcmin^
d+ix^
d;
7736 ixamax^
d=ixcmax^
d+ix^
d;
7737 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7739 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7742 end subroutine get_ambipolar_electric_field
7748 integer,
intent(in) :: ixo^
l
7758 do ix^db=ixomin^db,ixomax^db\}
7760 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7761 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7762 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7763 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7764 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7765 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7768 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7769 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7770 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7771 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7814 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7815 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7816 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7818 double precision :: adummy(ixis^s,1:3)
7824 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7827 integer,
intent(in) :: ixi^
l, ixo^
l
7828 double precision,
intent(in) :: w(ixi^s,1:nw)
7829 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7830 double precision,
intent(out):: rfactor(ixi^s)
7832 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7836 rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/(2.d0+3.d0*
he_abundance)
7838 end subroutine rfactor_from_temperature_ionization
7840 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7842 integer,
intent(in) :: ixi^
l, ixo^
l
7843 double precision,
intent(in) :: w(ixi^s,1:nw)
7844 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7845 double precision,
intent(out):: rfactor(ixi^s)
7849 end subroutine rfactor_from_constant_ionization
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
subroutine cak_get_dt(w, ixil, ixol, dtnew, dxd, x)
Check time step for total radiation contribution.
subroutine cak_init(phys_gamma)
Initialize the module.
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)
Module for flux conservation near refinement boundaries.
subroutine, public store_flux(igrid, fc, idimlim, nwfluxin)
subroutine, public store_edge(igrid, ixil, fe, idimlim)
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)
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 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.
logical b0fieldalloccoarse
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision unit_mass
Physical scaling factor for mass.
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 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.
double precision, dimension(:), allocatable, parameter d
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.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
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
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(w, 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 ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
subroutine ionization_degree_init()
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
integer, public, protected c_
logical, public, protected mhd_gravity
Whether gravity is added.
logical, public has_equi_rho0
whether split off equilibrium density
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.
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)
logical, public, protected mhd_hyperbolic_thermal_conduction
Whether thermal conduction is used.
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.
logical, public, protected mhd_partial_ionization
Whether plasma is partially ionized.
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
double precision, public, protected rr
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
double precision, public mhd_gamma
The adiabatic index.
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 q_
Index of the heat flux q.
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 CAK radiation line force is activated.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
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.
double precision, public hypertc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
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.
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.
procedure(sub_get_pthermal), pointer, public mhd_get_rfactor
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.
logical, public, protected mhd_htc_sat
Whether saturation is considered for hyperbolic TC.
logical, public, protected mhd_glm
Whether GLM-MHD is used to control div B.
logical, public clean_initial_divb
clean initial divB
procedure(sub_convert), pointer, public mhd_to_conserved
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.
procedure(sub_get_pthermal), pointer, public mhd_get_temperature
integer, public equi_pe0_
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+)
procedure(sub_convert), pointer, public mhd_to_primitive
logical, public has_equi_pe0
whether split off equilibrium thermal pressure
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.
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_hydrodynamic_e
Whether hydrodynamic energy is solved instead of total energy.
subroutine, public mhd_phys_init()
logical, public, protected mhd_trac
Whether TRAC method is used.
logical, public, protected eq_state_units
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.
procedure(sub_get_pthermal), pointer, public mhd_get_pthermal
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)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected mhd_4th_order
MHD fourth order.
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 for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
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)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
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 1) in amrvac.par in sts_list set the following parameters which have...
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)
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 viscosity_init(phys_wider_stencil)
Initialize the module.
subroutine viscosity_get_dt(w, ixil, ixol, dtnew, dxd, x)
procedure(sub_add_source), pointer, public viscosity_add_source
The data structure that contains information about a tree node/grid block.