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)
2571 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2572 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2573 double precision,
intent(inout) :: cmax(ixi^s)
2575 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2581 {
do ix^db=ixomin^db,ixomax^db \}
2592 cfast2=b2*inv_rho+cmax(ix^
d)
2593 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2594 if(avmincs2<zero) avmincs2=zero
2595 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2599 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2601 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2604 {
do ix^db=ixomin^db,ixomax^db \}
2614 b2=(^
c&w(ix^d,
b^
c_)**2+)
2615 cfast2=b2*inv_rho+cmax(ix^d)
2616 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2617 if(avmincs2<zero) avmincs2=zero
2618 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2622 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2624 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2628 end subroutine mhd_get_cmax_origin_noe
2631 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2634 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2635 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2636 double precision,
intent(inout):: cmax(ixi^s)
2638 double precision :: csound, avmincs2, idim_alfven_speed2
2639 double precision :: inv_rho, alfven_speed2, gamma2
2642 {
do ix^db=ixomin^db,ixomax^db \}
2643 inv_rho=1.d0/w(ix^
d,
rho_)
2644 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2645 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2646 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2649 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2652 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2653 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2654 if(avmincs2<zero) avmincs2=zero
2656 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2657 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2660 end subroutine mhd_get_cmax_semirelati
2663 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2666 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2667 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2668 double precision,
intent(inout):: cmax(ixi^s)
2670 double precision :: csound, avmincs2, idim_alfven_speed2
2671 double precision :: inv_rho, alfven_speed2, gamma2
2674 {
do ix^db=ixomin^db,ixomax^db \}
2675 inv_rho=1.d0/w(ix^
d,
rho_)
2676 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2677 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2678 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2680 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2683 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2684 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2685 if(avmincs2<zero) avmincs2=zero
2687 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2688 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2691 end subroutine mhd_get_cmax_semirelati_noe
2693 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2696 integer,
intent(in) :: ixi^
l, ixo^
l
2697 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2698 double precision,
intent(inout) :: a2max(
ndim)
2699 double precision :: a2(ixi^s,
ndim,nw)
2700 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2705 hxo^
l=ixo^
l-
kr(i,^
d);
2706 gxo^
l=hxo^
l-
kr(i,^
d);
2707 jxo^
l=ixo^
l+
kr(i,^
d);
2708 kxo^
l=jxo^
l+
kr(i,^
d);
2709 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2710 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2711 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2713 end subroutine mhd_get_a2max
2716 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2719 integer,
intent(in) :: ixi^
l,ixo^
l
2720 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2722 double precision,
intent(inout) :: w(ixi^s,1:nw)
2723 double precision,
intent(out) :: tco_local,tmax_local
2725 double precision,
parameter :: trac_delta=0.25d0
2726 double precision :: te(ixi^s),lts(ixi^s)
2727 double precision,
dimension(1:ndim) :: bdir, bunitvec
2728 double precision,
dimension(ixI^S,1:ndim) :: gradt
2729 double precision :: ltrc,ltrp,altr
2730 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2731 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2734 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2737 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2740 tmax_local=maxval(te(ixo^s))
2748 do ix1=ixomin1,ixomax1
2749 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2750 if(lts(ix1)>trac_delta)
then
2751 tco_local=max(tco_local,te(ix1))
2763 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2764 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2765 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2766 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2768 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2779 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2786 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2791 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2792 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2796 if(bdir(1)/=0.d0)
then
2797 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2799 block%special_values(3)=0.d0
2801 if(bdir(2)/=0.d0)
then
2802 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2804 block%special_values(4)=0.d0
2808 if(bdir(1)/=0.d0)
then
2809 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2810 (bdir(3)/bdir(1))**2)
2812 block%special_values(3)=0.d0
2814 if(bdir(2)/=0.d0)
then
2815 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2816 (bdir(3)/bdir(2))**2)
2818 block%special_values(4)=0.d0
2820 if(bdir(3)/=0.d0)
then
2821 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2822 (bdir(2)/bdir(3))**2)
2824 block%special_values(5)=0.d0
2829 block%special_values(1)=zero
2830 {
do ix^db=ixomin^db,ixomax^db\}
2832 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2834 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2837 if(bdir(1)/=0.d0)
then
2838 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2842 if(bdir(2)/=0.d0)
then
2843 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2848 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2849 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2852 if(bdir(1)/=0.d0)
then
2853 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2857 if(bdir(2)/=0.d0)
then
2858 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2862 if(bdir(3)/=0.d0)
then
2863 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2868 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2869 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2871 if(lts(ix^d)>trac_delta)
then
2872 block%special_values(1)=max(block%special_values(1),te(ix^d))
2875 block%special_values(2)=tmax_local
2894 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2895 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2896 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2900 {
do ix^db=ixpmin^db,ixpmax^db\}
2901 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2903 if(bdir(1)/=0.d0)
then
2904 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2908 if(bdir(2)/=0.d0)
then
2909 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2915 if(bdir(1)/=0.d0)
then
2916 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2920 if(bdir(2)/=0.d0)
then
2921 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2925 if(bdir(3)/=0.d0)
then
2926 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2932 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2934 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2935 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2938 {
do ix^db=ixpmin^db,ixpmax^db\}
2940 if(w(ix^d,iw_mag(1))/=0.d0)
then
2941 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)
2945 if(w(ix^d,iw_mag(2))/=0.d0)
then
2946 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)
2952 if(w(ix^d,iw_mag(1))/=0.d0)
then
2953 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+&
2954 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
2958 if(w(ix^d,iw_mag(2))/=0.d0)
then
2959 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+&
2960 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
2964 if(w(ix^d,iw_mag(3))/=0.d0)
then
2965 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+&
2966 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
2972 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2974 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2975 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2981 {
do ix^db=ixpmin^db,ixpmax^db\}
2983 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
2984 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
2985 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2988 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
2989 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
2990 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
2991 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2997 call mpistop(
"unknown mhd_trac_type")
3000 end subroutine mhd_get_tcutoff
3003 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
3006 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3007 double precision,
intent(in) :: wprim(ixi^s, nw)
3008 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3009 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
3011 double precision :: csound(ixi^s,
ndim)
3012 double precision,
allocatable :: tmp(:^
d&)
3013 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
3017 allocate(tmp(ixa^s))
3019 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
3020 csound(ixa^s,id)=tmp(ixa^s)
3023 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
3024 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
3025 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
3026 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))
3030 ixamax^
d=ixcmax^
d+
kr(id,^
d);
3031 ixamin^
d=ixcmin^
d+
kr(id,^
d);
3032 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)))
3033 ixamax^
d=ixcmax^
d-
kr(id,^
d);
3034 ixamin^
d=ixcmin^
d-
kr(id,^
d);
3035 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)))
3040 ixamax^
d=jxcmax^
d+
kr(id,^
d);
3041 ixamin^
d=jxcmin^
d+
kr(id,^
d);
3042 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)))
3043 ixamax^
d=jxcmax^
d-
kr(id,^
d);
3044 ixamin^
d=jxcmin^
d-
kr(id,^
d);
3045 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)))
3049 end subroutine mhd_get_h_speed
3052 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3055 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3056 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3057 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3058 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3059 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3060 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3061 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3063 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3064 double precision :: umean, dmean, tmp1, tmp2, tmp3
3071 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3072 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3073 if(
present(cmin))
then
3074 {
do ix^db=ixomin^db,ixomax^db\}
3075 tmp1=sqrt(wlp(ix^
d,
rho_))
3076 tmp2=sqrt(wrp(ix^
d,
rho_))
3077 tmp3=1.d0/(tmp1+tmp2)
3078 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3079 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3080 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3081 cmin(ix^
d,1)=umean-dmean
3082 cmax(ix^
d,1)=umean+dmean
3084 if(h_correction)
then
3085 {
do ix^db=ixomin^db,ixomax^db\}
3086 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3087 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3091 {
do ix^db=ixomin^db,ixomax^db\}
3092 tmp1=sqrt(wlp(ix^d,
rho_))
3093 tmp2=sqrt(wrp(ix^d,
rho_))
3094 tmp3=1.d0/(tmp1+tmp2)
3095 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3096 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3097 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3098 cmax(ix^d,1)=abs(umean)+dmean
3102 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3103 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3104 if(
present(cmin))
then
3105 {
do ix^db=ixomin^db,ixomax^db\}
3106 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3107 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
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 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3120 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3121 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3122 if(
present(cmin))
then
3123 {
do ix^db=ixomin^db,ixomax^db\}
3124 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3125 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3126 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3128 if(h_correction)
then
3129 {
do ix^db=ixomin^db,ixomax^db\}
3130 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3131 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3135 {
do ix^db=ixomin^db,ixomax^db\}
3136 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3137 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3142 end subroutine mhd_get_cbounds
3145 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3148 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3149 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3150 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3151 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3152 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3153 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3154 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3156 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3161 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3162 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3164 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3165 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3167 if(
present(cmin))
then
3168 {
do ix^db=ixomin^db,ixomax^db\}
3169 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3170 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3171 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3174 {
do ix^db=ixomin^db,ixomax^db\}
3175 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3176 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3180 end subroutine mhd_get_cbounds_semirelati
3183 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3186 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3187 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3188 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3189 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3190 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3191 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3192 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3194 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3195 double precision :: umean, dmean, tmp1, tmp2, tmp3
3202 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3203 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3204 if(
present(cmin))
then
3205 {
do ix^db=ixomin^db,ixomax^db\}
3208 tmp3=1.d0/(tmp1+tmp2)
3209 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3210 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3211 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3212 cmin(ix^
d,1)=umean-dmean
3213 cmax(ix^
d,1)=umean+dmean
3215 if(h_correction)
then
3216 {
do ix^db=ixomin^db,ixomax^db\}
3217 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3218 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3222 {
do ix^db=ixomin^db,ixomax^db\}
3225 tmp3=1.d0/(tmp1+tmp2)
3226 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3227 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3228 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3229 cmax(ix^d,1)=abs(umean)+dmean
3233 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3234 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3235 if(
present(cmin))
then
3236 {
do ix^db=ixomin^db,ixomax^db\}
3237 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3238 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
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 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3251 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3252 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3253 if(
present(cmin))
then
3254 {
do ix^db=ixomin^db,ixomax^db\}
3255 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3256 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3257 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3259 if(h_correction)
then
3260 {
do ix^db=ixomin^db,ixomax^db\}
3261 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3262 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3266 {
do ix^db=ixomin^db,ixomax^db\}
3267 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3268 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3273 end subroutine mhd_get_cbounds_split_rho
3276 subroutine mhd_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3279 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3280 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3281 double precision,
intent(in) :: cmax(ixi^s)
3282 double precision,
intent(in),
optional :: cmin(ixi^s)
3283 type(ct_velocity),
intent(inout):: vcts
3285 end subroutine mhd_get_ct_velocity_average
3287 subroutine mhd_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3290 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3291 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3292 double precision,
intent(in) :: cmax(ixi^s)
3293 double precision,
intent(in),
optional :: cmin(ixi^s)
3294 type(ct_velocity),
intent(inout):: vcts
3296 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3298 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3300 end subroutine mhd_get_ct_velocity_contact
3302 subroutine mhd_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3305 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3306 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3307 double precision,
intent(in) :: cmax(ixi^s)
3308 double precision,
intent(in),
optional :: cmin(ixi^s)
3309 type(ct_velocity),
intent(inout):: vcts
3311 integer :: idime,idimn
3313 if(.not.
allocated(vcts%vbarC))
then
3314 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3315 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3318 if(
present(cmin))
then
3319 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3320 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3322 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3323 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3326 idimn=mod(idim,
ndir)+1
3327 idime=mod(idim+1,
ndir)+1
3329 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3330 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3331 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3332 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3333 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3335 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3336 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3337 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3338 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3339 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3341 end subroutine mhd_get_ct_velocity_hll
3344 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3347 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3348 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3349 double precision,
intent(out):: csound(ixo^s)
3351 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3358 {
do ix^db=ixomin^db,ixomax^db \}
3359 inv_rho=1.d0/w(ix^
d,
rho_)
3366 cfast2=b2*inv_rho+csound(ix^
d)
3367 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3369 if(avmincs2<zero) avmincs2=zero
3370 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3372 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3376 {
do ix^db=ixomin^db,ixomax^db \}
3377 inv_rho=1.d0/w(ix^d,
rho_)
3383 b2=(^
c&w(ix^d,
b^
c_)**2+)
3384 cfast2=b2*inv_rho+csound(ix^d)
3385 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3386 if(avmincs2<zero) avmincs2=zero
3387 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3389 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3394 end subroutine mhd_get_csound_prim
3397 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3400 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3401 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3402 double precision,
intent(out):: csound(ixo^s)
3404 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3411 {
do ix^db=ixomin^db,ixomax^db \}
3418 cfast2=b2*inv_rho+csound(ix^
d)
3419 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3421 if(avmincs2<zero) avmincs2=zero
3422 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3424 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3428 {
do ix^db=ixomin^db,ixomax^db \}
3434 b2=(^
c&w(ix^d,
b^
c_)**2+)
3435 cfast2=b2*inv_rho+csound(ix^d)
3436 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3437 if(avmincs2<zero) avmincs2=zero
3438 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3440 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3445 end subroutine mhd_get_csound_prim_split
3448 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3451 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3453 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3454 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3456 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3459 {
do ix^db=ixomin^db,ixomax^db\}
3460 inv_rho = 1.d0/w(ix^
d,
rho_)
3463 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3464 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3465 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3466 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3469 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3470 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3471 if(avmincs2<zero) avmincs2=zero
3473 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3476 end subroutine mhd_get_csound_semirelati
3479 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3482 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3484 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3485 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3487 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3490 {
do ix^db=ixomin^db,ixomax^db\}
3491 inv_rho = 1.d0/w(ix^
d,
rho_)
3494 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3495 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3496 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3497 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3500 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3501 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3502 if(avmincs2<zero) avmincs2=zero
3504 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3507 end subroutine mhd_get_csound_semirelati_noe
3510 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3513 integer,
intent(in) :: ixi^
l, ixo^
l
3514 double precision,
intent(in) :: w(ixi^s,nw)
3515 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3516 double precision,
intent(out):: pth(ixi^s)
3524 end subroutine mhd_get_pthermal_noe
3527 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3531 integer,
intent(in) :: ixi^
l, ixo^
l
3532 double precision,
intent(in) :: w(ixi^s,nw)
3533 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3534 double precision,
intent(out):: pth(ixi^s)
3538 {
do ix^db= ixomin^db,ixomax^db\}
3542 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3547 if(check_small_values.and..not.fix_small_values)
then
3548 {
do ix^db= ixomin^db,ixomax^db\}
3549 if(pth(ix^d)<small_pressure)
then
3550 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3551 " encountered when call mhd_get_pthermal_inte"
3552 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3553 write(*,*)
"Location: ", x(ix^d,:)
3554 write(*,*)
"Cell number: ", ix^d
3556 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3559 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3560 write(*,*)
"Saving status at the previous time step"
3566 end subroutine mhd_get_pthermal_inte
3569 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3573 integer,
intent(in) :: ixi^
l, ixo^
l
3574 double precision,
intent(in) :: w(ixi^s,nw)
3575 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3576 double precision,
intent(out):: pth(ixi^s)
3580 {
do ix^db=ixomin^db,ixomax^db\}
3583 +(^
c&w(ix^
d,
b^
c_)**2+)))
3585 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3586 +(^
c&w(ix^
d,
b^
c_)**2+)))
3594 if(check_small_values.and..not.fix_small_values)
then
3595 {
do ix^db=ixomin^db,ixomax^db\}
3596 if(pth(ix^d)<small_pressure)
then
3597 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3598 " encountered when call mhd_get_pthermal"
3599 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3600 write(*,*)
"Location: ", x(ix^d,:)
3601 write(*,*)
"Cell number: ", ix^d
3603 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3606 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3607 write(*,*)
"Saving status at the previous time step"
3613 end subroutine mhd_get_pthermal_origin
3616 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3620 integer,
intent(in) :: ixi^
l, ixo^
l
3621 double precision,
intent(in) :: w(ixi^s,nw)
3622 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3623 double precision,
intent(out):: pth(ixi^s)
3625 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3628 {
do ix^db=ixomin^db,ixomax^db\}
3629 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3630 if(b2>smalldouble)
then
3638 inv_rho=1.d0/w(ix^
d,
rho_)
3640 b2=b2*inv_rho*inv_squared_c
3642 gamma2=1.d0/(1.d0+b2)
3644 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3648 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3649 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3650 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3654 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3660 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3661 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3662 +(^
c&w(ix^
d,
b^
c_)**2+)&
3663 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3667 if(check_small_values.and..not.fix_small_values)
then
3668 {
do ix^db=ixomin^db,ixomax^db\}
3669 if(pth(ix^d)<small_pressure)
then
3670 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3671 " encountered when call mhd_get_pthermal_semirelati"
3672 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3673 write(*,*)
"Location: ", x(ix^d,:)
3674 write(*,*)
"Cell number: ", ix^d
3676 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3679 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3680 write(*,*)
"Saving status at the previous time step"
3686 end subroutine mhd_get_pthermal_semirelati
3689 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3693 integer,
intent(in) :: ixi^
l, ixo^
l
3694 double precision,
intent(in) :: w(ixi^s,nw)
3695 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3696 double precision,
intent(out):: pth(ixi^s)
3700 {
do ix^db= ixomin^db,ixomax^db\}
3701 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3704 if(check_small_values.and..not.fix_small_values)
then
3705 {
do ix^db= ixomin^db,ixomax^db\}
3706 if(pth(ix^d)<small_pressure)
then
3707 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3708 " encountered when call mhd_get_pthermal_hde"
3709 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3710 write(*,*)
"Location: ", x(ix^d,:)
3711 write(*,*)
"Cell number: ", ix^d
3713 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3716 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3717 write(*,*)
"Saving status at the previous time step"
3723 end subroutine mhd_get_pthermal_hde
3726 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3728 integer,
intent(in) :: ixi^
l, ixo^
l
3729 double precision,
intent(in) :: w(ixi^s, 1:nw)
3730 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3731 double precision,
intent(out):: res(ixi^s)
3732 res(ixo^s) = w(ixo^s,
te_)
3733 end subroutine mhd_get_temperature_from_te
3736 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3738 integer,
intent(in) :: ixi^
l, ixo^
l
3739 double precision,
intent(in) :: w(ixi^s, 1:nw)
3740 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3741 double precision,
intent(out):: res(ixi^s)
3743 double precision :: r(ixi^s)
3746 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3747 end subroutine mhd_get_temperature_from_eint
3750 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3752 integer,
intent(in) :: ixi^
l, ixo^
l
3753 double precision,
intent(in) :: w(ixi^s, 1:nw)
3754 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3755 double precision,
intent(out):: res(ixi^s)
3757 double precision :: r(ixi^s),rho(ixi^s)
3762 res(ixo^s)=res(ixo^s)/(r(ixo^s)*rho(ixo^s))
3764 end subroutine mhd_get_temperature_from_etot
3766 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3768 integer,
intent(in) :: ixi^
l, ixo^
l
3769 double precision,
intent(in) :: w(ixi^s, 1:nw)
3770 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3771 double precision,
intent(out):: res(ixi^s)
3773 double precision :: r(ixi^s)
3779 end subroutine mhd_get_temperature_from_eint_with_equi
3781 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3783 integer,
intent(in) :: ixi^
l, ixo^
l
3784 double precision,
intent(in) :: w(ixi^s, 1:nw)
3785 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3786 double precision,
intent(out):: res(ixi^s)
3788 double precision :: r(ixi^s)
3794 end subroutine mhd_get_temperature_equi
3796 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3798 integer,
intent(in) :: ixi^
l, ixo^
l
3799 double precision,
intent(in) :: w(ixi^s, 1:nw)
3800 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3801 double precision,
intent(out):: res(ixi^s)
3803 end subroutine mhd_get_rho_equi
3805 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3807 integer,
intent(in) :: ixi^
l, ixo^
l
3808 double precision,
intent(in) :: w(ixi^s, 1:nw)
3809 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3810 double precision,
intent(out):: res(ixi^s)
3812 end subroutine mhd_get_pe_equi
3815 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3819 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3821 double precision,
intent(in) :: wc(ixi^s,nw)
3823 double precision,
intent(in) :: w(ixi^s,nw)
3824 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3825 double precision,
intent(out) :: f(ixi^s,nwflux)
3827 double precision :: vhall(ixi^s,1:
ndir)
3828 double precision :: ptotal
3832 {
do ix^db=ixomin^db,ixomax^db\}
3845 {
do ix^db=ixomin^db,ixomax^db\}
3849 ^
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_)\
3850 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3852 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3855 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3856 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3858 ^
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_)\
3862 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3863 {
do ix^db=ixomin^db,ixomax^db\}
3864 if(total_energy)
then
3866 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3867 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3870 ^
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))\
3874 {
do ix^db=ixomin^db,ixomax^db\}
3875 f(ix^d,mag(idim))=w(ix^d,
psi_)
3877 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3882 {
do ix^db=ixomin^db,ixomax^db\}
3888 {
do ix^db=ixomin^db,ixomax^db\}
3889 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)
3894 end subroutine mhd_get_flux
3897 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3901 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3903 double precision,
intent(in) :: wc(ixi^s,nw)
3905 double precision,
intent(in) :: w(ixi^s,nw)
3906 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3907 double precision,
intent(out) :: f(ixi^s,nwflux)
3909 double precision :: vhall(ixi^s,1:
ndir)
3912 {
do ix^db=ixomin^db,ixomax^db\}
3923 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3924 {
do ix^db=ixomin^db,ixomax^db\}
3926 ^
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))\
3930 {
do ix^db=ixomin^db,ixomax^db\}
3931 f(ix^d,mag(idim))=w(ix^d,
psi_)
3933 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3938 {
do ix^db=ixomin^db,ixomax^db\}
3943 end subroutine mhd_get_flux_noe
3946 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3950 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3952 double precision,
intent(in) :: wc(ixi^s,nw)
3954 double precision,
intent(in) :: w(ixi^s,nw)
3955 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3956 double precision,
intent(out) :: f(ixi^s,nwflux)
3958 double precision :: vhall(ixi^s,1:
ndir)
3961 {
do ix^db=ixomin^db,ixomax^db\}
3974 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3975 {
do ix^db=ixomin^db,ixomax^db\}
3977 ^
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))\
3981 {
do ix^db=ixomin^db,ixomax^db\}
3982 f(ix^d,mag(idim))=w(ix^d,
psi_)
3984 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3989 {
do ix^db=ixomin^db,ixomax^db\}
3995 {
do ix^db=ixomin^db,ixomax^db\}
3996 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)
4001 end subroutine mhd_get_flux_hde
4004 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
4008 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4010 double precision,
intent(in) :: wc(ixi^s,nw)
4012 double precision,
intent(in) :: w(ixi^s,nw)
4013 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4014 double precision,
intent(out) :: f(ixi^s,nwflux)
4016 double precision :: vhall(ixi^s,1:
ndir)
4017 double precision :: ptotal, btotal(ixo^s,1:
ndir)
4020 {
do ix^db=ixomin^db,ixomax^db\}
4028 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
4032 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
4036 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
4037 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4039 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
4043 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4046 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
4053 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
4054 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
4059 {
do ix^db=ixomin^db,ixomax^db\}
4060 f(ix^d,mag(idim))=w(ix^d,
psi_)
4062 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4067 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4068 {
do ix^db=ixomin^db,ixomax^db\}
4070 ^
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))\
4071 if(total_energy)
then
4073 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
4074 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
4080 {
do ix^db=ixomin^db,ixomax^db\}
4085 {
do ix^db=ixomin^db,ixomax^db\}
4086 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
4091 end subroutine mhd_get_flux_split
4094 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
4098 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4100 double precision,
intent(in) :: wc(ixi^s,nw)
4102 double precision,
intent(in) :: w(ixi^s,nw)
4103 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4104 double precision,
intent(out) :: f(ixi^s,nwflux)
4106 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
4109 {
do ix^db=ixomin^db,ixomax^db\}
4114 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4115 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4116 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4121 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4126 e2=(^
c&e(ix^
d,^
c)**2+)
4133 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4134 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4135 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4138 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4139 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4152 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4154 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)
4161 {
do ix^db=ixomin^db,ixomax^db\}
4162 f(ix^d,mag(idim))=w(ix^d,
psi_)
4164 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4169 {
do ix^db=ixomin^db,ixomax^db\}
4174 {
do ix^db=ixomin^db,ixomax^db\}
4175 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)
4180 end subroutine mhd_get_flux_semirelati
4182 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4186 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4188 double precision,
intent(in) :: wc(ixi^s,nw)
4190 double precision,
intent(in) :: w(ixi^s,nw)
4191 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4192 double precision,
intent(out) :: f(ixi^s,nwflux)
4194 double precision :: e(ixo^s,1:
ndir),e2
4197 {
do ix^db=ixomin^db,ixomax^db\}
4202 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4203 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4204 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4205 e2=(^
c&e(ix^
d,^
c)**2+)
4210 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4219 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4228 {
do ix^db=ixomin^db,ixomax^db\}
4229 f(ix^d,mag(idim))=w(ix^d,
psi_)
4231 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4236 {
do ix^db=ixomin^db,ixomax^db\}
4241 end subroutine mhd_get_flux_semirelati_noe
4247 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4249 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4250 double precision,
intent(in) :: qdt
4251 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4252 double precision,
intent(inout) :: w(ixi^s,1:nw)
4253 double precision :: tmp(ixi^s)
4254 double precision :: jxbxb(ixi^s,1:3)
4256 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4259 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4261 end subroutine add_source_ambipolar_internal_energy
4263 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4266 integer,
intent(in) :: ixi^
l, ixo^
l
4267 double precision,
intent(in) :: w(ixi^s,nw)
4268 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4269 double precision,
intent(out) :: res(:^
d&,:)
4271 double precision :: btot(ixi^s,1:3)
4272 double precision :: current(ixi^s,7-2*
ndir:3)
4273 double precision :: tmp(ixi^s),b2(ixi^s)
4274 integer :: idir, idirmin
4283 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4286 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4289 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4290 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4292 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4295 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4297 end subroutine mhd_get_jxbxb
4304 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4308 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4309 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4310 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4311 double precision,
intent(in) :: my_dt
4312 logical,
intent(in) :: fix_conserve_at_step
4314 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4315 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4316 double precision :: fe(ixi^s,
sdim:3)
4317 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4318 integer :: i, ixa^
l, ie_
4324 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4334 btot(ixa^s,1:3)=0.d0
4340 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4343 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4344 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4346 wres(ixo^s,
e_)=-tmp2(ixo^s)
4352 ff(ixa^s,1) = tmp(ixa^s,2)
4353 ff(ixa^s,2) = -tmp(ixa^s,1)
4355 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4356 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4357 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4360 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4362 ixamin^
d=ixomin^
d-1;
4363 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4372 ff(ixa^s,2) = tmp(ixa^s,3)
4373 ff(ixa^s,3) = -tmp(ixa^s,2)
4374 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4375 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4377 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4379 ff(ixa^s,1) = -tmp(ixa^s,3)
4381 ff(ixa^s,3) = tmp(ixa^s,1)
4382 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4383 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4384 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4388 ff(ixa^s,1) = tmp(ixa^s,2)
4389 ff(ixa^s,2) = -tmp(ixa^s,1)
4391 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4392 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4393 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4398 if(fix_conserve_at_step)
then
4399 fluxall=my_dt*fluxall
4406 end subroutine sts_set_source_ambipolar
4409 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4412 integer,
intent(in) :: ixi^
l, ixo^
l
4413 double precision,
intent(in) :: w(ixi^s,1:nw)
4414 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4416 double precision,
intent(in) :: ecc(ixi^s,1:3)
4417 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4418 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4420 integer :: hxc^
l,ixc^
l,ixa^
l
4421 integer :: idim1,idim2,idir,ix^
d
4427 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4429 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4430 ixamin^
d=ixcmin^
d+ix^
d;
4431 ixamax^
d=ixcmax^
d+ix^
d;
4432 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4434 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4440 ixcmin^d=ixomin^d-1;
4448 hxc^l=ixc^l-kr(idim2,^d);
4450 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4451 +lvc(idim1,idim2,idir)&
4456 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4459 end subroutine update_faces_ambipolar
4463 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4466 integer,
intent(in) :: ixi^
l, ixo^
l
4467 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4468 double precision,
intent(out) :: src(ixi^s)
4470 double precision :: ffc(ixi^s,1:
ndim)
4471 double precision :: dxinv(
ndim)
4472 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4478 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4480 ixbmin^
d=ixcmin^
d+ix^
d;
4481 ixbmax^
d=ixcmax^
d+ix^
d;
4484 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4486 ff(ixi^s,1:ndim)=0.d0
4488 ixb^l=ixo^l-kr(idims,^d);
4489 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4491 if({ ix^d==0 .and. ^d==idims | .or.})
then
4492 ixbmin^d=ixcmin^d-ix^d;
4493 ixbmax^d=ixcmax^d-ix^d;
4494 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4497 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4500 if(slab_uniform)
then
4502 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4503 ixb^l=ixo^l-kr(idims,^d);
4504 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4508 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4509 ixb^l=ixo^l-kr(idims,^d);
4510 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4512 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4514 end subroutine get_flux_on_cell_face
4518 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4521 integer,
intent(in) :: ixi^
l, ixo^
l
4522 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4523 double precision,
intent(in) :: w(ixi^s,1:nw)
4524 double precision :: dtnew
4526 double precision :: coef
4527 double precision :: dxarr(
ndim)
4528 double precision :: tmp(ixi^s)
4533 coef = maxval(abs(tmp(ixo^s)))
4540 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4542 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4545 end function get_ambipolar_dt
4553 integer,
intent(in) :: ixi^
l, ixo^
l
4554 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4555 double precision,
intent(inout) :: res(ixi^s)
4556 double precision :: tmp(ixi^s)
4557 double precision :: rho(ixi^s)
4566 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4570 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4577 integer,
intent(in) :: ixi^
l, ixo^
l
4578 double precision,
intent(in) :: qdt,dtfactor
4579 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4580 double precision,
intent(inout) :: w(ixi^s,1:nw)
4581 logical,
intent(in) :: qsourcesplit
4582 logical,
intent(inout) :: active
4589 if (.not. qsourcesplit)
then
4593 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4597 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4603 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4609 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4613 if (abs(
mhd_eta)>smalldouble)
then
4615 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4620 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4626 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4630 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4637 select case (type_divb)
4642 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4645 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4648 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4649 case (divb_janhunen)
4651 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4652 case (divb_lindejanhunen)
4654 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4655 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4656 case (divb_lindepowel)
4658 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4659 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4660 case (divb_lindeglm)
4662 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4663 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4664 case (divb_multigrid)
4669 call mpistop(
'Unknown divB fix')
4676 w,x,qsourcesplit,active,
rc_fl)
4686 w,x,gravity_energy,qsourcesplit,active)
4695 if(.not.qsourcesplit)
then
4697 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4701 end subroutine mhd_add_source
4703 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4707 integer,
intent(in) :: ixi^
l, ixo^
l
4708 double precision,
intent(in) :: qdt,dtfactor
4709 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4710 double precision,
intent(inout) :: w(ixi^s,1:nw)
4711 double precision :: divv(ixi^s)
4727 end subroutine add_pe0_divv
4729 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4731 integer,
intent(in) :: ixi^
l,ixo^
l
4732 double precision,
intent(in) :: qdt
4733 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4734 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4735 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4737 double precision :: r(ixi^s),te(ixi^s),rho_loc(ixi^s),pth_loc(ixi^s)
4738 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,bdir(
ndim),bunitvec(
ndim)
4742 {
do ix^db=iximin^db,iximax^db\}
4746 rho_loc(ix^
d)=wctprim(ix^
d,
rho_)
4751 pth_loc(ix^
d)=wctprim(ix^
d,
p_)
4753 te(ix^
d)=pth_loc(ix^
d)/(r(ix^
d)*rho_loc(ix^
d))
4759 do ix1=ixomin1,ixomax1
4761 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4763 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4766 sigma_t7=sigma_t5*te(ix^d)
4770 sigma_t7=sigma_t5*te(ix^d)
4772 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)
4774 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4775 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4776 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4778 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4779 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4784 do ix2=ixomin2,ixomax2
4785 do ix1=ixomin1,ixomax1
4787 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4789 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4792 sigma_t7=sigma_t5*te(ix^d)
4796 sigma_t7=sigma_t5*te(ix^d)
4799 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4801 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4803 if(bdir(1)/=0.d0)
then
4804 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
4808 if(bdir(2)/=0.d0)
then
4809 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
4813 sigmat5_bgradt=sigma_t5*(&
4814 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)&
4815 +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))
4817 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4818 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4819 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4821 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4822 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4828 do ix3=ixomin3,ixomax3
4829 do ix2=ixomin2,ixomax2
4830 do ix1=ixomin1,ixomax1
4832 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4834 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4837 sigma_t7=sigma_t5*te(ix^d)
4841 sigma_t7=sigma_t5*te(ix^d)
4844 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4846 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4848 if(bdir(1)/=0.d0)
then
4849 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
4853 if(bdir(2)/=0.d0)
then
4854 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
4858 if(bdir(3)/=0.d0)
then
4859 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
4863 sigmat5_bgradt=sigma_t5*(&
4864 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)&
4865 +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)&
4866 +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))
4868 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4869 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4870 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4872 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4873 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4879 end subroutine add_hypertc_source
4882 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4884 integer,
intent(in) :: ixi^
l, ixo^
l
4885 double precision,
intent(in) :: w(ixi^s,1:nw)
4886 double precision,
intent(inout) :: jxb(ixi^s,3)
4887 double precision :: a(ixi^s,3),
b(ixi^s,3)
4889 double precision :: current(ixi^s,7-2*
ndir:3)
4890 integer :: idir, idirmin
4895 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4899 b(ixo^s, idir) = w(ixo^s,mag(idir))
4908 a(ixo^s,idir)=current(ixo^s,idir)
4912 end subroutine get_lorentz_force
4916 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4918 integer,
intent(in) :: ixi^
l, ixo^
l
4919 double precision,
intent(in) :: w(ixi^s, nw)
4920 double precision,
intent(out) :: gamma_a2(ixo^s)
4921 double precision :: rho(ixi^s)
4927 rho(ixo^s) = w(ixo^s,
rho_)
4930 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4931 end subroutine mhd_gamma2_alfven
4935 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4937 integer,
intent(in) :: ixi^
l, ixo^
l
4938 double precision,
intent(in) :: w(ixi^s, nw)
4939 double precision :: gamma_a(ixo^s)
4941 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4942 gamma_a = sqrt(gamma_a)
4943 end function mhd_gamma_alfven
4947 integer,
intent(in) :: ixi^
l, ixo^
l
4948 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4949 double precision,
intent(out) :: rho(ixi^s)
4954 rho(ixo^s) = w(ixo^s,
rho_)
4960 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4963 integer,
intent(in) :: ixi^
l,ixo^
l, ie
4964 double precision,
intent(inout) :: w(ixi^s,1:nw)
4965 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4966 character(len=*),
intent(in) :: subname
4968 double precision :: rho(ixi^s)
4970 logical :: flag(ixi^s,1:nw)
4974 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4975 flag(ixo^s,ie)=.true.
4977 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4979 if(any(flag(ixo^s,ie)))
then
4983 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4986 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4992 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4995 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
5001 end subroutine mhd_handle_small_ei
5003 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
5007 integer,
intent(in) :: ixi^
l, ixo^
l
5008 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5009 double precision,
intent(inout) :: w(ixi^s,1:nw)
5011 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
5020 end subroutine mhd_update_temperature
5023 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
5026 integer,
intent(in) :: ixi^
l, ixo^
l
5027 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5028 double precision,
intent(inout) :: w(ixi^s,1:nw)
5030 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
5042 a(ixo^s,idir)=
block%J0(ixo^s,idir)
5047 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5050 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5056 if(total_energy)
then
5059 b(ixo^s,:)=wct(ixo^s,mag(:))
5068 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5071 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5075 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5079 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5084 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5089 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5091 end subroutine add_source_b0split
5094 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5098 integer,
intent(in) :: ixi^
l, ixo^
l
5099 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5100 double precision,
intent(inout) :: w(ixi^s,1:nw)
5101 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5103 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5104 integer :: idir, idirmin, ix^
d
5108 {
do ix^db=iximin^db,iximax^db\}
5110 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5111 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5112 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5114 call divvector(e,ixi^l,ixo^l,dive)
5116 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5119 {
do ix^db=ixomin^db,ixomax^db\}
5120 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5121 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5122 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5123 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5124 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5125 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5129 end subroutine add_source_semirelativistic
5132 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5136 integer,
intent(in) :: ixi^
l, ixo^
l
5137 double precision,
intent(in) :: qdt
5138 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5139 double precision,
intent(inout) :: w(ixi^s,1:nw)
5140 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5142 double precision :: divv(ixi^s), tmp
5154 {
do ix^db=ixomin^db,ixomax^db\}
5156 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5157 if(w(ix^
d,
e_)<small_e)
then
5162 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
5165 if(fix_small_values)
then
5166 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5168 end subroutine add_source_internal_e
5171 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5176 integer,
intent(in) :: ixi^
l, ixo^
l
5177 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5178 double precision,
intent(inout) :: w(ixi^s,1:nw)
5179 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5181 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5182 double precision :: current(ixi^s,7-2*
ndir:3)
5183 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5184 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5185 integer :: idir, idirmin, idims, ix^
d
5190 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5198 j(ixo^s,idir)=current(ixo^s,idir)
5216 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
5222 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5229 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)
5233 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5237 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5238 tmp(ixo^s)=sqrt(b2(ixo^s))
5239 where(tmp(ixo^s)>smalldouble)
5240 tmp(ixo^s)=1.d0/tmp(ixo^s)
5246 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5251 {
do ix^db=ixomin^db,ixomax^db\}
5253 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5255 b2(ix^
d)=vaoc/(1.d0+vaoc)
5258 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5261 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5265 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5268 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5269 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5272 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5275 end subroutine add_source_hydrodynamic_e
5281 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5286 integer,
intent(in) :: ixi^
l, ixo^
l
5287 double precision,
intent(in) :: qdt
5288 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5289 double precision,
intent(inout) :: w(ixi^s,1:nw)
5290 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5291 integer :: lxo^
l, kxo^
l
5293 double precision :: tmp(ixi^s),tmp2(ixi^s)
5296 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5297 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5306 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5307 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5314 gradeta(ixo^s,1:
ndim)=zero
5320 gradeta(ixo^s,idim)=tmp(ixo^s)
5327 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5334 tmp2(ixi^s)=bf(ixi^s,idir)
5336 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5337 jxo^
l=ixo^
l+
kr(idim,^
d);
5338 hxo^
l=ixo^
l-
kr(idim,^
d);
5339 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5340 tmp(ixo^s)=tmp(ixo^s)+&
5341 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5346 tmp2(ixi^s)=bf(ixi^s,idir)
5348 jxo^
l=ixo^
l+
kr(idim,^
d);
5349 hxo^
l=ixo^
l-
kr(idim,^
d);
5350 tmp(ixo^s)=tmp(ixo^s)+&
5351 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5356 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5360 do jdir=1,
ndim;
do kdir=idirmin,3
5361 if (
lvc(idir,jdir,kdir)/=0)
then
5362 if (
lvc(idir,jdir,kdir)==1)
then
5363 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5365 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5372 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5373 if(total_energy)
then
5374 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5380 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5383 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5385 end subroutine add_source_res1
5389 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5394 integer,
intent(in) :: ixi^
l, ixo^
l
5395 double precision,
intent(in) :: qdt
5396 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5397 double precision,
intent(inout) :: w(ixi^s,1:nw)
5400 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5401 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5402 integer :: ixa^
l,idir,idirmin,idirmin1
5406 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5407 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5417 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5422 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5431 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5434 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5439 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5441 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5443 if(total_energy)
then
5446 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5447 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5450 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5454 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5455 end subroutine add_source_res2
5459 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5463 integer,
intent(in) :: ixi^
l, ixo^
l
5464 double precision,
intent(in) :: qdt
5465 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5466 double precision,
intent(inout) :: w(ixi^s,1:nw)
5468 double precision :: current(ixi^s,7-2*
ndir:3)
5469 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5470 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5473 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5474 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5477 tmpvec(ixa^s,1:
ndir)=zero
5479 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5483 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5486 tmpvec(ixa^s,1:
ndir)=zero
5487 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5491 tmpvec2(ixa^s,1:
ndir)=zero
5492 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5495 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5498 if(total_energy)
then
5501 tmpvec2(ixa^s,1:
ndir)=zero
5502 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5503 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5504 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5505 end do;
end do;
end do
5507 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5508 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5511 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5513 end subroutine add_source_hyperres
5515 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5522 integer,
intent(in) :: ixi^
l, ixo^
l
5523 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5524 double precision,
intent(inout) :: w(ixi^s,1:nw)
5526 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5547 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5550 if(total_energy)
then
5559 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5568 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5572 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5574 end subroutine add_source_glm
5577 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5580 integer,
intent(in) :: ixi^
l, ixo^
l
5581 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5582 double precision,
intent(inout) :: w(ixi^s,1:nw)
5584 double precision :: divb(ixi^s), ba(1:
ndir)
5585 integer :: idir, ix^
d
5591 {
do ix^db=ixomin^db,ixomax^db\}
5596 if (total_energy)
then
5602 {
do ix^db=ixomin^db,ixomax^db\}
5604 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5606 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5607 if (total_energy)
then
5609 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5614 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5616 end subroutine add_source_powel
5618 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5623 integer,
intent(in) :: ixi^
l, ixo^
l
5624 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5625 double precision,
intent(inout) :: w(ixi^s,1:nw)
5627 double precision :: divb(ixi^s)
5628 integer :: idir, ix^
d
5633 {
do ix^db=ixomin^db,ixomax^db\}
5638 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5640 end subroutine add_source_janhunen
5642 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5647 integer,
intent(in) :: ixi^
l, ixo^
l
5648 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5649 double precision,
intent(inout) :: w(ixi^s,1:nw)
5651 double precision :: divb(ixi^s),graddivb(ixi^s)
5652 integer :: idim, idir, ixp^
l, i^
d, iside
5653 logical,
dimension(-1:1^D&) :: leveljump
5661 if(i^
d==0|.and.) cycle
5662 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5663 leveljump(i^
d)=.true.
5665 leveljump(i^
d)=.false.
5674 i^dd=kr(^dd,^d)*(2*iside-3);
5675 if (leveljump(i^dd))
then
5677 ixpmin^d=ixomin^d-i^d
5679 ixpmax^d=ixomax^d-i^d
5690 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5692 {
do i^db=ixpmin^db,ixpmax^db\}
5694 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5696 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5698 if (typedivbdiff==
'all' .and. total_energy)
then
5700 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5705 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5707 end subroutine add_source_linde
5714 integer,
intent(in) :: ixi^
l, ixo^
l
5715 double precision,
intent(in) :: w(ixi^s,1:nw)
5716 double precision :: divb(ixi^s), dsurface(ixi^s)
5718 double precision :: invb(ixo^s)
5719 integer :: ixa^
l,idims
5721 call get_divb(w,ixi^
l,ixo^
l,divb)
5723 where(invb(ixo^s)/=0.d0)
5724 invb(ixo^s)=1.d0/invb(ixo^s)
5727 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5729 ixamin^
d=ixomin^
d-1;
5730 ixamax^
d=ixomax^
d-1;
5731 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5733 ixa^
l=ixo^
l-
kr(idims,^
d);
5734 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5736 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5737 block%dvolume(ixo^s)/dsurface(ixo^s)
5748 integer,
intent(in) :: ixo^
l, ixi^
l
5749 double precision,
intent(in) :: w(ixi^s,1:nw)
5750 integer,
intent(out) :: idirmin
5753 double precision :: current(ixi^s,7-2*
ndir:3)
5754 integer :: idir, idirmin0
5760 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5761 block%J0(ixo^s,idirmin0:3)
5765 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5773 integer,
intent(in) :: ixi^
l, ixo^
l
5774 double precision,
intent(inout) :: dtnew
5775 double precision,
intent(in) ::
dx^
d
5776 double precision,
intent(in) :: w(ixi^s,1:nw)
5777 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5779 double precision :: dxarr(
ndim)
5780 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5781 integer :: idirmin,idim
5795 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5798 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5824 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5831 end subroutine mhd_get_dt
5834 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5839 integer,
intent(in) :: ixi^
l, ixo^
l
5840 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5841 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5843 double precision :: tmp,tmp1,invr,cot
5845 integer :: mr_,mphi_
5846 integer :: br_,bphi_
5849 br_=mag(1); bphi_=mag(1)-1+
phi_
5854 {
do ix^db=ixomin^db,ixomax^db\}
5857 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5862 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5867 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5868 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5869 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5870 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5871 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5873 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5874 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5875 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5878 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5883 {
do ix^db=ixomin^db,ixomax^db\}
5885 if(local_timestep)
then
5886 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5891 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5897 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5900 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5901 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5905 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5911 cot=1.d0/tan(x(ix^d,2))
5915 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5916 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5918 if(.not.stagger_grid)
then
5919 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5921 tmp=tmp+wprim(ix^d,
psi_)*cot
5923 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5928 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5929 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5930 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5932 if(.not.stagger_grid)
then
5933 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5935 tmp=tmp+wprim(ix^d,
psi_)*cot
5937 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5940 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5941 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5942 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5943 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5944 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5946 if(.not.stagger_grid)
then
5947 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5948 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5949 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5950 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5951 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5958 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5961 end subroutine mhd_add_source_geom
5964 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5969 integer,
intent(in) :: ixi^
l, ixo^
l
5970 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5971 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5973 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
5975 integer :: mr_,mphi_
5976 integer :: br_,bphi_
5979 br_=mag(1); bphi_=mag(1)-1+
phi_
5984 {
do ix^db=ixomin^db,ixomax^db\}
5987 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5998 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
5999 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
6000 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6005 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6011 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
6012 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
6013 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
6014 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6015 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6016 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
6018 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6019 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6020 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6023 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
6024 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
6029 {
do ix^db=ixomin^db,ixomax^db\}
6031 if(local_timestep)
then
6032 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
6038 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
6039 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
6040 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6044 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6051 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6057 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
6060 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
6061 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
6062 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
6066 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
6072 cot=1.d0/tan(x(ix^d,2))
6076 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
6077 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
6079 if(.not.stagger_grid)
then
6080 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6082 tmp=tmp+wprim(ix^d,
psi_)*cot
6084 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6090 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6091 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
6092 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6093 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
6095 if(.not.stagger_grid)
then
6096 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6098 tmp=tmp+wprim(ix^d,
psi_)*cot
6100 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6103 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6104 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6105 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6106 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
6107 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6108 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6109 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
6111 if(.not.stagger_grid)
then
6112 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6113 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6114 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6115 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6116 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6123 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6126 end subroutine mhd_add_source_geom_semirelati
6129 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6134 integer,
intent(in) :: ixi^
l, ixo^
l
6135 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6136 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6138 double precision :: tmp,tmp1,tmp2,invr,cot
6140 integer :: mr_,mphi_
6141 integer :: br_,bphi_
6144 br_=mag(1); bphi_=mag(1)-1+
phi_
6149 {
do ix^db=ixomin^db,ixomax^db\}
6152 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6157 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6162 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6163 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6164 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6165 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6166 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6168 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6169 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6170 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6173 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6178 {
do ix^db=ixomin^db,ixomax^db\}
6180 if(local_timestep)
then
6181 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6185 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6186 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6189 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6193 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6194 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6195 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6197 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6198 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6203 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6209 cot=1.d0/tan(x(ix^d,2))
6214 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6215 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6216 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6218 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6219 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6222 if(.not.stagger_grid)
then
6224 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6225 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6227 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6230 tmp=tmp+wprim(ix^d,
psi_)*cot
6232 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6238 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6239 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6240 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6241 +(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)
6243 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6244 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6245 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6248 if(.not.stagger_grid)
then
6250 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6251 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6253 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6256 tmp=tmp+wprim(ix^d,
psi_)*cot
6258 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6262 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6263 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6264 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6265 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6266 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6267 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6268 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6269 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6270 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6272 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6273 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6274 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6275 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6276 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6279 if(.not.stagger_grid)
then
6281 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6282 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6283 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6284 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6285 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6286 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6287 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6288 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6289 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6291 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6292 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6293 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6294 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6295 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6303 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6306 end subroutine mhd_add_source_geom_split
6311 integer,
intent(in) :: ixi^
l, ixo^
l
6312 double precision,
intent(in) :: w(ixi^s, nw)
6313 double precision :: mge(ixo^s)
6316 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6318 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6322 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6325 integer,
intent(in) :: ixi^
l, ixo^
l
6326 double precision,
intent(in) :: w(ixi^s,nw)
6327 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6328 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6330 double precision :: current(ixi^s,7-2*
ndir:3)
6331 double precision :: rho(ixi^s)
6332 integer :: idir, idirmin, ix^
d
6337 do idir = idirmin,
ndir
6338 {
do ix^db=ixomin^db,ixomax^db\}
6339 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6343 end subroutine mhd_getv_hall
6345 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6348 integer,
intent(in) :: ixi^
l, ixo^
l
6349 double precision,
intent(in) :: w(ixi^s,nw)
6350 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6351 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6354 double precision :: current(ixi^s,7-2*
ndir:3)
6355 integer :: idir, idirmin
6362 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6363 do idir = idirmin, 3
6367 end subroutine mhd_get_jambi
6369 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6372 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6373 double precision,
intent(in) :: qt
6374 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6375 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6378 double precision :: db(ixo^s), dpsi(ixo^s)
6382 {
do ix^db=ixomin^db,ixomax^db\}
6383 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6384 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6385 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6386 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6395 {
do ix^db=ixomin^db,ixomax^db\}
6396 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6397 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6398 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6399 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6400 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6402 if(total_energy)
then
6403 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6404 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6406 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6408 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6411 if(total_energy)
then
6412 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6413 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6418 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6420 end subroutine mhd_modify_wlr
6422 subroutine mhd_boundary_adjust(igrid,psb)
6424 integer,
intent(in) :: igrid
6427 integer :: ib, idims, iside, ixo^
l, i^
d
6436 i^
d=
kr(^
d,idims)*(2*iside-3);
6437 if (neighbor_type(i^
d,igrid)/=1) cycle
6438 ib=(idims-1)*2+iside
6456 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6461 end subroutine mhd_boundary_adjust
6463 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6466 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6467 double precision,
intent(inout) :: w(ixg^s,1:nw)
6468 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6470 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6471 integer :: ix^
d,ixf^
l
6484 do ix1=ixfmax1,ixfmin1,-1
6485 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6486 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6487 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6490 do ix1=ixfmax1,ixfmin1,-1
6491 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6492 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6493 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6494 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6495 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6496 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6497 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6511 do ix1=ixfmax1,ixfmin1,-1
6512 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6513 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6514 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6515 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6516 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6517 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6520 do ix1=ixfmax1,ixfmin1,-1
6521 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6522 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6523 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6524 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6525 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6526 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6527 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6528 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6529 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6530 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6531 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6532 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6533 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6534 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6535 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6536 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6537 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6538 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6552 do ix1=ixfmin1,ixfmax1
6553 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6554 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6555 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6558 do ix1=ixfmin1,ixfmax1
6559 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6560 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6561 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6562 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6563 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6564 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6565 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6579 do ix1=ixfmin1,ixfmax1
6580 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6581 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6582 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6583 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6584 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6585 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6588 do ix1=ixfmin1,ixfmax1
6589 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6590 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6591 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6592 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6593 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6594 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6595 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6596 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6597 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6598 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6599 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6600 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6601 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6602 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6603 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6604 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6605 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6606 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6620 do ix2=ixfmax2,ixfmin2,-1
6621 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6622 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6623 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6626 do ix2=ixfmax2,ixfmin2,-1
6627 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6628 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6629 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6630 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6631 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6632 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6633 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6647 do ix2=ixfmax2,ixfmin2,-1
6648 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6649 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6650 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6651 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6652 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6653 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6656 do ix2=ixfmax2,ixfmin2,-1
6657 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6658 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6659 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6660 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6661 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6662 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6663 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6664 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6665 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6666 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6667 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6668 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6669 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6670 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6671 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6672 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6673 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6674 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6688 do ix2=ixfmin2,ixfmax2
6689 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6690 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6691 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6694 do ix2=ixfmin2,ixfmax2
6695 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6696 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6697 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6698 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6699 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6700 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6701 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6715 do ix2=ixfmin2,ixfmax2
6716 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6717 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6718 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6719 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6720 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6721 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6724 do ix2=ixfmin2,ixfmax2
6725 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6726 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6727 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6728 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6729 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6730 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6731 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6732 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6733 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6734 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6735 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6736 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6737 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6738 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6739 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6740 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6741 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6742 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6759 do ix3=ixfmax3,ixfmin3,-1
6760 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6761 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6762 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6763 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6764 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6765 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6768 do ix3=ixfmax3,ixfmin3,-1
6769 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6770 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6771 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6772 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6773 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6774 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6775 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6776 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6777 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6778 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6779 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6780 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6781 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6782 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6783 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6784 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6785 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6786 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6801 do ix3=ixfmin3,ixfmax3
6802 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6803 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6804 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6805 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6806 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6807 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6810 do ix3=ixfmin3,ixfmax3
6811 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6812 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6813 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6814 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6815 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6816 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6817 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6818 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6819 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6820 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6821 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6822 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6823 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6824 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6825 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6826 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6827 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6828 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6834 call mpistop(
"Special boundary is not defined for this region")
6837 end subroutine fixdivb_boundary
6846 double precision,
intent(in) :: qdt
6847 double precision,
intent(in) :: qt
6848 logical,
intent(inout) :: active
6851 integer,
parameter :: max_its = 50
6852 double precision :: residual_it(max_its), max_divb
6853 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6854 double precision :: res
6855 double precision,
parameter :: max_residual = 1
d-3
6856 double precision,
parameter :: residual_reduction = 1
d-10
6857 integer :: iigrid, igrid
6858 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6861 mg%operator_type = mg_laplacian
6869 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6870 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6873 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6874 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6876 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6877 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6880 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6881 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6885 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6886 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6887 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6895 do iigrid = 1, igridstail
6896 igrid = igrids(iigrid);
6899 lvl =
mg%boxes(id)%lvl
6900 nc =
mg%box_size_lvl(lvl)
6906 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6908 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6909 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6914 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6917 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6918 if (
mype == 0) print *,
"iteration vs residual"
6921 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6922 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6923 if (residual_it(n) < residual_reduction * max_divb)
exit
6925 if (
mype == 0 .and. n > max_its)
then
6926 print *,
"divb_multigrid warning: not fully converged"
6927 print *,
"current amplitude of divb: ", residual_it(max_its)
6928 print *,
"multigrid smallest grid: ", &
6929 mg%domain_size_lvl(:,
mg%lowest_lvl)
6930 print *,
"note: smallest grid ideally has <= 8 cells"
6931 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6932 print *,
"note: dx/dy/dz should be similar"
6936 call mg_fas_vcycle(
mg, max_res=res)
6937 if (res < max_residual)
exit
6939 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6944 do iigrid = 1, igridstail
6945 igrid = igrids(iigrid);
6954 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6958 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6960 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
6962 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6975 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6976 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6979 if(total_energy)
then
6981 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6984 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6994 subroutine mhd_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
6998 integer,
intent(in) :: ixi^
l, ixo^
l
6999 double precision,
intent(in) :: qt,qdt
7001 double precision,
intent(in) :: wp(ixi^s,1:nw)
7002 type(state) :: sct, s
7003 type(ct_velocity) :: vcts
7004 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7005 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7007 double precision :: circ(ixi^s,1:
ndim)
7009 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7010 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
7011 integer :: idim1,idim2,idir,iwdim1,iwdim2
7013 associate(bfaces=>s%ws,x=>s%x)
7020 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7027 i1kr^
d=
kr(idim1,^
d);
7030 i2kr^
d=
kr(idim2,^
d);
7033 if (
lvc(idim1,idim2,idir)==1)
then
7035 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7037 {
do ix^db=ixcmin^db,ixcmax^db\}
7038 fe(ix^
d,idir)=quarter*&
7039 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
7040 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
7042 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
7047 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
7055 if(
associated(usr_set_electric_field)) &
7056 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7058 circ(ixi^s,1:ndim)=zero
7063 ixcmin^d=ixomin^d-kr(idim1,^d);
7065 ixa^l=ixc^l-kr(idim2,^d);
7068 if(lvc(idim1,idim2,idir)==1)
then
7070 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7073 else if(lvc(idim1,idim2,idir)==-1)
then
7075 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7081 {
do ix^db=ixcmin^db,ixcmax^db\}
7083 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7085 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7092 end subroutine mhd_update_faces_average
7095 subroutine mhd_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7100 integer,
intent(in) :: ixi^
l, ixo^
l
7101 double precision,
intent(in) :: qt, qdt
7103 double precision,
intent(in) :: wp(ixi^s,1:nw)
7104 type(state) :: sct, s
7105 type(ct_velocity) :: vcts
7106 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7107 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7109 double precision :: circ(ixi^s,1:
ndim)
7111 double precision :: ecc(ixi^s,
sdim:3)
7112 double precision :: ein(ixi^s,
sdim:3)
7114 double precision :: el(ixi^s),er(ixi^s)
7116 double precision :: elc,erc
7118 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7120 double precision :: jce(ixi^s,
sdim:3)
7122 double precision :: xs(ixgs^t,1:
ndim)
7123 double precision :: gradi(ixgs^t)
7124 integer :: ixc^
l,ixa^
l
7125 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7127 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7130 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7136 {
do ix^db=iximin^db,iximax^db\}
7139 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_)
7140 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_)
7141 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_)
7144 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7151 {
do ix^db=iximin^db,iximax^db\}
7154 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7155 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7156 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7159 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7173 i1kr^d=kr(idim1,^d);
7176 i2kr^d=kr(idim2,^d);
7179 if (lvc(idim1,idim2,idir)==1)
then
7181 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7184 {
do ix^db=ixcmin^db,ixcmax^db\}
7185 fe(ix^d,idir)=quarter*&
7186 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7187 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7192 ixamax^d=ixcmax^d+i1kr^d;
7193 {
do ix^db=ixamin^db,ixamax^db\}
7194 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7195 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7198 do ix^db=ixcmin^db,ixcmax^db\}
7199 if(vnorm(ix^d,idim1)>0.d0)
then
7201 else if(vnorm(ix^d,idim1)<0.d0)
then
7202 elc=el({ix^d+i1kr^d})
7204 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7206 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7208 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7209 erc=er({ix^d+i1kr^d})
7211 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7213 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7218 ixamax^d=ixcmax^d+i2kr^d;
7219 {
do ix^db=ixamin^db,ixamax^db\}
7220 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7221 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7224 do ix^db=ixcmin^db,ixcmax^db\}
7225 if(vnorm(ix^d,idim2)>0.d0)
then
7227 else if(vnorm(ix^d,idim2)<0.d0)
then
7228 elc=el({ix^d+i2kr^d})
7230 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7232 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7234 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7235 erc=er({ix^d+i2kr^d})
7237 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7239 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7243 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7248 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7262 if (lvc(idim1,idim2,idir)==0) cycle
7264 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7265 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7268 xs(ixa^s,:)=x(ixa^s,:)
7269 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7270 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7271 if (lvc(idim1,idim2,idir)==1)
then
7272 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7274 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7281 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7283 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7287 {
do ix^db=ixomin^db,ixomax^db\}
7288 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7289 +ein(ix1,ix2-1,ix3-1,idir))
7290 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7291 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7293 else if(idir==2)
then
7294 {
do ix^db=ixomin^db,ixomax^db\}
7295 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7296 +ein(ix1-1,ix2,ix3-1,idir))
7297 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7298 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7301 {
do ix^db=ixomin^db,ixomax^db\}
7302 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7303 +ein(ix1-1,ix2-1,ix3,idir))
7304 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7305 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7311 {
do ix^db=ixomin^db,ixomax^db\}
7312 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7313 +ein(ix1-1,ix2-1,idir))
7314 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7315 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7320 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7326 if(
associated(usr_set_electric_field)) &
7327 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7329 circ(ixi^s,1:ndim)=zero
7334 ixcmin^d=ixomin^d-kr(idim1,^d);
7336 ixa^l=ixc^l-kr(idim2,^d);
7339 if(lvc(idim1,idim2,idir)==1)
then
7341 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7344 else if(lvc(idim1,idim2,idir)==-1)
then
7346 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7352 {
do ix^db=ixcmin^db,ixcmax^db\}
7354 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7356 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7363 end subroutine mhd_update_faces_contact
7366 subroutine mhd_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7371 integer,
intent(in) :: ixi^
l, ixo^
l
7372 double precision,
intent(in) :: qt, qdt
7374 double precision,
intent(in) :: wp(ixi^s,1:nw)
7375 type(state) :: sct, s
7376 type(ct_velocity) :: vcts
7377 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7378 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7380 double precision :: vtill(ixi^s,2)
7381 double precision :: vtilr(ixi^s,2)
7382 double precision :: bfacetot(ixi^s,
ndim)
7383 double precision :: btill(ixi^s,
ndim)
7384 double precision :: btilr(ixi^s,
ndim)
7385 double precision :: cp(ixi^s,2)
7386 double precision :: cm(ixi^s,2)
7387 double precision :: circ(ixi^s,1:
ndim)
7389 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7390 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7391 integer :: idim1,idim2,idir,ix^
d
7393 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7394 cbarmax=>vcts%cbarmax)
7407 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7423 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7427 idim2=mod(idir+1,3)+1
7429 jxc^
l=ixc^
l+
kr(idim1,^
d);
7430 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7434 vtill(ixi^s,2),vtilr(ixi^s,2))
7437 vtill(ixi^s,1),vtilr(ixi^s,1))
7443 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7444 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7446 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7447 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7450 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7453 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7457 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7458 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7460 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7461 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7465 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7466 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7467 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7468 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7469 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7470 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7471 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7472 /(cp(ixc^s,2)+cm(ixc^s,2))
7475 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7479 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7493 circ(ixi^s,1:
ndim)=zero
7498 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7502 if(
lvc(idim1,idim2,idir)/=0)
then
7503 hxc^
l=ixc^
l-
kr(idim2,^
d);
7505 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7506 +
lvc(idim1,idim2,idir)&
7512 {
do ix^db=ixcmin^db,ixcmax^db\}
7514 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7516 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7522 end subroutine mhd_update_faces_hll
7525 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
7530 integer,
intent(in) :: ixi^
l, ixo^
l
7532 double precision,
intent(in) :: wp(ixi^s,1:nw)
7533 type(state),
intent(in) :: sct, s
7535 double precision :: jce(ixi^s,
sdim:3)
7538 double precision :: jcc(ixi^s,7-2*
ndir:3)
7540 double precision :: xs(ixgs^t,1:
ndim)
7542 double precision :: eta(ixi^s)
7543 double precision :: gradi(ixgs^t)
7544 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7546 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7552 if (
lvc(idim1,idim2,idir)==0) cycle
7554 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7555 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7558 xs(ixb^s,:)=x(ixb^s,:)
7559 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7560 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7561 if (
lvc(idim1,idim2,idir)==1)
then
7562 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7564 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7571 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7579 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7580 jcc(ixc^s,idir)=0.d0
7582 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7583 ixamin^
d=ixcmin^
d+ix^
d;
7584 ixamax^
d=ixcmax^
d+ix^
d;
7585 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7587 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7588 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7593 end subroutine get_resistive_electric_field
7596 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7599 integer,
intent(in) :: ixi^
l, ixo^
l
7600 double precision,
intent(in) :: w(ixi^s,1:nw)
7601 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7602 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7604 double precision :: jxbxb(ixi^s,1:3)
7605 integer :: idir,ixa^
l,ixc^
l,ix^
d
7608 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7615 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7618 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7619 ixamin^
d=ixcmin^
d+ix^
d;
7620 ixamax^
d=ixcmax^
d+ix^
d;
7621 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7623 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7626 end subroutine get_ambipolar_electric_field
7632 integer,
intent(in) :: ixo^
l
7642 do ix^db=ixomin^db,ixomax^db\}
7644 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7645 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7646 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7647 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7648 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7649 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7652 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7653 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7654 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7655 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7698 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7699 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7700 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7702 double precision :: adummy(ixis^s,1:3)
7708 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7711 integer,
intent(in) :: ixi^
l, ixo^
l
7712 double precision,
intent(in) :: w(ixi^s,1:nw)
7713 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7714 double precision,
intent(out):: rfactor(ixi^s)
7716 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7720 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)
7722 end subroutine rfactor_from_temperature_ionization
7724 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7726 integer,
intent(in) :: ixi^
l, ixo^
l
7727 double precision,
intent(in) :: w(ixi^s,1:nw)
7728 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7729 double precision,
intent(out):: rfactor(ixi^s)
7733 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(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.