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.
167 logical :: compactres = .false.
181 logical :: total_energy = .true.
185 logical :: gravity_energy
187 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
189 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
191 character(len=std_len) :: typedivbdiff =
'all'
202 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
204 integer,
intent(in) :: ixi^
l, ixo^
l
205 double precision,
intent(in) :: x(ixi^s,1:
ndim)
206 double precision,
intent(in) :: w(ixi^s,1:nw)
207 double precision,
intent(inout) :: res(ixi^s)
208 end subroutine mask_subroutine
246 subroutine mhd_read_params(files)
249 character(len=*),
intent(in) :: files(:)
265 do n = 1,
size(files)
266 open(
unitpar, file=trim(files(n)), status=
"old")
267 read(
unitpar, mhd_list,
end=111)
271 end subroutine mhd_read_params
274 subroutine mhd_write_info(fh)
276 integer,
intent(in) :: fh
279 integer,
parameter :: n_par = 1
280 double precision :: values(n_par)
281 integer,
dimension(MPI_STATUS_SIZE) :: st
282 character(len=name_len) :: names(n_par)
284 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
288 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
289 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
290 end subroutine mhd_write_info
317 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
321 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_internal_e=T'
327 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_hydrodynamic_e=T'
331 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_hydrodynamic_e=T'
335 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_hydrodynamic_e=T'
346 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
350 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
354 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
358 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
362 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
366 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
370 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
374 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
378 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho_and_p=F when mhd_energy=F'
384 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
390 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
394 if(
mype==0)
write(*,*)
'WARNING: turn off hyperbolic TC when using parabolic TC'
418 phys_total_energy=total_energy
421 gravity_energy=.false.
423 gravity_energy=.true.
426 gravity_energy=.false.
432 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
437 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
445 type_divb = divb_none
448 type_divb = divb_multigrid
450 mg%operator_type = mg_laplacian
457 case (
'powel',
'powell')
458 type_divb = divb_powel
460 type_divb = divb_janhunen
462 type_divb = divb_linde
463 case (
'lindejanhunen')
464 type_divb = divb_lindejanhunen
466 type_divb = divb_lindepowel
470 type_divb = divb_lindeglm
475 call mpistop(
'Unknown divB fix')
480 allocate(start_indices(number_species),stop_indices(number_species))
487 mom(:) = var_set_momentum(
ndir)
493 e_ = var_set_energy()
502 mag(:) = var_set_bfield(
ndir)
506 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
522 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
527 te_ = var_set_auxvar(
'Te',
'Te')
536 stop_indices(1)=nwflux
565 allocate(iw_vector(nvector))
566 iw_vector(1) =
mom(1) - 1
567 iw_vector(2) = mag(1) - 1
570 if (.not.
allocated(flux_type))
then
571 allocate(flux_type(
ndir, nwflux))
572 flux_type = flux_default
573 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
574 call mpistop(
"phys_check error: flux_type has wrong shape")
577 if(nwflux>mag(
ndir))
then
579 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
584 flux_type(:,
psi_)=flux_special
586 flux_type(idir,mag(idir))=flux_special
590 flux_type(idir,mag(idir))=flux_tvdlf
596 phys_get_dt => mhd_get_dt
599 phys_get_cmax => mhd_get_cmax_semirelati
601 phys_get_cmax => mhd_get_cmax_semirelati_noe
605 phys_get_cmax => mhd_get_cmax_origin
607 phys_get_cmax => mhd_get_cmax_origin_noe
610 phys_get_a2max => mhd_get_a2max
611 phys_get_tcutoff => mhd_get_tcutoff
612 phys_get_h_speed => mhd_get_h_speed
614 phys_get_cbounds => mhd_get_cbounds_split_rho
616 phys_get_cbounds => mhd_get_cbounds_semirelati
618 phys_get_cbounds => mhd_get_cbounds
621 phys_to_primitive => mhd_to_primitive_hde
623 phys_to_conserved => mhd_to_conserved_hde
627 phys_to_primitive => mhd_to_primitive_semirelati
629 phys_to_conserved => mhd_to_conserved_semirelati
632 phys_to_primitive => mhd_to_primitive_semirelati_noe
634 phys_to_conserved => mhd_to_conserved_semirelati_noe
639 phys_to_primitive => mhd_to_primitive_split_rho
641 phys_to_conserved => mhd_to_conserved_split_rho
644 phys_to_primitive => mhd_to_primitive_inte
646 phys_to_conserved => mhd_to_conserved_inte
649 phys_to_primitive => mhd_to_primitive_origin
651 phys_to_conserved => mhd_to_conserved_origin
654 phys_to_primitive => mhd_to_primitive_origin_noe
656 phys_to_conserved => mhd_to_conserved_origin_noe
661 phys_get_flux => mhd_get_flux_hde
664 phys_get_flux => mhd_get_flux_semirelati
666 phys_get_flux => mhd_get_flux_semirelati_noe
670 phys_get_flux => mhd_get_flux_split
672 phys_get_flux => mhd_get_flux
674 phys_get_flux => mhd_get_flux_noe
679 phys_add_source_geom => mhd_add_source_geom_semirelati
681 phys_add_source_geom => mhd_add_source_geom_split
683 phys_add_source_geom => mhd_add_source_geom
685 phys_add_source => mhd_add_source
686 phys_check_params => mhd_check_params
687 phys_write_info => mhd_write_info
690 phys_handle_small_values => mhd_handle_small_values_inte
691 mhd_handle_small_values => mhd_handle_small_values_inte
692 phys_check_w => mhd_check_w_inte
694 phys_handle_small_values => mhd_handle_small_values_hde
695 mhd_handle_small_values => mhd_handle_small_values_hde
696 phys_check_w => mhd_check_w_hde
698 phys_handle_small_values => mhd_handle_small_values_semirelati
699 mhd_handle_small_values => mhd_handle_small_values_semirelati
700 phys_check_w => mhd_check_w_semirelati
702 phys_handle_small_values => mhd_handle_small_values_split
703 mhd_handle_small_values => mhd_handle_small_values_split
704 phys_check_w => mhd_check_w_split
706 phys_handle_small_values => mhd_handle_small_values_origin
707 mhd_handle_small_values => mhd_handle_small_values_origin
708 phys_check_w => mhd_check_w_origin
710 phys_handle_small_values => mhd_handle_small_values_noe
711 mhd_handle_small_values => mhd_handle_small_values_noe
712 phys_check_w => mhd_check_w_noe
716 phys_get_pthermal => mhd_get_pthermal_inte
719 phys_get_pthermal => mhd_get_pthermal_hde
722 phys_get_pthermal => mhd_get_pthermal_semirelati
725 phys_get_pthermal => mhd_get_pthermal_origin
728 phys_get_pthermal => mhd_get_pthermal_noe
733 phys_set_equi_vars => set_equi_vars_grid
736 if(type_divb==divb_glm)
then
737 phys_modify_wlr => mhd_modify_wlr
743 phys_update_temperature => mhd_update_temperature
768 transverse_ghost_cells = 1
769 phys_get_ct_velocity => mhd_get_ct_velocity_average
770 phys_update_faces => mhd_update_faces_average
772 transverse_ghost_cells = 1
773 phys_get_ct_velocity => mhd_get_ct_velocity_contact
774 phys_update_faces => mhd_update_faces_contact
776 transverse_ghost_cells = 2
777 phys_get_ct_velocity => mhd_get_ct_velocity_hll
778 phys_update_faces => mhd_update_faces_hll
780 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
783 phys_modify_wlr => mhd_modify_wlr
785 phys_boundary_adjust => mhd_boundary_adjust
794 call mhd_physical_units()
812 if(phys_internal_e)
then
814 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
816 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
819 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
822 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
824 tc_fl%has_equi = .true.
825 tc_fl%get_temperature_equi => mhd_get_temperature_equi
826 tc_fl%get_rho_equi => mhd_get_rho_equi
828 tc_fl%has_equi = .false.
831 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
859 rc_fl%has_equi = .true.
860 rc_fl%get_rho_equi => mhd_get_rho_equi
861 rc_fl%get_pthermal_equi => mhd_get_pe_equi
863 rc_fl%has_equi = .false.
871 phys_te_images => mhd_te_images
887 if (particles_eta < zero) particles_eta =
mhd_eta
888 if (particles_etah < zero) particles_eta =
mhd_etah
890 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
891 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
905 phys_wider_stencil = 2
907 phys_wider_stencil = 1
915 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
927 phys_wider_stencil = 2
929 phys_wider_stencil = 1
943 subroutine mhd_te_images
948 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
950 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
952 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
954 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
957 call mpistop(
"Error in synthesize emission: Unknown convert_type")
959 end subroutine mhd_te_images
965 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
969 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
970 double precision,
intent(in) :: x(ixi^s,1:
ndim)
971 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
972 double precision,
intent(in) :: my_dt
973 logical,
intent(in) :: fix_conserve_at_step
975 end subroutine mhd_sts_set_source_tc_mhd
977 subroutine mhd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
981 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
982 double precision,
intent(in) :: x(ixi^s,1:
ndim)
983 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
984 double precision,
intent(in) :: my_dt
985 logical,
intent(in) :: fix_conserve_at_step
987 end subroutine mhd_sts_set_source_tc_hd
989 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
996 integer,
intent(in) :: ixi^
l, ixo^
l
997 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
998 double precision,
intent(in) :: w(ixi^s,1:nw)
999 double precision :: dtnew
1002 end function mhd_get_tc_dt_mhd
1004 function mhd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
1011 integer,
intent(in) :: ixi^
l, ixo^
l
1012 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
1013 double precision,
intent(in) :: w(ixi^s,1:nw)
1014 double precision :: dtnew
1017 end function mhd_get_tc_dt_hd
1019 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
1022 integer,
intent(in) :: ixi^
l,ixo^
l
1023 double precision,
intent(inout) :: w(ixi^s,1:nw)
1024 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1025 integer,
intent(in) :: step
1026 character(len=140) :: error_msg
1028 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1029 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1030 end subroutine mhd_tc_handle_small_e
1033 subroutine tc_params_read_mhd(fl)
1035 type(tc_fluid),
intent(inout) :: fl
1037 double precision :: tc_k_para=0d0
1038 double precision :: tc_k_perp=0d0
1041 logical :: tc_perpendicular=.false.
1042 logical :: tc_saturate=.false.
1043 character(len=std_len) :: tc_slope_limiter=
"MC"
1045 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1049 read(
unitpar, tc_list,
end=111)
1053 fl%tc_perpendicular = tc_perpendicular
1054 fl%tc_saturate = tc_saturate
1055 fl%tc_k_para = tc_k_para
1056 fl%tc_k_perp = tc_k_perp
1057 select case(tc_slope_limiter)
1059 fl%tc_slope_limiter = 0
1062 fl%tc_slope_limiter = 1
1065 fl%tc_slope_limiter = 2
1068 fl%tc_slope_limiter = 3
1071 fl%tc_slope_limiter = 4
1073 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod, superbee, koren")
1075 end subroutine tc_params_read_mhd
1079 subroutine rc_params_read(fl)
1082 type(rc_fluid),
intent(inout) :: fl
1084 double precision :: cfrac=0.1d0
1087 double precision :: rad_cut_hgt=0.5d0
1088 double precision :: rad_cut_dey=0.15d0
1091 integer :: ncool = 4000
1093 logical :: tfix=.false.
1095 logical :: rc_split=.false.
1096 logical :: rad_cut=.false.
1098 character(len=std_len) :: coolcurve=
'JCcorona'
1100 character(len=std_len) :: coolmethod=
'exact'
1102 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1106 read(
unitpar, rc_list,
end=111)
1111 fl%coolcurve=coolcurve
1112 fl%coolmethod=coolmethod
1115 fl%rc_split=rc_split
1118 fl%rad_cut_hgt=rad_cut_hgt
1119 fl%rad_cut_dey=rad_cut_dey
1120 end subroutine rc_params_read
1124 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1127 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1128 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1130 double precision :: delx(ixi^s,1:
ndim)
1131 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1132 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1138 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1142 hxo^
l=ixo^
l-
kr(idims,^
d);
1148 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1151 xshift^
d=half*(one-
kr(^
d,idims));
1158 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1162 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1165 end subroutine set_equi_vars_grid_faces
1168 subroutine set_equi_vars_grid(igrid)
1172 integer,
intent(in) :: igrid
1178 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1180 end subroutine set_equi_vars_grid
1183 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1185 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1186 double precision,
intent(in) :: w(ixi^s, 1:nw)
1187 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1188 double precision :: wnew(ixo^s, 1:nwc)
1195 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1201 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1205 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1209 if(
b0field .and. total_energy)
then
1210 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1211 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1215 end function convert_vars_splitting
1217 subroutine mhd_check_params
1225 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1226 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1230 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1231 inv_gamma_1=1.d0/gamma_1
1236 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1241 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1246 end subroutine mhd_check_params
1248 subroutine mhd_physical_units()
1250 double precision :: mp,kb,miu0,c_lightspeed
1251 double precision :: a,
b
1262 c_lightspeed=const_c
1403 end subroutine mhd_physical_units
1405 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1408 logical,
intent(in) :: primitive
1409 logical,
intent(inout) :: flag(ixi^s,1:nw)
1410 integer,
intent(in) :: ixi^
l, ixo^
l
1411 double precision,
intent(in) :: w(ixi^s,nw)
1413 double precision :: tmp,b2,
b(ixo^s,1:
ndir)
1414 double precision :: v(ixo^s,1:
ndir),gamma2,inv_rho
1425 {
do ix^db=ixomin^db,ixomax^db \}
1426 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1429 {
do ix^db=ixomin^db,ixomax^db \}
1430 b2=(^
c&w(ix^d,
b^
c_)**2+)
1431 if(b2>smalldouble)
then
1436 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1437 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1438 inv_rho = 1d0/w(ix^d,
rho_)
1440 b2=b2*inv_rho*inv_squared_c
1442 gamma2=1.d0/(1.d0+b2)
1444 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1447 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1448 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1449 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1454 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1460 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1461 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1462 if(tmp<small_e) flag(ix^d,
e_)=.true.
1468 end subroutine mhd_check_w_semirelati
1470 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1473 logical,
intent(in) :: primitive
1474 integer,
intent(in) :: ixi^
l, ixo^
l
1475 double precision,
intent(in) :: w(ixi^s,nw)
1476 logical,
intent(inout) :: flag(ixi^s,1:nw)
1481 {
do ix^db=ixomin^db,ixomax^db\}
1487 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1491 end subroutine mhd_check_w_origin
1493 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1496 logical,
intent(in) :: primitive
1497 integer,
intent(in) :: ixi^
l, ixo^
l
1498 double precision,
intent(in) :: w(ixi^s,nw)
1499 logical,
intent(inout) :: flag(ixi^s,1:nw)
1501 double precision :: tmp
1505 {
do ix^db=ixomin^db,ixomax^db\}
1511 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1516 end subroutine mhd_check_w_split
1518 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1521 logical,
intent(in) :: primitive
1522 integer,
intent(in) :: ixi^
l, ixo^
l
1523 double precision,
intent(in) :: w(ixi^s,nw)
1524 logical,
intent(inout) :: flag(ixi^s,1:nw)
1529 {
do ix^db=ixomin^db,ixomax^db\}
1533 end subroutine mhd_check_w_noe
1535 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1538 logical,
intent(in) :: primitive
1539 integer,
intent(in) :: ixi^
l, ixo^
l
1540 double precision,
intent(in) :: w(ixi^s,nw)
1541 logical,
intent(inout) :: flag(ixi^s,1:nw)
1546 {
do ix^db=ixomin^db,ixomax^db\}
1551 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1555 end subroutine mhd_check_w_inte
1557 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1560 logical,
intent(in) :: primitive
1561 integer,
intent(in) :: ixi^
l, ixo^
l
1562 double precision,
intent(in) :: w(ixi^s,nw)
1563 logical,
intent(inout) :: flag(ixi^s,1:nw)
1568 {
do ix^db=ixomin^db,ixomax^db\}
1573 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1577 end subroutine mhd_check_w_hde
1580 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1582 integer,
intent(in) :: ixi^
l, ixo^
l
1583 double precision,
intent(inout) :: w(ixi^s, nw)
1584 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1588 {
do ix^db=ixomin^db,ixomax^db\}
1590 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1592 +(^
c&w(ix^
d,
b^
c_)**2+))
1597 end subroutine mhd_to_conserved_origin
1600 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1602 integer,
intent(in) :: ixi^
l, ixo^
l
1603 double precision,
intent(inout) :: w(ixi^s, nw)
1604 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1608 {
do ix^db=ixomin^db,ixomax^db\}
1613 end subroutine mhd_to_conserved_origin_noe
1616 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1618 integer,
intent(in) :: ixi^
l, ixo^
l
1619 double precision,
intent(inout) :: w(ixi^s, nw)
1620 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1624 {
do ix^db=ixomin^db,ixomax^db\}
1626 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1632 end subroutine mhd_to_conserved_hde
1635 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1637 integer,
intent(in) :: ixi^
l, ixo^
l
1638 double precision,
intent(inout) :: w(ixi^s, nw)
1639 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1643 {
do ix^db=ixomin^db,ixomax^db\}
1645 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1650 end subroutine mhd_to_conserved_inte
1653 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1655 integer,
intent(in) :: ixi^
l, ixo^
l
1656 double precision,
intent(inout) :: w(ixi^s, nw)
1657 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1659 double precision :: rho
1662 {
do ix^db=ixomin^db,ixomax^db\}
1665 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1666 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1667 +(^
c&w(ix^
d,
b^
c_)**2+))
1672 end subroutine mhd_to_conserved_split_rho
1675 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1677 integer,
intent(in) :: ixi^
l, ixo^
l
1678 double precision,
intent(inout) :: w(ixi^s, nw)
1679 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1681 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1684 {
do ix^db=ixomin^db,ixomax^db\}
1686 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1687 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1688 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1689 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1690 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1691 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1696 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1697 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1698 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1706 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1710 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1712 +(^
c&w(ix^
d,
b^
c_)**2+)&
1713 +(^
c&e(ix^
d,^
c)**2+)*inv_squared_c)
1721 end subroutine mhd_to_conserved_semirelati
1723 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1725 integer,
intent(in) :: ixi^
l, ixo^
l
1726 double precision,
intent(inout) :: w(ixi^s, nw)
1727 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1729 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1732 {
do ix^db=ixomin^db,ixomax^db\}
1734 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1735 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1736 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1737 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1738 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1739 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1744 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1745 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1746 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1756 end subroutine mhd_to_conserved_semirelati_noe
1759 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1761 integer,
intent(in) :: ixi^
l, ixo^
l
1762 double precision,
intent(inout) :: w(ixi^s, nw)
1763 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1765 double precision :: inv_rho
1770 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1773 {
do ix^db=ixomin^db,ixomax^db\}
1774 inv_rho = 1.d0/w(ix^
d,
rho_)
1778 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1780 +(^
c&w(ix^
d,
b^
c_)**2+)))
1783 end subroutine mhd_to_primitive_origin
1786 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1788 integer,
intent(in) :: ixi^
l, ixo^
l
1789 double precision,
intent(inout) :: w(ixi^s, nw)
1790 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1792 double precision :: inv_rho
1797 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1800 {
do ix^db=ixomin^db,ixomax^db\}
1801 inv_rho = 1.d0/w(ix^
d,
rho_)
1806 end subroutine mhd_to_primitive_origin_noe
1809 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1811 integer,
intent(in) :: ixi^
l, ixo^
l
1812 double precision,
intent(inout) :: w(ixi^s, nw)
1813 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1815 double precision :: inv_rho
1820 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1823 {
do ix^db=ixomin^db,ixomax^db\}
1824 inv_rho = 1d0/w(ix^
d,
rho_)
1831 end subroutine mhd_to_primitive_hde
1834 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1836 integer,
intent(in) :: ixi^
l, ixo^
l
1837 double precision,
intent(inout) :: w(ixi^s, nw)
1838 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1840 double precision :: inv_rho
1845 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1848 {
do ix^db=ixomin^db,ixomax^db\}
1850 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1852 inv_rho = 1.d0/w(ix^
d,
rho_)
1856 end subroutine mhd_to_primitive_inte
1859 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1861 integer,
intent(in) :: ixi^
l, ixo^
l
1862 double precision,
intent(inout) :: w(ixi^s, nw)
1863 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1865 double precision :: inv_rho
1870 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1873 {
do ix^db=ixomin^db,ixomax^db\}
1878 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1880 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1883 end subroutine mhd_to_primitive_split_rho
1886 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1888 integer,
intent(in) :: ixi^
l, ixo^
l
1889 double precision,
intent(inout) :: w(ixi^s, nw)
1890 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1892 double precision ::
b(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
1897 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1900 {
do ix^db=ixomin^db,ixomax^db\}
1901 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1902 if(b2>smalldouble)
then
1910 inv_rho=1.d0/w(ix^
d,
rho_)
1912 b2=b2*inv_rho*inv_squared_c
1914 gamma2=1.d0/(1.d0+b2)
1916 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1920 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1924 b(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1925 b(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1926 b(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1930 b(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1936 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1938 +(^
c&w(ix^
d,
b^
c_)**2+)&
1939 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
1943 end subroutine mhd_to_primitive_semirelati
1946 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1948 integer,
intent(in) :: ixi^
l, ixo^
l
1949 double precision,
intent(inout) :: w(ixi^s, nw)
1950 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1952 double precision ::
b(ixo^s,1:
ndir),tmp,b2,gamma2,inv_rho
1953 integer :: ix^
d, idir
1957 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1960 {
do ix^db=ixomin^db,ixomax^db\}
1961 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1962 if(b2>smalldouble)
then
1970 inv_rho=1.d0/w(ix^
d,
rho_)
1972 b2=b2*inv_rho*inv_squared_c
1974 gamma2=1.d0/(1.d0+b2)
1976 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1979 end subroutine mhd_to_primitive_semirelati_noe
1984 integer,
intent(in) :: ixi^
l, ixo^
l
1985 double precision,
intent(inout) :: w(ixi^s, nw)
1986 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1991 {
do ix^db=ixomin^db,ixomax^db\}
1994 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1996 +(^
c&w(ix^
d,
b^
c_)**2+))
1999 {
do ix^db=ixomin^db,ixomax^db\}
2001 w(ix^d,
e_)=w(ix^d,
e_)&
2002 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2003 +(^
c&w(ix^d,
b^
c_)**2+))
2010 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
2012 integer,
intent(in) :: ixi^
l, ixo^
l
2013 double precision,
intent(inout) :: w(ixi^s, nw)
2014 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2018 {
do ix^db=ixomin^db,ixomax^db\}
2024 end subroutine mhd_ei_to_e_hde
2027 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
2029 integer,
intent(in) :: ixi^
l, ixo^
l
2030 double precision,
intent(inout) :: w(ixi^s, nw)
2031 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2033 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2034 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2036 end subroutine mhd_ei_to_e_semirelati
2041 integer,
intent(in) :: ixi^
l, ixo^
l
2042 double precision,
intent(inout) :: w(ixi^s, nw)
2043 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2048 {
do ix^db=ixomin^db,ixomax^db\}
2051 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2053 +(^
c&w(ix^
d,
b^
c_)**2+))
2056 {
do ix^db=ixomin^db,ixomax^db\}
2058 w(ix^d,
e_)=w(ix^d,
e_)&
2059 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2060 +(^
c&w(ix^d,
b^
c_)**2+))
2064 if(fix_small_values)
then
2065 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2071 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2073 integer,
intent(in) :: ixi^
l, ixo^
l
2074 double precision,
intent(inout) :: w(ixi^s, nw)
2075 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2079 {
do ix^db=ixomin^db,ixomax^db\}
2085 if(fix_small_values)
then
2086 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2089 end subroutine mhd_e_to_ei_hde
2092 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2094 integer,
intent(in) :: ixi^
l, ixo^
l
2095 double precision,
intent(inout) :: w(ixi^s, nw)
2096 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2098 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2099 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2101 end subroutine mhd_e_to_ei_semirelati
2103 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2106 logical,
intent(in) :: primitive
2107 integer,
intent(in) :: ixi^
l,ixo^
l
2108 double precision,
intent(inout) :: w(ixi^s,1:nw)
2109 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2110 character(len=*),
intent(in) :: subname
2112 double precision ::
b(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2113 double precision :: tmp, b2, gamma2, inv_rho
2115 logical :: flag(ixi^s,1:nw)
2124 {
do ix^db=ixomin^db,ixomax^db\}
2125 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2126 if(b2>smalldouble)
then
2133 inv_rho=1.d0/w(ix^
d,
rho_)
2135 b2=b2*inv_rho*inv_squared_c
2137 gamma2=1.d0/(1.d0+b2)
2139 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
2142 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2143 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2144 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2148 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2154 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2155 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2156 +(^
c&w(ix^
d,
b^
c_)**2+)&
2157 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
2164 select case (small_values_method)
2166 {
do ix^db=ixomin^db,ixomax^db\}
2167 if(flag(ix^d,
rho_))
then
2168 w(ix^d,
rho_) = small_density
2169 ^
c&w(ix^d,
m^
c_)=0.d0\
2173 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2175 if(flag(ix^d,
e_))
then
2176 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2177 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2184 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2187 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2189 w(ixo^s,
e_)=pressure(ixo^s)
2190 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2191 {
do ix^db=ixomin^db,ixomax^db\}
2192 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2193 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2198 if(.not.primitive)
then
2200 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2201 w(ixo^s,
e_)=pressure(ixo^s)
2203 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2207 end subroutine mhd_handle_small_values_semirelati
2209 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2212 logical,
intent(in) :: primitive
2213 integer,
intent(in) :: ixi^
l,ixo^
l
2214 double precision,
intent(inout) :: w(ixi^s,1:nw)
2215 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2216 character(len=*),
intent(in) :: subname
2219 logical :: flag(ixi^s,1:nw)
2221 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2226 {
do ix^db=ixomin^db,ixomax^db\}
2230 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2242 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2244 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2247 {
do ix^db=iximin^db,iximax^db\}
2248 w(ix^d,
e_)=w(ix^d,
e_)&
2249 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2251 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2253 {
do ix^db=iximin^db,iximax^db\}
2254 w(ix^d,
e_)=w(ix^d,
e_)&
2255 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2259 if(.not.primitive)
then
2261 {
do ix^db=ixomin^db,ixomax^db\}
2263 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2264 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2267 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2271 end subroutine mhd_handle_small_values_origin
2273 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2276 logical,
intent(in) :: primitive
2277 integer,
intent(in) :: ixi^
l,ixo^
l
2278 double precision,
intent(inout) :: w(ixi^s,1:nw)
2279 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2280 character(len=*),
intent(in) :: subname
2282 double precision :: rho
2284 logical :: flag(ixi^s,1:nw)
2286 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2291 {
do ix^db=ixomin^db,ixomax^db\}
2296 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2303 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2309 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2311 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2314 {
do ix^db=iximin^db,iximax^db\}
2316 w(ix^d,
e_)=w(ix^d,
e_)&
2317 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2319 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2321 {
do ix^db=iximin^db,iximax^db\}
2323 w(ix^d,
e_)=w(ix^d,
e_)&
2324 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2328 if(.not.primitive)
then
2330 {
do ix^db=ixomin^db,ixomax^db\}
2332 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2333 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2334 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2337 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2341 end subroutine mhd_handle_small_values_split
2343 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2346 logical,
intent(in) :: primitive
2347 integer,
intent(in) :: ixi^
l,ixo^
l
2348 double precision,
intent(inout) :: w(ixi^s,1:nw)
2349 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2350 character(len=*),
intent(in) :: subname
2353 logical :: flag(ixi^s,1:nw)
2355 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2360 {
do ix^db=ixomin^db,ixomax^db\}
2361 if(flag(ix^
d,
rho_))
then
2363 ^
c&w(ix^
d,
m^
c_)=0.d0\
2368 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2373 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2375 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2377 if(.not.primitive)
then
2379 {
do ix^db=ixomin^db,ixomax^db\}
2381 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2384 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2388 end subroutine mhd_handle_small_values_inte
2390 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2393 logical,
intent(in) :: primitive
2394 integer,
intent(in) :: ixi^
l,ixo^
l
2395 double precision,
intent(inout) :: w(ixi^s,1:nw)
2396 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2397 character(len=*),
intent(in) :: subname
2400 logical :: flag(ixi^s,1:nw)
2402 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2407 {
do ix^db=ixomin^db,ixomax^db\}
2411 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2417 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2419 if(.not.primitive)
then
2421 {
do ix^db=ixomin^db,ixomax^db\}
2425 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2429 end subroutine mhd_handle_small_values_noe
2431 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2434 logical,
intent(in) :: primitive
2435 integer,
intent(in) :: ixi^
l,ixo^
l
2436 double precision,
intent(inout) :: w(ixi^s,1:nw)
2437 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2438 character(len=*),
intent(in) :: subname
2441 logical :: flag(ixi^s,1:nw)
2443 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2448 {
do ix^db=ixomin^db,ixomax^db\}
2449 if(flag(ix^
d,
rho_))
then
2451 ^
c&w(ix^
d,
m^
c_)=0.d0\
2456 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2461 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2463 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2465 if(.not.primitive)
then
2467 {
do ix^db=ixomin^db,ixomax^db\}
2469 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2472 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2476 end subroutine mhd_handle_small_values_hde
2482 integer,
intent(in) :: ixi^
l, ixo^
l
2483 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2484 double precision,
intent(out) :: v(ixi^s,
ndir)
2486 double precision :: rho(ixi^s)
2491 rho(ixo^s)=1.d0/rho(ixo^s)
2494 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2500 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2503 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2504 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2505 double precision,
intent(inout) :: cmax(ixi^s)
2507 double precision :: rho, inv_rho, ploc, cfast2, avmincs2, b2, kmax
2513 {
do ix^db=ixomin^db,ixomax^db \}
2526 cfast2=b2*inv_rho+cmax(ix^
d)
2527 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2528 if(avmincs2<zero) avmincs2=zero
2529 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2533 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2535 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2538 {
do ix^db=ixomin^db,ixomax^db \}
2541 ploc=(w(ix^d,
p_)+block%equi_vars(ix^d,
equi_pe0_,b0i))
2550 b2=(^
c&w(ix^d,
b^
c_)**2+)
2551 cfast2=b2*inv_rho+cmax(ix^d)
2552 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2553 if(avmincs2<zero) avmincs2=zero
2554 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2558 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2560 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2564 end subroutine mhd_get_cmax_origin
2567 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
2576 double precision :: adiabs(ixo^s), gammas(ixo^s)
2591 {
do ix^db=ixomin^db,ixomax^db \}
2595 cmax(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*rho**(gammas(ix^
d)-1.d0)
2597 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2598 cfast2=b2*inv_rho+cmax(ix^
d)
2599 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*w(ix^
d,mag(idim))**2*inv_rho
2600 if(avmincs2<zero) avmincs2=zero
2601 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2605 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2607 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2610 end subroutine mhd_get_cmax_origin_noe
2613 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2616 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2617 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2618 double precision,
intent(inout):: cmax(ixi^s)
2620 double precision :: csound, avmincs2, idim_alfven_speed2
2621 double precision :: inv_rho, alfven_speed2, gamma2
2624 {
do ix^db=ixomin^db,ixomax^db \}
2625 inv_rho=1.d0/w(ix^
d,
rho_)
2626 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2627 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2628 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2631 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2634 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2635 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2636 if(avmincs2<zero) avmincs2=zero
2638 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2639 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2642 end subroutine mhd_get_cmax_semirelati
2645 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2649 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2650 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2651 double precision,
intent(inout):: cmax(ixi^s)
2653 double precision :: adiabs(ixo^s), gammas(ixo^s)
2654 double precision :: csound, avmincs2, idim_alfven_speed2
2655 double precision :: inv_rho, alfven_speed2, gamma2
2669 {
do ix^db=ixomin^db,ixomax^db \}
2670 inv_rho=1.d0/w(ix^
d,
rho_)
2671 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2672 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2673 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2674 csound=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
2675 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2678 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2679 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2680 if(avmincs2<zero) avmincs2=zero
2682 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2683 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2686 end subroutine mhd_get_cmax_semirelati_noe
2688 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2691 integer,
intent(in) :: ixi^
l, ixo^
l
2692 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2693 double precision,
intent(inout) :: a2max(
ndim)
2694 double precision :: a2(ixi^s,
ndim,nw)
2695 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2700 hxo^
l=ixo^
l-
kr(i,^
d);
2701 gxo^
l=hxo^
l-
kr(i,^
d);
2702 jxo^
l=ixo^
l+
kr(i,^
d);
2703 kxo^
l=jxo^
l+
kr(i,^
d);
2704 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2705 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2706 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2708 end subroutine mhd_get_a2max
2711 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2714 integer,
intent(in) :: ixi^
l,ixo^
l
2715 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2717 double precision,
intent(inout) :: w(ixi^s,1:nw)
2718 double precision,
intent(out) :: tco_local,tmax_local
2720 double precision,
parameter :: trac_delta=0.25d0
2721 double precision :: te(ixi^s),lts(ixi^s)
2722 double precision,
dimension(1:ndim) :: bdir, bunitvec
2723 double precision,
dimension(ixI^S,1:ndim) :: gradt
2724 double precision :: ltrc,ltrp,altr
2725 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2726 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2729 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2732 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2735 tmax_local=maxval(te(ixo^s))
2743 do ix1=ixomin1,ixomax1
2744 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2745 if(lts(ix1)>trac_delta)
then
2746 tco_local=max(tco_local,te(ix1))
2758 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2759 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2760 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2761 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2763 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2774 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2781 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2786 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2787 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2791 if(bdir(1)/=0.d0)
then
2792 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2794 block%special_values(3)=0.d0
2796 if(bdir(2)/=0.d0)
then
2797 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2799 block%special_values(4)=0.d0
2803 if(bdir(1)/=0.d0)
then
2804 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2805 (bdir(3)/bdir(1))**2)
2807 block%special_values(3)=0.d0
2809 if(bdir(2)/=0.d0)
then
2810 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2811 (bdir(3)/bdir(2))**2)
2813 block%special_values(4)=0.d0
2815 if(bdir(3)/=0.d0)
then
2816 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2817 (bdir(2)/bdir(3))**2)
2819 block%special_values(5)=0.d0
2824 block%special_values(1)=zero
2825 {
do ix^db=ixomin^db,ixomax^db\}
2827 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2829 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2832 if(bdir(1)/=0.d0)
then
2833 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2837 if(bdir(2)/=0.d0)
then
2838 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2843 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2844 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2847 if(bdir(1)/=0.d0)
then
2848 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2852 if(bdir(2)/=0.d0)
then
2853 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2857 if(bdir(3)/=0.d0)
then
2858 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2863 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2864 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2866 if(lts(ix^d)>trac_delta)
then
2867 block%special_values(1)=max(block%special_values(1),te(ix^d))
2870 block%special_values(2)=tmax_local
2889 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2890 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2891 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2895 {
do ix^db=ixpmin^db,ixpmax^db\}
2896 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2898 if(bdir(1)/=0.d0)
then
2899 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2903 if(bdir(2)/=0.d0)
then
2904 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2910 if(bdir(1)/=0.d0)
then
2911 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2915 if(bdir(2)/=0.d0)
then
2916 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2920 if(bdir(3)/=0.d0)
then
2921 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2927 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2929 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2930 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2933 {
do ix^db=ixpmin^db,ixpmax^db\}
2935 if(w(ix^d,iw_mag(1))/=0.d0)
then
2936 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)
2940 if(w(ix^d,iw_mag(2))/=0.d0)
then
2941 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)
2947 if(w(ix^d,iw_mag(1))/=0.d0)
then
2948 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+&
2949 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
2953 if(w(ix^d,iw_mag(2))/=0.d0)
then
2954 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+&
2955 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
2959 if(w(ix^d,iw_mag(3))/=0.d0)
then
2960 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+&
2961 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
2967 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2969 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2970 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2976 {
do ix^db=ixpmin^db,ixpmax^db\}
2978 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
2979 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
2980 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2983 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
2984 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
2985 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
2986 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2992 call mpistop(
"unknown mhd_trac_type")
2995 end subroutine mhd_get_tcutoff
2998 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
3001 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3002 double precision,
intent(in) :: wprim(ixi^s, nw)
3003 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3004 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
3006 double precision :: csound(ixi^s,
ndim)
3007 double precision,
allocatable :: tmp(:^
d&)
3008 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
3012 allocate(tmp(ixa^s))
3014 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
3015 csound(ixa^s,id)=tmp(ixa^s)
3018 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
3019 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
3020 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
3021 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))
3025 ixamax^
d=ixcmax^
d+
kr(id,^
d);
3026 ixamin^
d=ixcmin^
d+
kr(id,^
d);
3027 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)))
3028 ixamax^
d=ixcmax^
d-
kr(id,^
d);
3029 ixamin^
d=ixcmin^
d-
kr(id,^
d);
3030 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)))
3035 ixamax^
d=jxcmax^
d+
kr(id,^
d);
3036 ixamin^
d=jxcmin^
d+
kr(id,^
d);
3037 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)))
3038 ixamax^
d=jxcmax^
d-
kr(id,^
d);
3039 ixamin^
d=jxcmin^
d-
kr(id,^
d);
3040 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)))
3044 end subroutine mhd_get_h_speed
3047 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3050 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3051 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3052 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3053 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3054 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3055 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3056 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3058 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3059 double precision :: umean, dmean, tmp1, tmp2, tmp3
3066 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3067 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3068 if(
present(cmin))
then
3069 {
do ix^db=ixomin^db,ixomax^db\}
3070 tmp1=sqrt(wlp(ix^
d,
rho_))
3071 tmp2=sqrt(wrp(ix^
d,
rho_))
3072 tmp3=1.d0/(tmp1+tmp2)
3073 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3074 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3075 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3076 cmin(ix^
d,1)=umean-dmean
3077 cmax(ix^
d,1)=umean+dmean
3079 if(h_correction)
then
3080 {
do ix^db=ixomin^db,ixomax^db\}
3081 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3082 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3086 {
do ix^db=ixomin^db,ixomax^db\}
3087 tmp1=sqrt(wlp(ix^d,
rho_))
3088 tmp2=sqrt(wrp(ix^d,
rho_))
3089 tmp3=1.d0/(tmp1+tmp2)
3090 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3091 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3092 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3093 cmax(ix^d,1)=abs(umean)+dmean
3097 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3098 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3099 if(
present(cmin))
then
3100 {
do ix^db=ixomin^db,ixomax^db\}
3101 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3102 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3104 if(h_correction)
then
3105 {
do ix^db=ixomin^db,ixomax^db\}
3106 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3107 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3111 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3115 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3116 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3117 if(
present(cmin))
then
3118 {
do ix^db=ixomin^db,ixomax^db\}
3119 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3120 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3121 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3123 if(h_correction)
then
3124 {
do ix^db=ixomin^db,ixomax^db\}
3125 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3126 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3130 {
do ix^db=ixomin^db,ixomax^db\}
3131 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3132 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3137 end subroutine mhd_get_cbounds
3140 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3143 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3144 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3145 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3146 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3147 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3148 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3149 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3151 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3156 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3157 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3159 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3160 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3162 if(
present(cmin))
then
3163 {
do ix^db=ixomin^db,ixomax^db\}
3164 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3165 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3166 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3169 {
do ix^db=ixomin^db,ixomax^db\}
3170 csoundl(ix^d)=max(csoundl(ix^d),csoundr(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)
3175 end subroutine mhd_get_cbounds_semirelati
3178 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3181 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3182 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3183 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3184 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3185 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3186 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3187 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3189 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3190 double precision :: umean, dmean, tmp1, tmp2, tmp3
3197 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3198 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3199 if(
present(cmin))
then
3200 {
do ix^db=ixomin^db,ixomax^db\}
3203 tmp3=1.d0/(tmp1+tmp2)
3204 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3205 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3206 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3207 cmin(ix^
d,1)=umean-dmean
3208 cmax(ix^
d,1)=umean+dmean
3210 if(h_correction)
then
3211 {
do ix^db=ixomin^db,ixomax^db\}
3212 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3213 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3217 {
do ix^db=ixomin^db,ixomax^db\}
3220 tmp3=1.d0/(tmp1+tmp2)
3221 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3222 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3223 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3224 cmax(ix^d,1)=abs(umean)+dmean
3228 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3229 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3230 if(
present(cmin))
then
3231 {
do ix^db=ixomin^db,ixomax^db\}
3232 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3233 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3235 if(h_correction)
then
3236 {
do ix^db=ixomin^db,ixomax^db\}
3237 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3238 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3242 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3246 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3247 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3248 if(
present(cmin))
then
3249 {
do ix^db=ixomin^db,ixomax^db\}
3250 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3251 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3252 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3254 if(h_correction)
then
3255 {
do ix^db=ixomin^db,ixomax^db\}
3256 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3257 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3261 {
do ix^db=ixomin^db,ixomax^db\}
3262 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3263 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3268 end subroutine mhd_get_cbounds_split_rho
3271 subroutine mhd_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3274 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3275 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3276 double precision,
intent(in) :: cmax(ixi^s)
3277 double precision,
intent(in),
optional :: cmin(ixi^s)
3278 type(ct_velocity),
intent(inout):: vcts
3280 end subroutine mhd_get_ct_velocity_average
3282 subroutine mhd_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3285 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3286 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3287 double precision,
intent(in) :: cmax(ixi^s)
3288 double precision,
intent(in),
optional :: cmin(ixi^s)
3289 type(ct_velocity),
intent(inout):: vcts
3291 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3293 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3295 end subroutine mhd_get_ct_velocity_contact
3297 subroutine mhd_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3300 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3301 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3302 double precision,
intent(in) :: cmax(ixi^s)
3303 double precision,
intent(in),
optional :: cmin(ixi^s)
3304 type(ct_velocity),
intent(inout):: vcts
3306 integer :: idime,idimn
3308 if(.not.
allocated(vcts%vbarC))
then
3309 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3310 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3313 if(
present(cmin))
then
3314 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3315 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3317 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3318 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3321 idimn=mod(idim,
ndir)+1
3322 idime=mod(idim+1,
ndir)+1
3324 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3325 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3326 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3327 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3328 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3330 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3331 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3332 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3333 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3334 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3336 end subroutine mhd_get_ct_velocity_hll
3339 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3343 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3344 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3345 double precision,
intent(out):: csound(ixo^s)
3347 double precision :: adiabs(ixo^s), gammas(ixo^s)
3348 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3368 {
do ix^db=ixomin^db,ixomax^db \}
3369 inv_rho=1.d0/w(ix^
d,
rho_)
3373 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3376 cfast2=b2*inv_rho+csound(ix^
d)
3377 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3379 if(avmincs2<zero) avmincs2=zero
3380 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3382 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3386 {
do ix^db=ixomin^db,ixomax^db \}
3387 inv_rho=1.d0/w(ix^d,
rho_)
3391 csound(ix^d)=gammas(ix^d)*adiabs(ix^d)*w(ix^d,
rho_)**(gammas(ix^d)-1.d0)
3393 b2=(^
c&w(ix^d,
b^
c_)**2+)
3394 cfast2=b2*inv_rho+csound(ix^d)
3395 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3396 if(avmincs2<zero) avmincs2=zero
3397 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3399 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3404 end subroutine mhd_get_csound_prim
3407 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3410 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3411 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3412 double precision,
intent(out):: csound(ixo^s)
3414 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3421 {
do ix^db=ixomin^db,ixomax^db \}
3428 cfast2=b2*inv_rho+csound(ix^
d)
3429 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3431 if(avmincs2<zero) avmincs2=zero
3432 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3434 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3438 {
do ix^db=ixomin^db,ixomax^db \}
3444 b2=(^
c&w(ix^d,
b^
c_)**2+)
3445 cfast2=b2*inv_rho+csound(ix^d)
3446 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3447 if(avmincs2<zero) avmincs2=zero
3448 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3450 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3455 end subroutine mhd_get_csound_prim_split
3458 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3461 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3463 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3464 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3466 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3469 {
do ix^db=ixomin^db,ixomax^db\}
3470 inv_rho = 1.d0/w(ix^
d,
rho_)
3473 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3474 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3475 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3476 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3479 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3480 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3481 if(avmincs2<zero) avmincs2=zero
3483 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3486 end subroutine mhd_get_csound_semirelati
3489 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3493 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3495 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3496 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3498 double precision :: adiabs(ixo^s), gammas(ixo^s)
3499 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3512 {
do ix^db=ixomin^db,ixomax^db\}
3513 inv_rho = 1.d0/w(ix^
d,
rho_)
3515 csound(ix^
d)=gammas(ix^
d)*adiabs(ix^
d)*w(ix^
d,
rho_)**(gammas(ix^
d)-1.d0)
3516 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3517 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3518 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3519 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3522 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3523 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3524 if(avmincs2<zero) avmincs2=zero
3526 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3529 end subroutine mhd_get_csound_semirelati_noe
3532 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3536 integer,
intent(in) :: ixi^
l, ixo^
l
3537 double precision,
intent(in) :: w(ixi^s,nw)
3538 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3539 double precision,
intent(out):: pth(ixi^s)
3541 double precision :: adiabs(ixo^s), gammas(ixo^s)
3554 {
do ix^db=ixomin^db,ixomax^db\}
3555 pth(ix^
d)=adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)
3558 end subroutine mhd_get_pthermal_noe
3561 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3565 integer,
intent(in) :: ixi^
l, ixo^
l
3566 double precision,
intent(in) :: w(ixi^s,nw)
3567 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3568 double precision,
intent(out):: pth(ixi^s)
3572 {
do ix^db= ixomin^db,ixomax^db\}
3576 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3581 if(check_small_values.and..not.fix_small_values)
then
3582 {
do ix^db= ixomin^db,ixomax^db\}
3583 if(pth(ix^d)<small_pressure)
then
3584 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3585 " encountered when call mhd_get_pthermal_inte"
3586 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3587 write(*,*)
"Location: ", x(ix^d,:)
3588 write(*,*)
"Cell number: ", ix^d
3590 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3593 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3594 write(*,*)
"Saving status at the previous time step"
3600 end subroutine mhd_get_pthermal_inte
3603 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3607 integer,
intent(in) :: ixi^
l, ixo^
l
3608 double precision,
intent(in) :: w(ixi^s,nw)
3609 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3610 double precision,
intent(out):: pth(ixi^s)
3614 {
do ix^db=ixomin^db,ixomax^db\}
3619 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3620 +(^
c&w(ix^
d,
b^
c_)**2+)))
3625 if(check_small_values.and..not.fix_small_values)
then
3626 {
do ix^db=ixomin^db,ixomax^db\}
3627 if(pth(ix^d)<small_pressure)
then
3628 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3629 " encountered when call mhd_get_pthermal"
3630 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3631 write(*,*)
"Location: ", x(ix^d,:)
3632 write(*,*)
"Cell number: ", ix^d
3634 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3637 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3638 write(*,*)
"Saving status at the previous time step"
3644 end subroutine mhd_get_pthermal_origin
3647 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3651 integer,
intent(in) :: ixi^
l, ixo^
l
3652 double precision,
intent(in) :: w(ixi^s,nw)
3653 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3654 double precision,
intent(out):: pth(ixi^s)
3656 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3659 {
do ix^db=ixomin^db,ixomax^db\}
3660 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3661 if(b2>smalldouble)
then
3669 inv_rho=1.d0/w(ix^
d,
rho_)
3671 b2=b2*inv_rho*inv_squared_c
3673 gamma2=1.d0/(1.d0+b2)
3675 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3679 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3680 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3681 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3685 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3691 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3692 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3693 +(^
c&w(ix^
d,
b^
c_)**2+)&
3694 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3698 if(check_small_values.and..not.fix_small_values)
then
3699 {
do ix^db=ixomin^db,ixomax^db\}
3700 if(pth(ix^d)<small_pressure)
then
3701 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3702 " encountered when call mhd_get_pthermal_semirelati"
3703 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3704 write(*,*)
"Location: ", x(ix^d,:)
3705 write(*,*)
"Cell number: ", ix^d
3707 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3710 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3711 write(*,*)
"Saving status at the previous time step"
3717 end subroutine mhd_get_pthermal_semirelati
3720 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3724 integer,
intent(in) :: ixi^
l, ixo^
l
3725 double precision,
intent(in) :: w(ixi^s,nw)
3726 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3727 double precision,
intent(out):: pth(ixi^s)
3731 {
do ix^db= ixomin^db,ixomax^db\}
3732 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3735 if(check_small_values.and..not.fix_small_values)
then
3736 {
do ix^db= ixomin^db,ixomax^db\}
3737 if(pth(ix^d)<small_pressure)
then
3738 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3739 " encountered when call mhd_get_pthermal_hde"
3740 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3741 write(*,*)
"Location: ", x(ix^d,:)
3742 write(*,*)
"Cell number: ", ix^d
3744 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3747 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3748 write(*,*)
"Saving status at the previous time step"
3754 end subroutine mhd_get_pthermal_hde
3757 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3759 integer,
intent(in) :: ixi^
l, ixo^
l
3760 double precision,
intent(in) :: w(ixi^s, 1:nw)
3761 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3762 double precision,
intent(out):: res(ixi^s)
3763 res(ixo^s) = w(ixo^s,
te_)
3764 end subroutine mhd_get_temperature_from_te
3767 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3769 integer,
intent(in) :: ixi^
l, ixo^
l
3770 double precision,
intent(in) :: w(ixi^s, 1:nw)
3771 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3772 double precision,
intent(out):: res(ixi^s)
3774 double precision :: r(ixi^s)
3777 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3778 end subroutine mhd_get_temperature_from_eint
3781 subroutine mhd_get_temperature_from_etot(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),rho(ixi^s)
3793 res(ixo^s)=res(ixo^s)/(r(ixo^s)*rho(ixo^s))
3795 end subroutine mhd_get_temperature_from_etot
3797 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3799 integer,
intent(in) :: ixi^
l, ixo^
l
3800 double precision,
intent(in) :: w(ixi^s, 1:nw)
3801 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3802 double precision,
intent(out):: res(ixi^s)
3804 double precision :: r(ixi^s)
3810 end subroutine mhd_get_temperature_from_eint_with_equi
3812 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3814 integer,
intent(in) :: ixi^
l, ixo^
l
3815 double precision,
intent(in) :: w(ixi^s, 1:nw)
3816 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3817 double precision,
intent(out):: res(ixi^s)
3819 double precision :: r(ixi^s)
3825 end subroutine mhd_get_temperature_equi
3827 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3829 integer,
intent(in) :: ixi^
l, ixo^
l
3830 double precision,
intent(in) :: w(ixi^s, 1:nw)
3831 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3832 double precision,
intent(out):: res(ixi^s)
3834 end subroutine mhd_get_rho_equi
3836 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3838 integer,
intent(in) :: ixi^
l, ixo^
l
3839 double precision,
intent(in) :: w(ixi^s, 1:nw)
3840 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3841 double precision,
intent(out):: res(ixi^s)
3843 end subroutine mhd_get_pe_equi
3846 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3850 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3852 double precision,
intent(in) :: wc(ixi^s,nw)
3854 double precision,
intent(in) :: w(ixi^s,nw)
3855 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3856 double precision,
intent(out) :: f(ixi^s,nwflux)
3858 double precision :: vhall(ixi^s,1:
ndir)
3859 double precision :: ptotal
3863 {
do ix^db=ixomin^db,ixomax^db\}
3876 {
do ix^db=ixomin^db,ixomax^db\}
3880 ^
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_)\
3881 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3883 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3886 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3887 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3889 ^
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_)\
3893 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3894 {
do ix^db=ixomin^db,ixomax^db\}
3895 if(total_energy)
then
3897 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3898 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3901 ^
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))\
3905 {
do ix^db=ixomin^db,ixomax^db\}
3906 f(ix^d,mag(idim))=w(ix^d,
psi_)
3908 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3913 {
do ix^db=ixomin^db,ixomax^db\}
3919 {
do ix^db=ixomin^db,ixomax^db\}
3920 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)
3925 end subroutine mhd_get_flux
3929 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3934 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3936 double precision,
intent(in) :: wc(ixi^s,nw)
3938 double precision,
intent(in) :: w(ixi^s,nw)
3939 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3940 double precision,
intent(out) :: f(ixi^s,nwflux)
3942 double precision :: vhall(ixi^s,1:
ndir)
3943 double precision :: adiabs(ixo^s), gammas(ixo^s)
3956 {
do ix^db=ixomin^db,ixomax^db\}
3962 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3967 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3968 {
do ix^db=ixomin^db,ixomax^db\}
3970 ^
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))\
3974 {
do ix^db=ixomin^db,ixomax^db\}
3975 f(ix^d,mag(idim))=w(ix^d,
psi_)
3977 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3982 {
do ix^db=ixomin^db,ixomax^db\}
3987 end subroutine mhd_get_flux_noe
3990 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3994 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3996 double precision,
intent(in) :: wc(ixi^s,nw)
3998 double precision,
intent(in) :: w(ixi^s,nw)
3999 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4000 double precision,
intent(out) :: f(ixi^s,nwflux)
4002 double precision :: vhall(ixi^s,1:
ndir)
4005 {
do ix^db=ixomin^db,ixomax^db\}
4018 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4019 {
do ix^db=ixomin^db,ixomax^db\}
4021 ^
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))\
4025 {
do ix^db=ixomin^db,ixomax^db\}
4026 f(ix^d,mag(idim))=w(ix^d,
psi_)
4028 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4033 {
do ix^db=ixomin^db,ixomax^db\}
4039 {
do ix^db=ixomin^db,ixomax^db\}
4040 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)
4045 end subroutine mhd_get_flux_hde
4048 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
4052 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4054 double precision,
intent(in) :: wc(ixi^s,nw)
4056 double precision,
intent(in) :: w(ixi^s,nw)
4057 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4058 double precision,
intent(out) :: f(ixi^s,nwflux)
4060 double precision :: vhall(ixi^s,1:
ndir)
4061 double precision :: ptotal, btotal(ixo^s,1:
ndir)
4064 {
do ix^db=ixomin^db,ixomax^db\}
4072 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
4076 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
4080 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
4081 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4083 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
4087 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4090 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
4097 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
4098 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
4103 {
do ix^db=ixomin^db,ixomax^db\}
4104 f(ix^d,mag(idim))=w(ix^d,
psi_)
4106 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4111 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4112 {
do ix^db=ixomin^db,ixomax^db\}
4114 ^
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))\
4115 if(total_energy)
then
4117 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
4118 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
4124 {
do ix^db=ixomin^db,ixomax^db\}
4129 {
do ix^db=ixomin^db,ixomax^db\}
4130 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
4135 end subroutine mhd_get_flux_split
4138 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
4142 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4144 double precision,
intent(in) :: wc(ixi^s,nw)
4146 double precision,
intent(in) :: w(ixi^s,nw)
4147 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4148 double precision,
intent(out) :: f(ixi^s,nwflux)
4150 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
4153 {
do ix^db=ixomin^db,ixomax^db\}
4158 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4159 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4160 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4165 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4170 e2=(^
c&e(ix^
d,^
c)**2+)
4177 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4178 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4179 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4182 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4183 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4196 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4198 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)
4205 {
do ix^db=ixomin^db,ixomax^db\}
4206 f(ix^d,mag(idim))=w(ix^d,
psi_)
4208 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4213 {
do ix^db=ixomin^db,ixomax^db\}
4218 {
do ix^db=ixomin^db,ixomax^db\}
4219 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)
4224 end subroutine mhd_get_flux_semirelati
4226 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4231 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4233 double precision,
intent(in) :: wc(ixi^s,nw)
4235 double precision,
intent(in) :: w(ixi^s,nw)
4236 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4237 double precision,
intent(out) :: f(ixi^s,nwflux)
4239 double precision :: adiabs(ixo^s), gammas(ixo^s)
4240 double precision :: e(ixo^s,1:
ndir),e2
4253 {
do ix^db=ixomin^db,ixomax^db\}
4258 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4259 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4260 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4261 e2=(^
c&e(ix^
d,^
c)**2+)
4266 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4275 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4277 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+adiabs(ix^
d)*w(ix^
d,
rho_)**gammas(ix^
d)+half*((^
c&w(ix^
d,
b^
c_)**2+)+e2*inv_squared_c)
4284 {
do ix^db=ixomin^db,ixomax^db\}
4285 f(ix^d,mag(idim))=w(ix^d,
psi_)
4287 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4292 {
do ix^db=ixomin^db,ixomax^db\}
4297 end subroutine mhd_get_flux_semirelati_noe
4303 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4305 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4306 double precision,
intent(in) :: qdt
4307 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4308 double precision,
intent(inout) :: w(ixi^s,1:nw)
4309 double precision :: tmp(ixi^s)
4310 double precision :: jxbxb(ixi^s,1:3)
4312 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4315 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4317 end subroutine add_source_ambipolar_internal_energy
4319 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4322 integer,
intent(in) :: ixi^
l, ixo^
l
4323 double precision,
intent(in) :: w(ixi^s,nw)
4324 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4325 double precision,
intent(out) :: res(:^
d&,:)
4327 double precision :: btot(ixi^s,1:3)
4328 double precision :: current(ixi^s,7-2*
ndir:3)
4329 double precision :: tmp(ixi^s),b2(ixi^s)
4330 integer :: idir, idirmin
4339 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4342 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4345 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4346 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4348 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4351 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4353 end subroutine mhd_get_jxbxb
4360 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4364 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4365 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4366 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4367 double precision,
intent(in) :: my_dt
4368 logical,
intent(in) :: fix_conserve_at_step
4370 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4371 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4372 double precision :: fe(ixi^s,
sdim:3)
4373 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4374 integer :: i, ixa^
l, ie_
4380 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4390 btot(ixa^s,1:3)=0.d0
4396 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4399 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4400 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4402 wres(ixo^s,
e_)=-tmp2(ixo^s)
4408 ff(ixa^s,1) = tmp(ixa^s,2)
4409 ff(ixa^s,2) = -tmp(ixa^s,1)
4411 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4412 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4413 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4416 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4418 ixamin^
d=ixomin^
d-1;
4419 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4428 ff(ixa^s,2) = tmp(ixa^s,3)
4429 ff(ixa^s,3) = -tmp(ixa^s,2)
4430 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4431 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4433 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4435 ff(ixa^s,1) = -tmp(ixa^s,3)
4437 ff(ixa^s,3) = tmp(ixa^s,1)
4438 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4439 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4440 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4444 ff(ixa^s,1) = tmp(ixa^s,2)
4445 ff(ixa^s,2) = -tmp(ixa^s,1)
4447 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4448 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4449 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4454 if(fix_conserve_at_step)
then
4455 fluxall=my_dt*fluxall
4462 end subroutine sts_set_source_ambipolar
4465 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4468 integer,
intent(in) :: ixi^
l, ixo^
l
4469 double precision,
intent(in) :: w(ixi^s,1:nw)
4470 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4472 double precision,
intent(in) :: ecc(ixi^s,1:3)
4473 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4474 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4476 integer :: hxc^
l,ixc^
l,ixa^
l
4477 integer :: idim1,idim2,idir,ix^
d
4483 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4485 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4486 ixamin^
d=ixcmin^
d+ix^
d;
4487 ixamax^
d=ixcmax^
d+ix^
d;
4488 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4490 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4496 ixcmin^d=ixomin^d-1;
4504 hxc^l=ixc^l-kr(idim2,^d);
4506 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4507 +lvc(idim1,idim2,idir)&
4512 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4515 end subroutine update_faces_ambipolar
4519 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4522 integer,
intent(in) :: ixi^
l, ixo^
l
4523 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4524 double precision,
intent(out) :: src(ixi^s)
4526 double precision :: ffc(ixi^s,1:
ndim)
4527 double precision :: dxinv(
ndim)
4528 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4534 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4536 ixbmin^
d=ixcmin^
d+ix^
d;
4537 ixbmax^
d=ixcmax^
d+ix^
d;
4540 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4542 ff(ixi^s,1:ndim)=0.d0
4544 ixb^l=ixo^l-kr(idims,^d);
4545 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4547 if({ ix^d==0 .and. ^d==idims | .or.})
then
4548 ixbmin^d=ixcmin^d-ix^d;
4549 ixbmax^d=ixcmax^d-ix^d;
4550 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4553 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4556 if(slab_uniform)
then
4558 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4559 ixb^l=ixo^l-kr(idims,^d);
4560 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4564 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4565 ixb^l=ixo^l-kr(idims,^d);
4566 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4568 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4570 end subroutine get_flux_on_cell_face
4574 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4577 integer,
intent(in) :: ixi^
l, ixo^
l
4578 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4579 double precision,
intent(in) :: w(ixi^s,1:nw)
4580 double precision :: dtnew
4582 double precision :: coef
4583 double precision :: dxarr(
ndim)
4584 double precision :: tmp(ixi^s)
4589 coef = maxval(abs(tmp(ixo^s)))
4596 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4598 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4601 end function get_ambipolar_dt
4609 integer,
intent(in) :: ixi^
l, ixo^
l
4610 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4611 double precision,
intent(inout) :: res(ixi^s)
4612 double precision :: tmp(ixi^s)
4613 double precision :: rho(ixi^s)
4622 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4626 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4633 integer,
intent(in) :: ixi^
l, ixo^
l
4634 double precision,
intent(in) :: qdt,dtfactor
4635 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4636 double precision,
intent(inout) :: w(ixi^s,1:nw)
4637 logical,
intent(in) :: qsourcesplit
4638 logical,
intent(inout) :: active
4645 if (.not. qsourcesplit)
then
4649 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4653 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4659 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4665 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4669 if (abs(
mhd_eta)>smalldouble)
then
4671 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4676 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4682 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4686 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4693 select case (type_divb)
4698 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4701 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4704 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4705 case (divb_janhunen)
4707 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4708 case (divb_lindejanhunen)
4710 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4711 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4712 case (divb_lindepowel)
4714 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4715 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4716 case (divb_lindeglm)
4718 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4719 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4720 case (divb_multigrid)
4725 call mpistop(
'Unknown divB fix')
4732 w,x,qsourcesplit,active,
rc_fl)
4742 w,x,gravity_energy,qsourcesplit,active)
4751 if(.not.qsourcesplit)
then
4753 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4757 end subroutine mhd_add_source
4759 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4763 integer,
intent(in) :: ixi^
l, ixo^
l
4764 double precision,
intent(in) :: qdt,dtfactor
4765 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4766 double precision,
intent(inout) :: w(ixi^s,1:nw)
4767 double precision :: divv(ixi^s)
4783 end subroutine add_pe0_divv
4785 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4787 integer,
intent(in) :: ixi^
l,ixo^
l
4788 double precision,
intent(in) :: qdt
4789 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4790 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4791 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4793 double precision :: r(ixi^s),te(ixi^s),rho_loc(ixi^s),pth_loc(ixi^s)
4794 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,bdir(
ndim),bunitvec(
ndim)
4798 {
do ix^db=iximin^db,iximax^db\}
4803 rho_loc(ix^
d)=wctprim(ix^
d,
rho_)
4804 pth_loc(ix^
d)=wctprim(ix^
d,
p_)
4806 te(ix^
d)=pth_loc(ix^
d)/(r(ix^
d)*rho_loc(ix^
d))
4812 do ix1=ixomin1,ixomax1
4814 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4816 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4819 sigma_t7=sigma_t5*te(ix^d)
4823 sigma_t7=sigma_t5*te(ix^d)
4825 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)
4827 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4828 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4829 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4831 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4832 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4837 do ix2=ixomin2,ixomax2
4838 do ix1=ixomin1,ixomax1
4840 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4842 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4845 sigma_t7=sigma_t5*te(ix^d)
4849 sigma_t7=sigma_t5*te(ix^d)
4852 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4854 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4856 if(bdir(1)/=0.d0)
then
4857 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
4861 if(bdir(2)/=0.d0)
then
4862 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
4866 sigmat5_bgradt=sigma_t5*(&
4867 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)&
4868 +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))
4870 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4871 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4872 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4874 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4875 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4881 do ix3=ixomin3,ixomax3
4882 do ix2=ixomin2,ixomax2
4883 do ix1=ixomin1,ixomax1
4885 if(te(ix^d)<block%wextra(ix^d,
tcoff_))
then
4887 sigma_t7=sigma_t5*block%wextra(ix^d,
tcoff_)
4890 sigma_t7=sigma_t5*te(ix^d)
4894 sigma_t7=sigma_t5*te(ix^d)
4897 ^d&bdir(^d)=wct({ix^d},mag(^d))+block%B0({ix^d},^d,0)\
4899 ^d&bdir(^d)=wct({ix^d},mag(^d))\
4901 if(bdir(1)/=0.d0)
then
4902 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
4906 if(bdir(2)/=0.d0)
then
4907 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
4911 if(bdir(3)/=0.d0)
then
4912 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
4916 sigmat5_bgradt=sigma_t5*(&
4917 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)&
4918 +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)&
4919 +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))
4921 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^d)*(
mhd_gamma*te(ix^d))**1.5d0)
4922 tau=max(4.d0*dt, f_sat*sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4923 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^d,
q_))/tau
4925 w(ix^d,
q_)=w(ix^d,
q_)-qdt*(sigmat5_bgradt+wct(ix^d,
q_))/&
4926 max(4.d0*dt, sigma_t7*courantpar**2/(pth_loc(ix^d)*inv_gamma_1*cmax_global**2))
4932 end subroutine add_hypertc_source
4935 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4937 integer,
intent(in) :: ixi^
l, ixo^
l
4938 double precision,
intent(in) :: w(ixi^s,1:nw)
4939 double precision,
intent(inout) :: jxb(ixi^s,3)
4940 double precision :: a(ixi^s,3),
b(ixi^s,3)
4942 double precision :: current(ixi^s,7-2*
ndir:3)
4943 integer :: idir, idirmin
4948 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4952 b(ixo^s, idir) = w(ixo^s,mag(idir))
4961 a(ixo^s,idir)=current(ixo^s,idir)
4965 end subroutine get_lorentz_force
4969 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4971 integer,
intent(in) :: ixi^
l, ixo^
l
4972 double precision,
intent(in) :: w(ixi^s, nw)
4973 double precision,
intent(out) :: gamma_a2(ixo^s)
4974 double precision :: rho(ixi^s)
4980 rho(ixo^s) = w(ixo^s,
rho_)
4983 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4984 end subroutine mhd_gamma2_alfven
4988 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4990 integer,
intent(in) :: ixi^
l, ixo^
l
4991 double precision,
intent(in) :: w(ixi^s, nw)
4992 double precision :: gamma_a(ixo^s)
4994 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4995 gamma_a = sqrt(gamma_a)
4996 end function mhd_gamma_alfven
5000 integer,
intent(in) :: ixi^
l, ixo^
l
5001 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
5002 double precision,
intent(out) :: rho(ixi^s)
5007 rho(ixo^s) = w(ixo^s,
rho_)
5013 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
5016 integer,
intent(in) :: ixi^
l,ixo^
l, ie
5017 double precision,
intent(inout) :: w(ixi^s,1:nw)
5018 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5019 character(len=*),
intent(in) :: subname
5021 double precision :: rho(ixi^s)
5023 logical :: flag(ixi^s,1:nw)
5027 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
5028 flag(ixo^s,ie)=.true.
5030 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
5032 if(any(flag(ixo^s,ie)))
then
5036 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
5039 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
5045 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
5048 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
5054 end subroutine mhd_handle_small_ei
5056 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
5060 integer,
intent(in) :: ixi^
l, ixo^
l
5061 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5062 double precision,
intent(inout) :: w(ixi^s,1:nw)
5064 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
5073 end subroutine mhd_update_temperature
5076 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
5079 integer,
intent(in) :: ixi^
l, ixo^
l
5080 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5081 double precision,
intent(inout) :: w(ixi^s,1:nw)
5083 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
5095 a(ixo^s,idir)=
block%J0(ixo^s,idir)
5100 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5103 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5109 if(total_energy)
then
5112 b(ixo^s,:)=wct(ixo^s,mag(:))
5121 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5124 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5128 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5132 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5137 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5142 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5144 end subroutine add_source_b0split
5147 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5151 integer,
intent(in) :: ixi^
l, ixo^
l
5152 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5153 double precision,
intent(inout) :: w(ixi^s,1:nw)
5154 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5156 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5157 integer :: idir, idirmin, ix^
d
5161 {
do ix^db=iximin^db,iximax^db\}
5163 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5164 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5165 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5167 call divvector(e,ixi^l,ixo^l,dive)
5169 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5172 {
do ix^db=ixomin^db,ixomax^db\}
5173 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5174 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5175 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5176 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5177 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5178 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5182 end subroutine add_source_semirelativistic
5185 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5189 integer,
intent(in) :: ixi^
l, ixo^
l
5190 double precision,
intent(in) :: qdt
5191 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5192 double precision,
intent(inout) :: w(ixi^s,1:nw)
5193 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5195 double precision :: divv(ixi^s), tmp
5207 {
do ix^db=ixomin^db,ixomax^db\}
5209 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5210 if(w(ix^
d,
e_)<small_e)
then
5215 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
5218 if(fix_small_values)
then
5219 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5221 end subroutine add_source_internal_e
5224 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5229 integer,
intent(in) :: ixi^
l, ixo^
l
5230 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5231 double precision,
intent(inout) :: w(ixi^s,1:nw)
5232 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5234 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5235 double precision :: current(ixi^s,7-2*
ndir:3)
5236 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5237 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5238 integer :: idir, idirmin, idims, ix^
d
5243 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5251 j(ixo^s,idir)=current(ixo^s,idir)
5269 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
5275 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5282 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)
5286 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5290 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5291 tmp(ixo^s)=sqrt(b2(ixo^s))
5292 where(tmp(ixo^s)>smalldouble)
5293 tmp(ixo^s)=1.d0/tmp(ixo^s)
5299 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5304 {
do ix^db=ixomin^db,ixomax^db\}
5306 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5308 b2(ix^
d)=vaoc/(1.d0+vaoc)
5311 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5314 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5318 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5321 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5322 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5325 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5328 end subroutine add_source_hydrodynamic_e
5334 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5339 integer,
intent(in) :: ixi^
l, ixo^
l
5340 double precision,
intent(in) :: qdt
5341 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5342 double precision,
intent(inout) :: w(ixi^s,1:nw)
5343 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5344 integer :: lxo^
l, kxo^
l
5346 double precision :: tmp(ixi^s),tmp2(ixi^s)
5349 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5350 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5359 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5360 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5367 gradeta(ixo^s,1:
ndim)=zero
5373 gradeta(ixo^s,idim)=tmp(ixo^s)
5380 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5387 tmp2(ixi^s)=bf(ixi^s,idir)
5389 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5390 jxo^
l=ixo^
l+
kr(idim,^
d);
5391 hxo^
l=ixo^
l-
kr(idim,^
d);
5392 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5393 tmp(ixo^s)=tmp(ixo^s)+&
5394 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5399 tmp2(ixi^s)=bf(ixi^s,idir)
5401 jxo^
l=ixo^
l+
kr(idim,^
d);
5402 hxo^
l=ixo^
l-
kr(idim,^
d);
5403 tmp(ixo^s)=tmp(ixo^s)+&
5404 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5409 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5413 do jdir=1,
ndim;
do kdir=idirmin,3
5414 if (
lvc(idir,jdir,kdir)/=0)
then
5415 if (
lvc(idir,jdir,kdir)==1)
then
5416 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5418 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5425 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5426 if(total_energy)
then
5427 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5433 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5436 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5438 end subroutine add_source_res1
5442 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5447 integer,
intent(in) :: ixi^
l, ixo^
l
5448 double precision,
intent(in) :: qdt
5449 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5450 double precision,
intent(inout) :: w(ixi^s,1:nw)
5453 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5454 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5455 integer :: ixa^
l,idir,idirmin,idirmin1
5459 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5460 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5470 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5475 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5484 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5487 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5492 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5494 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5496 if(total_energy)
then
5499 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5500 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5503 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5507 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5508 end subroutine add_source_res2
5512 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5516 integer,
intent(in) :: ixi^
l, ixo^
l
5517 double precision,
intent(in) :: qdt
5518 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5519 double precision,
intent(inout) :: w(ixi^s,1:nw)
5521 double precision :: current(ixi^s,7-2*
ndir:3)
5522 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5523 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5526 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5527 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5530 tmpvec(ixa^s,1:
ndir)=zero
5532 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5536 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5539 tmpvec(ixa^s,1:
ndir)=zero
5540 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5544 tmpvec2(ixa^s,1:
ndir)=zero
5545 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5548 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5551 if(total_energy)
then
5554 tmpvec2(ixa^s,1:
ndir)=zero
5555 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5556 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5557 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5558 end do;
end do;
end do
5560 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5561 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5564 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5566 end subroutine add_source_hyperres
5568 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5575 integer,
intent(in) :: ixi^
l, ixo^
l
5576 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5577 double precision,
intent(inout) :: w(ixi^s,1:nw)
5579 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5600 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5603 if(total_energy)
then
5612 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5621 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5625 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5627 end subroutine add_source_glm
5630 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5633 integer,
intent(in) :: ixi^
l, ixo^
l
5634 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5635 double precision,
intent(inout) :: w(ixi^s,1:nw)
5637 double precision :: divb(ixi^s), ba(1:
ndir)
5638 integer :: idir, ix^
d
5644 {
do ix^db=ixomin^db,ixomax^db\}
5649 if (total_energy)
then
5655 {
do ix^db=ixomin^db,ixomax^db\}
5657 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5659 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5660 if (total_energy)
then
5662 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5667 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5669 end subroutine add_source_powel
5671 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5676 integer,
intent(in) :: ixi^
l, ixo^
l
5677 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5678 double precision,
intent(inout) :: w(ixi^s,1:nw)
5680 double precision :: divb(ixi^s)
5681 integer :: idir, ix^
d
5686 {
do ix^db=ixomin^db,ixomax^db\}
5691 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5693 end subroutine add_source_janhunen
5695 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5700 integer,
intent(in) :: ixi^
l, ixo^
l
5701 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5702 double precision,
intent(inout) :: w(ixi^s,1:nw)
5704 double precision :: divb(ixi^s),graddivb(ixi^s)
5705 integer :: idim, idir, ixp^
l, i^
d, iside
5706 logical,
dimension(-1:1^D&) :: leveljump
5714 if(i^
d==0|.and.) cycle
5715 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5716 leveljump(i^
d)=.true.
5718 leveljump(i^
d)=.false.
5727 i^dd=kr(^dd,^d)*(2*iside-3);
5728 if (leveljump(i^dd))
then
5730 ixpmin^d=ixomin^d-i^d
5732 ixpmax^d=ixomax^d-i^d
5743 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5745 {
do i^db=ixpmin^db,ixpmax^db\}
5747 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5749 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5751 if (typedivbdiff==
'all' .and. total_energy)
then
5753 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5758 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5760 end subroutine add_source_linde
5767 integer,
intent(in) :: ixi^
l, ixo^
l
5768 double precision,
intent(in) :: w(ixi^s,1:nw)
5769 double precision :: divb(ixi^s), dsurface(ixi^s)
5771 double precision :: invb(ixo^s)
5772 integer :: ixa^
l,idims
5774 call get_divb(w,ixi^
l,ixo^
l,divb)
5776 where(invb(ixo^s)/=0.d0)
5777 invb(ixo^s)=1.d0/invb(ixo^s)
5780 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5782 ixamin^
d=ixomin^
d-1;
5783 ixamax^
d=ixomax^
d-1;
5784 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5786 ixa^
l=ixo^
l-
kr(idims,^
d);
5787 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5789 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5790 block%dvolume(ixo^s)/dsurface(ixo^s)
5801 integer,
intent(in) :: ixo^
l, ixi^
l
5802 double precision,
intent(in) :: w(ixi^s,1:nw)
5803 integer,
intent(out) :: idirmin
5806 double precision :: current(ixi^s,7-2*
ndir:3)
5807 integer :: idir, idirmin0
5813 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5814 block%J0(ixo^s,idirmin0:3)
5818 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5826 integer,
intent(in) :: ixi^
l, ixo^
l
5827 double precision,
intent(inout) :: dtnew
5828 double precision,
intent(in) ::
dx^
d
5829 double precision,
intent(in) :: w(ixi^s,1:nw)
5830 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5832 double precision :: dxarr(
ndim)
5833 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5834 integer :: idirmin,idim
5848 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5851 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5877 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5884 end subroutine mhd_get_dt
5887 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5893 integer,
intent(in) :: ixi^
l, ixo^
l
5894 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5895 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5897 double precision :: adiabs(ixo^s), gammas(ixo^s)
5898 double precision :: tmp,tmp1,invr,cot
5900 integer :: mr_,mphi_
5901 integer :: br_,bphi_
5904 br_=mag(1); bphi_=mag(1)-1+
phi_
5921 {
do ix^db=ixomin^db,ixomax^db\}
5924 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5929 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5931 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5934 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5935 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5936 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5937 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5938 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5940 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5941 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5942 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5945 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5950 {
do ix^db=ixomin^db,ixomax^db\}
5952 if(local_timestep)
then
5953 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5958 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5960 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5964 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5967 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5968 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5972 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5978 cot=1.d0/tan(x(ix^d,2))
5982 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5983 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5985 if(.not.stagger_grid)
then
5986 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5988 tmp=tmp+wprim(ix^d,
psi_)*cot
5990 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5995 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5996 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5997 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5999 if(.not.stagger_grid)
then
6000 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6002 tmp=tmp+wprim(ix^d,
psi_)*cot
6004 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6007 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6008 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6009 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6010 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6011 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6013 if(.not.stagger_grid)
then
6014 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6015 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6016 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6017 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6018 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6025 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6028 end subroutine mhd_add_source_geom
6031 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6037 integer,
intent(in) :: ixi^
l, ixo^
l
6038 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6039 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6041 double precision :: adiabs(ixo^s), gammas(ixo^s)
6042 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
6044 integer :: mr_,mphi_
6045 integer :: br_,bphi_
6048 br_=mag(1); bphi_=mag(1)-1+
phi_
6065 {
do ix^db=ixomin^db,ixomax^db\}
6068 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6075 tmp=adiabs(ix^
d)*wprim(ix^
d,
rho_)**gammas(ix^
d)
6079 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
6080 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
6081 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6086 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
6092 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
6093 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
6094 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
6095 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6096 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6097 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
6099 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6100 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6101 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6104 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
6105 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
6110 {
do ix^db=ixomin^db,ixomax^db\}
6112 if(local_timestep)
then
6113 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
6119 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
6120 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
6121 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6125 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6132 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6134 tmp1=adiabs(ix^d)*wprim(ix^d,
rho_)**gammas(ix^d)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6138 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
6141 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
6142 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
6143 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
6147 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
6153 cot=1.d0/tan(x(ix^d,2))
6157 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
6158 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
6160 if(.not.stagger_grid)
then
6161 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6163 tmp=tmp+wprim(ix^d,
psi_)*cot
6165 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6171 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6172 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
6173 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6174 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
6176 if(.not.stagger_grid)
then
6177 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6179 tmp=tmp+wprim(ix^d,
psi_)*cot
6181 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6184 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6185 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6186 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6187 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
6188 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6189 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6190 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
6192 if(.not.stagger_grid)
then
6193 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6194 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6195 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6196 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6197 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6204 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6207 end subroutine mhd_add_source_geom_semirelati
6210 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6215 integer,
intent(in) :: ixi^
l, ixo^
l
6216 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6217 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6219 double precision :: tmp,tmp1,tmp2,invr,cot
6221 integer :: mr_,mphi_
6222 integer :: br_,bphi_
6225 br_=mag(1); bphi_=mag(1)-1+
phi_
6230 {
do ix^db=ixomin^db,ixomax^db\}
6233 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6237 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6239 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6240 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6241 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6242 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6243 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6245 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6246 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6247 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6250 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6255 {
do ix^db=ixomin^db,ixomax^db\}
6257 if(local_timestep)
then
6258 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6262 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6263 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6266 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6270 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6271 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6272 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6274 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6275 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6280 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6286 cot=1.d0/tan(x(ix^d,2))
6291 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6292 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6293 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6295 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6296 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6299 if(.not.stagger_grid)
then
6301 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6302 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6304 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6307 tmp=tmp+wprim(ix^d,
psi_)*cot
6309 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6315 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6316 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6317 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6318 +(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)
6320 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6321 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6322 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6325 if(.not.stagger_grid)
then
6327 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6328 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6330 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6333 tmp=tmp+wprim(ix^d,
psi_)*cot
6335 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6339 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6340 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6341 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6342 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6343 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6344 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6345 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6346 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6347 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6349 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6350 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6351 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6352 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6353 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6356 if(.not.stagger_grid)
then
6358 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6359 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6360 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6361 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6362 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6363 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6364 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6365 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6366 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6368 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6369 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6370 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6371 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6372 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6380 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6383 end subroutine mhd_add_source_geom_split
6388 integer,
intent(in) :: ixi^
l, ixo^
l
6389 double precision,
intent(in) :: w(ixi^s, nw)
6390 double precision :: mge(ixo^s)
6393 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6395 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6399 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6402 integer,
intent(in) :: ixi^
l, ixo^
l
6403 double precision,
intent(in) :: w(ixi^s,nw)
6404 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6405 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6407 double precision :: current(ixi^s,7-2*
ndir:3)
6408 double precision :: rho(ixi^s)
6409 integer :: idir, idirmin, ix^
d
6414 do idir = idirmin,
ndir
6415 {
do ix^db=ixomin^db,ixomax^db\}
6416 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6420 end subroutine mhd_getv_hall
6422 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6425 integer,
intent(in) :: ixi^
l, ixo^
l
6426 double precision,
intent(in) :: w(ixi^s,nw)
6427 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6428 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6431 double precision :: current(ixi^s,7-2*
ndir:3)
6432 integer :: idir, idirmin
6439 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6440 do idir = idirmin, 3
6444 end subroutine mhd_get_jambi
6446 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6449 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6450 double precision,
intent(in) :: qt
6451 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6452 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6455 double precision :: db(ixo^s), dpsi(ixo^s)
6459 {
do ix^db=ixomin^db,ixomax^db\}
6460 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6461 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6462 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6463 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6472 {
do ix^db=ixomin^db,ixomax^db\}
6473 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6474 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6475 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6476 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6477 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6479 if(total_energy)
then
6480 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6481 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6483 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6485 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6488 if(total_energy)
then
6489 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6490 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6495 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6497 end subroutine mhd_modify_wlr
6499 subroutine mhd_boundary_adjust(igrid,psb)
6501 integer,
intent(in) :: igrid
6504 integer :: ib, idims, iside, ixo^
l, i^
d
6513 i^
d=
kr(^
d,idims)*(2*iside-3);
6514 if (neighbor_type(i^
d,igrid)/=1) cycle
6515 ib=(idims-1)*2+iside
6533 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6538 end subroutine mhd_boundary_adjust
6540 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6543 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6544 double precision,
intent(inout) :: w(ixg^s,1:nw)
6545 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6547 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6548 integer :: ix^
d,ixf^
l
6561 do ix1=ixfmax1,ixfmin1,-1
6562 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6563 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6564 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6567 do ix1=ixfmax1,ixfmin1,-1
6568 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6569 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6570 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6571 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6572 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6573 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6574 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6588 do ix1=ixfmax1,ixfmin1,-1
6589 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6590 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6591 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6592 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6593 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6594 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6597 do ix1=ixfmax1,ixfmin1,-1
6598 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6599 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6600 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6601 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6602 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6603 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6604 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6605 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6606 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6607 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6608 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6609 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6610 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6611 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6612 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6613 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6614 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6615 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6629 do ix1=ixfmin1,ixfmax1
6630 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6631 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6632 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6635 do ix1=ixfmin1,ixfmax1
6636 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6637 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6638 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6639 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6640 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6641 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6642 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6656 do ix1=ixfmin1,ixfmax1
6657 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6658 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6659 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6660 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6661 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6662 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6665 do ix1=ixfmin1,ixfmax1
6666 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6667 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6668 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6669 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6670 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6671 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6672 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6673 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6674 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6675 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6676 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6677 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6678 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6679 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6680 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6681 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6682 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6683 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6697 do ix2=ixfmax2,ixfmin2,-1
6698 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6699 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6700 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6703 do ix2=ixfmax2,ixfmin2,-1
6704 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6705 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6706 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6707 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6708 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6709 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6710 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6724 do ix2=ixfmax2,ixfmin2,-1
6725 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6726 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6727 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6728 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6729 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6730 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6733 do ix2=ixfmax2,ixfmin2,-1
6734 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6735 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6736 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6737 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6738 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6739 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6740 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6741 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6742 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6743 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6744 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6745 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6746 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6747 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6748 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6749 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6750 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6751 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6765 do ix2=ixfmin2,ixfmax2
6766 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6767 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6768 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6771 do ix2=ixfmin2,ixfmax2
6772 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6773 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6774 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6775 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6776 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6777 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6778 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6792 do ix2=ixfmin2,ixfmax2
6793 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6794 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6795 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6796 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6797 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6798 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6801 do ix2=ixfmin2,ixfmax2
6802 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6803 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6804 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6805 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6806 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6807 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6808 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6809 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6810 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6811 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6812 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6813 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6814 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6815 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6816 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6817 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6818 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6819 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6836 do ix3=ixfmax3,ixfmin3,-1
6837 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6838 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6839 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6840 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6841 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6842 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6845 do ix3=ixfmax3,ixfmin3,-1
6846 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6847 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6848 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6849 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6850 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6851 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6852 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6853 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6854 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6855 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6856 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6857 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6858 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6859 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6860 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6861 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6862 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6863 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6878 do ix3=ixfmin3,ixfmax3
6879 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6880 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6881 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6882 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6883 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6884 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6887 do ix3=ixfmin3,ixfmax3
6888 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6889 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6890 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6891 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6892 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6893 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6894 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6895 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6896 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6897 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6898 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6899 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6900 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6901 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6902 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6903 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6904 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6905 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6911 call mpistop(
"Special boundary is not defined for this region")
6914 end subroutine fixdivb_boundary
6923 double precision,
intent(in) :: qdt
6924 double precision,
intent(in) :: qt
6925 logical,
intent(inout) :: active
6928 integer,
parameter :: max_its = 50
6929 double precision :: residual_it(max_its), max_divb
6930 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6931 double precision :: res
6932 double precision,
parameter :: max_residual = 1
d-3
6933 double precision,
parameter :: residual_reduction = 1
d-10
6934 integer :: iigrid, igrid
6935 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6938 mg%operator_type = mg_laplacian
6946 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6947 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6950 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6951 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6953 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6954 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6957 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6958 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6962 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6963 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6964 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6972 do iigrid = 1, igridstail
6973 igrid = igrids(iigrid);
6976 lvl =
mg%boxes(id)%lvl
6977 nc =
mg%box_size_lvl(lvl)
6983 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6985 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6986 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6991 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6994 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6995 if (
mype == 0) print *,
"iteration vs residual"
6998 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6999 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
7000 if (residual_it(n) < residual_reduction * max_divb)
exit
7002 if (
mype == 0 .and. n > max_its)
then
7003 print *,
"divb_multigrid warning: not fully converged"
7004 print *,
"current amplitude of divb: ", residual_it(max_its)
7005 print *,
"multigrid smallest grid: ", &
7006 mg%domain_size_lvl(:,
mg%lowest_lvl)
7007 print *,
"note: smallest grid ideally has <= 8 cells"
7008 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
7009 print *,
"note: dx/dy/dz should be similar"
7013 call mg_fas_vcycle(
mg, max_res=res)
7014 if (res < max_residual)
exit
7016 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
7021 do iigrid = 1, igridstail
7022 igrid = igrids(iigrid);
7031 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
7035 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
7037 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
7039 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
7052 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
7053 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
7056 if(total_energy)
then
7058 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
7061 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
7071 subroutine mhd_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7075 integer,
intent(in) :: ixi^
l, ixo^
l
7076 double precision,
intent(in) :: qt,qdt
7078 double precision,
intent(in) :: wp(ixi^s,1:nw)
7079 type(state) :: sct, s
7080 type(ct_velocity) :: vcts
7081 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7082 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7084 double precision :: circ(ixi^s,1:
ndim)
7086 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7087 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
7088 integer :: idim1,idim2,idir,iwdim1,iwdim2
7090 associate(bfaces=>s%ws,x=>s%x)
7097 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7104 i1kr^
d=
kr(idim1,^
d);
7107 i2kr^
d=
kr(idim2,^
d);
7110 if (
lvc(idim1,idim2,idir)==1)
then
7112 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7114 {
do ix^db=ixcmin^db,ixcmax^db\}
7115 fe(ix^
d,idir)=quarter*&
7116 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
7117 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
7119 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
7124 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
7132 if(
associated(usr_set_electric_field)) &
7133 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7135 circ(ixi^s,1:ndim)=zero
7140 ixcmin^d=ixomin^d-kr(idim1,^d);
7142 ixa^l=ixc^l-kr(idim2,^d);
7145 if(lvc(idim1,idim2,idir)==1)
then
7147 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7150 else if(lvc(idim1,idim2,idir)==-1)
then
7152 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7158 {
do ix^db=ixcmin^db,ixcmax^db\}
7160 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7162 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7169 end subroutine mhd_update_faces_average
7172 subroutine mhd_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7177 integer,
intent(in) :: ixi^
l, ixo^
l
7178 double precision,
intent(in) :: qt, qdt
7180 double precision,
intent(in) :: wp(ixi^s,1:nw)
7181 type(state) :: sct, s
7182 type(ct_velocity) :: vcts
7183 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7184 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7186 double precision :: circ(ixi^s,1:
ndim)
7188 double precision :: ecc(ixi^s,
sdim:3)
7189 double precision :: ein(ixi^s,
sdim:3)
7191 double precision :: el(ixi^s),er(ixi^s)
7193 double precision :: elc,erc
7195 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7197 double precision :: jce(ixi^s,
sdim:3)
7199 double precision :: xs(ixgs^t,1:
ndim)
7200 double precision :: gradi(ixgs^t)
7201 integer :: ixc^
l,ixa^
l
7202 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7204 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7207 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7213 {
do ix^db=iximin^db,iximax^db\}
7216 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_)
7217 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_)
7218 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_)
7221 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7228 {
do ix^db=iximin^db,iximax^db\}
7231 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7232 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7233 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7236 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7250 i1kr^d=kr(idim1,^d);
7253 i2kr^d=kr(idim2,^d);
7256 if (lvc(idim1,idim2,idir)==1)
then
7258 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7261 {
do ix^db=ixcmin^db,ixcmax^db\}
7262 fe(ix^d,idir)=quarter*&
7263 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7264 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7269 ixamax^d=ixcmax^d+i1kr^d;
7270 {
do ix^db=ixamin^db,ixamax^db\}
7271 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7272 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7275 do ix^db=ixcmin^db,ixcmax^db\}
7276 if(vnorm(ix^d,idim1)>0.d0)
then
7278 else if(vnorm(ix^d,idim1)<0.d0)
then
7279 elc=el({ix^d+i1kr^d})
7281 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7283 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7285 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7286 erc=er({ix^d+i1kr^d})
7288 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7290 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7295 ixamax^d=ixcmax^d+i2kr^d;
7296 {
do ix^db=ixamin^db,ixamax^db\}
7297 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7298 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7301 do ix^db=ixcmin^db,ixcmax^db\}
7302 if(vnorm(ix^d,idim2)>0.d0)
then
7304 else if(vnorm(ix^d,idim2)<0.d0)
then
7305 elc=el({ix^d+i2kr^d})
7307 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7309 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7311 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7312 erc=er({ix^d+i2kr^d})
7314 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7316 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7320 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7325 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7339 if (lvc(idim1,idim2,idir)==0) cycle
7341 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7342 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7345 xs(ixa^s,:)=x(ixa^s,:)
7346 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7347 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7348 if (lvc(idim1,idim2,idir)==1)
then
7349 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7351 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7358 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7360 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7364 {
do ix^db=ixomin^db,ixomax^db\}
7365 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7366 +ein(ix1,ix2-1,ix3-1,idir))
7367 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7368 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7370 else if(idir==2)
then
7371 {
do ix^db=ixomin^db,ixomax^db\}
7372 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7373 +ein(ix1-1,ix2,ix3-1,idir))
7374 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7375 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7378 {
do ix^db=ixomin^db,ixomax^db\}
7379 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7380 +ein(ix1-1,ix2-1,ix3,idir))
7381 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7382 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7388 {
do ix^db=ixomin^db,ixomax^db\}
7389 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7390 +ein(ix1-1,ix2-1,idir))
7391 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7392 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7397 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7403 if(
associated(usr_set_electric_field)) &
7404 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7406 circ(ixi^s,1:ndim)=zero
7411 ixcmin^d=ixomin^d-kr(idim1,^d);
7413 ixa^l=ixc^l-kr(idim2,^d);
7416 if(lvc(idim1,idim2,idir)==1)
then
7418 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7421 else if(lvc(idim1,idim2,idir)==-1)
then
7423 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7429 {
do ix^db=ixcmin^db,ixcmax^db\}
7431 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7433 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7440 end subroutine mhd_update_faces_contact
7443 subroutine mhd_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7448 integer,
intent(in) :: ixi^
l, ixo^
l
7449 double precision,
intent(in) :: qt, qdt
7451 double precision,
intent(in) :: wp(ixi^s,1:nw)
7452 type(state) :: sct, s
7453 type(ct_velocity) :: vcts
7454 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7455 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7457 double precision :: vtill(ixi^s,2)
7458 double precision :: vtilr(ixi^s,2)
7459 double precision :: bfacetot(ixi^s,
ndim)
7460 double precision :: btill(ixi^s,
ndim)
7461 double precision :: btilr(ixi^s,
ndim)
7462 double precision :: cp(ixi^s,2)
7463 double precision :: cm(ixi^s,2)
7464 double precision :: circ(ixi^s,1:
ndim)
7466 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7467 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7468 integer :: idim1,idim2,idir,ix^
d
7470 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7471 cbarmax=>vcts%cbarmax)
7484 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7500 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7504 idim2=mod(idir+1,3)+1
7506 jxc^
l=ixc^
l+
kr(idim1,^
d);
7507 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7511 vtill(ixi^s,2),vtilr(ixi^s,2))
7514 vtill(ixi^s,1),vtilr(ixi^s,1))
7520 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7521 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7523 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7524 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7527 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7530 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7534 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7535 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7537 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7538 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7542 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7543 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7544 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7545 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7546 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7547 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7548 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7549 /(cp(ixc^s,2)+cm(ixc^s,2))
7552 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7556 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7570 circ(ixi^s,1:
ndim)=zero
7575 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7579 if(
lvc(idim1,idim2,idir)/=0)
then
7580 hxc^
l=ixc^
l-
kr(idim2,^
d);
7582 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7583 +
lvc(idim1,idim2,idir)&
7589 {
do ix^db=ixcmin^db,ixcmax^db\}
7591 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7593 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7599 end subroutine mhd_update_faces_hll
7602 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
7607 integer,
intent(in) :: ixi^
l, ixo^
l
7609 double precision,
intent(in) :: wp(ixi^s,1:nw)
7610 type(state),
intent(in) :: sct, s
7612 double precision :: jce(ixi^s,
sdim:3)
7615 double precision :: jcc(ixi^s,7-2*
ndir:3)
7617 double precision :: xs(ixgs^t,1:
ndim)
7619 double precision :: eta(ixi^s)
7620 double precision :: gradi(ixgs^t)
7621 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7623 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7629 if (
lvc(idim1,idim2,idir)==0) cycle
7631 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7632 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7635 xs(ixb^s,:)=x(ixb^s,:)
7636 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7637 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7638 if (
lvc(idim1,idim2,idir)==1)
then
7639 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7641 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7648 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7656 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7657 jcc(ixc^s,idir)=0.d0
7659 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7660 ixamin^
d=ixcmin^
d+ix^
d;
7661 ixamax^
d=ixcmax^
d+ix^
d;
7662 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7664 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7665 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7670 end subroutine get_resistive_electric_field
7673 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7676 integer,
intent(in) :: ixi^
l, ixo^
l
7677 double precision,
intent(in) :: w(ixi^s,1:nw)
7678 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7679 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7681 double precision :: jxbxb(ixi^s,1:3)
7682 integer :: idir,ixa^
l,ixc^
l,ix^
d
7685 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7692 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7695 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7696 ixamin^
d=ixcmin^
d+ix^
d;
7697 ixamax^
d=ixcmax^
d+ix^
d;
7698 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7700 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7703 end subroutine get_ambipolar_electric_field
7709 integer,
intent(in) :: ixo^
l
7719 do ix^db=ixomin^db,ixomax^db\}
7721 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7722 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7723 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7724 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7725 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7726 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7729 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7730 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7731 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7732 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7775 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7776 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7777 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7779 double precision :: adummy(ixis^s,1:3)
7785 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7788 integer,
intent(in) :: ixi^
l, ixo^
l
7789 double precision,
intent(in) :: w(ixi^s,1:nw)
7790 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7791 double precision,
intent(out):: rfactor(ixi^s)
7793 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7797 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)
7799 end subroutine rfactor_from_temperature_ionization
7801 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7803 integer,
intent(in) :: ixi^
l, ixo^
l
7804 double precision,
intent(in) :: w(ixi^s,1:nw)
7805 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7806 double precision,
intent(out):: rfactor(ixi^s)
7810 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, 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.
logical, public has_equi_rho_and_p
whether split off equilibrium density and pressure
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixil, ixol)
Compute 2 times total magnetic energy.
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, protected mhd_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
logical, public, protected mhd_particles
Whether particles module is added.
integer, public, protected b
subroutine, public mhd_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
subroutine, public get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
double precision, public mhd_etah
Hall resistivity.
subroutine, public mhd_get_v(w, x, ixil, ixol, v)
Calculate v vector.
double precision, public mhd_eta_ambi
The MHD ambipolar coefficient.
logical, public, protected mhd_hydrodynamic_e
Whether hydrodynamic energy is solved instead of total energy.
subroutine, public mhd_phys_init()
logical, public, protected mhd_trac
Whether TRAC method is used.
logical, public, protected eq_state_units
type(rc_fluid), allocatable, public rc_fl
type of fluid for radiative cooling
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
integer, public, protected rho_
Index of the density (in the w array)
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
integer, public, protected tweight_
logical, public, protected mhd_ambipolar_sts
Whether Ambipolar term is implemented using supertimestepping.
procedure(sub_get_pthermal), pointer, public mhd_get_pthermal
subroutine, public mhd_ei_to_e(ixil, ixol, w, x)
Transform internal energy to total energy.
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected mhd_4th_order
MHD fourth order.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine, public mhd_get_rho(w, x, ixil, ixol, rho)
integer, public, protected psi_
Indices of the GLM psi.
logical, public mhd_equi_thermal
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_add_source(qdt, dtfactor, ixil, ixol, wct, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rotating_frame_init()
Initialize the module.
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (mhd implementation) Note: for multi-D MHD (1D MHD will use HD f...
double precision function, public get_tc_dt_hd(w, ixil, ixol, dxd, x, fl)
Get the explicit timestep for the TC (hd implementation) Note: also used in 1D MHD (or for neutrals i...
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_hd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
procedure(special_resistivity), pointer usr_special_resistivity
procedure(set_adiab), pointer usr_set_adiab
procedure(set_adiab), pointer usr_set_gamma
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
The module add viscous source terms and check time step.
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
subroutine viscosity_get_dt(w, ixil, ixol, dtnew, dxd, x)
procedure(sub_add_source), pointer, public viscosity_add_source
The data structure that contains information about a tree node/grid block.