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
133 logical,
public,
protected ::
mhd_glm = .false.
168 logical :: compactres = .false.
182 logical :: total_energy = .true.
186 logical :: gravity_energy
188 logical :: gravity_rhov = .false.
190 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
192 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
194 character(len=std_len) :: typedivbdiff =
'all'
205 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
207 integer,
intent(in) :: ixi^
l, ixo^
l
208 double precision,
intent(in) :: x(ixi^s,1:
ndim)
209 double precision,
intent(in) :: w(ixi^s,1:nw)
210 double precision,
intent(inout) :: res(ixi^s)
211 end subroutine mask_subroutine
248 subroutine mhd_read_params(files)
251 character(len=*),
intent(in) :: files(:)
267 do n = 1,
size(files)
268 open(
unitpar, file=trim(files(n)), status=
"old")
269 read(
unitpar, mhd_list,
end=111)
273 end subroutine mhd_read_params
276 subroutine mhd_write_info(fh)
278 integer,
intent(in) :: fh
281 integer,
parameter :: n_par = 1
282 double precision :: values(n_par)
283 integer,
dimension(MPI_STATUS_SIZE) :: st
284 character(len=name_len) :: names(n_par)
286 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
290 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
291 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
292 end subroutine mhd_write_info
319 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
330 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
334 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
338 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
342 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
346 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
350 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
354 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
358 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
362 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho0=F when mhd_energy=F'
366 if(
mype==0)
write(*,*)
'WARNING: set has_equi_pe0=F when mhd_energy=F'
372 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
378 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
403 phys_total_energy=total_energy
406 gravity_energy=.false.
408 gravity_energy=.true.
417 gravity_energy=.false.
423 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
428 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
437 type_divb = divb_none
440 type_divb = divb_multigrid
442 mg%operator_type = mg_laplacian
449 case (
'powel',
'powell')
450 type_divb = divb_powel
452 type_divb = divb_janhunen
454 type_divb = divb_linde
455 case (
'lindejanhunen')
456 type_divb = divb_lindejanhunen
458 type_divb = divb_lindepowel
462 type_divb = divb_lindeglm
467 call mpistop(
'Unknown divB fix')
470 allocate(start_indices(number_species),stop_indices(number_species))
477 mom(:) = var_set_momentum(
ndir)
483 e_ = var_set_energy()
492 mag(:) = var_set_bfield(
ndir)
496 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
512 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
525 te_ = var_set_auxvar(
'Te',
'Te')
534 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
742 mhd_get_rfactor=>rfactor_from_temperature_ionization
743 phys_update_temperature => mhd_update_temperature
747 mhd_get_rfactor=>rfactor_from_constant_ionization
770 phys_get_ct_velocity => mhd_get_ct_velocity
771 phys_update_faces => mhd_update_faces
773 phys_modify_wlr => mhd_modify_wlr
775 phys_boundary_adjust => mhd_boundary_adjust
784 call mhd_physical_units()
790 call mpistop(
"thermal conduction needs mhd_energy=T")
793 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
796 call mpistop(
"radiative cooling needs mhd_energy=T")
807 if(phys_internal_e)
then
809 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
811 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
815 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot_with_equi
817 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
821 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
823 tc_fl%has_equi = .true.
824 tc_fl%get_temperature_equi => mhd_get_temperature_equi
825 tc_fl%get_rho_equi => mhd_get_rho_equi
827 tc_fl%has_equi = .false.
830 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
854 rc_fl%get_var_Rfactor => mhd_get_rfactor
858 rc_fl%has_equi = .true.
859 rc_fl%get_rho_equi => mhd_get_rho_equi
860 rc_fl%get_pthermal_equi => mhd_get_pe_equi
862 rc_fl%has_equi = .false.
868 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
870 phys_te_images => mhd_te_images
886 if (particles_eta < zero) particles_eta =
mhd_eta
887 if (particles_etah < zero) particles_eta =
mhd_etah
889 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
890 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
904 phys_wider_stencil = 2
906 phys_wider_stencil = 1
914 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
926 phys_wider_stencil = 2
928 phys_wider_stencil = 1
942 subroutine mhd_te_images
947 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
949 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
951 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
953 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
956 call mpistop(
"Error in synthesize emission: Unknown convert_type")
958 end subroutine mhd_te_images
964 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
968 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
969 double precision,
intent(in) :: x(ixi^s,1:
ndim)
970 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
971 double precision,
intent(in) :: my_dt
972 logical,
intent(in) :: fix_conserve_at_step
974 end subroutine mhd_sts_set_source_tc_mhd
976 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
983 integer,
intent(in) :: ixi^
l, ixo^
l
984 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
985 double precision,
intent(in) :: w(ixi^s,1:nw)
986 double precision :: dtnew
989 end function mhd_get_tc_dt_mhd
991 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
994 integer,
intent(in) :: ixi^
l,ixo^
l
995 double precision,
intent(inout) :: w(ixi^s,1:nw)
996 double precision,
intent(in) :: x(ixi^s,1:
ndim)
997 integer,
intent(in) :: step
998 character(len=140) :: error_msg
1000 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1001 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1002 end subroutine mhd_tc_handle_small_e
1005 subroutine tc_params_read_mhd(fl)
1007 type(tc_fluid),
intent(inout) :: fl
1009 double precision :: tc_k_para=0d0
1010 double precision :: tc_k_perp=0d0
1013 logical :: tc_perpendicular=.false.
1014 logical :: tc_saturate=.false.
1015 character(len=std_len) :: tc_slope_limiter=
"MC"
1017 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1021 read(
unitpar, tc_list,
end=111)
1025 fl%tc_perpendicular = tc_perpendicular
1026 fl%tc_saturate = tc_saturate
1027 fl%tc_k_para = tc_k_para
1028 fl%tc_k_perp = tc_k_perp
1029 select case(tc_slope_limiter)
1031 fl%tc_slope_limiter = 0
1034 fl%tc_slope_limiter = 1
1037 fl%tc_slope_limiter = 2
1040 fl%tc_slope_limiter = 3
1043 fl%tc_slope_limiter = 4
1045 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1047 end subroutine tc_params_read_mhd
1051 subroutine rc_params_read(fl)
1054 type(rc_fluid),
intent(inout) :: fl
1056 double precision :: cfrac=0.1d0
1059 double precision :: rad_cut_hgt=0.5d0
1060 double precision :: rad_cut_dey=0.15d0
1063 integer :: ncool = 4000
1065 logical :: tfix=.false.
1067 logical :: rc_split=.false.
1068 logical :: rad_cut=.false.
1070 character(len=std_len) :: coolcurve=
'JCcorona'
1072 character(len=std_len) :: coolmethod=
'exact'
1074 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1078 read(
unitpar, rc_list,
end=111)
1083 fl%coolcurve=coolcurve
1084 fl%coolmethod=coolmethod
1087 fl%rc_split=rc_split
1090 fl%rad_cut_hgt=rad_cut_hgt
1091 fl%rad_cut_dey=rad_cut_dey
1092 end subroutine rc_params_read
1096 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1099 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1100 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1102 double precision :: delx(ixi^s,1:
ndim)
1103 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1104 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1110 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1114 hxo^
l=ixo^
l-
kr(idims,^
d);
1120 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1123 xshift^
d=half*(one-
kr(^
d,idims));
1130 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1134 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1137 end subroutine set_equi_vars_grid_faces
1140 subroutine set_equi_vars_grid(igrid)
1144 integer,
intent(in) :: igrid
1150 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1152 end subroutine set_equi_vars_grid
1155 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1157 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1158 double precision,
intent(in) :: w(ixi^s, 1:nw)
1159 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1160 double precision :: wnew(ixo^s, 1:nwc)
1167 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1173 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1177 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1181 if(
b0field .and. total_energy)
then
1182 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1183 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1187 end function convert_vars_splitting
1189 subroutine mhd_check_params
1197 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1198 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1202 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1203 inv_gamma_1=1.d0/gamma_1
1208 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1213 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1218 end subroutine mhd_check_params
1220 subroutine mhd_physical_units()
1222 double precision :: mp,kb,miu0,c_lightspeed
1223 double precision :: a,
b
1234 c_lightspeed=const_c
1375 end subroutine mhd_physical_units
1377 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1380 logical,
intent(in) :: primitive
1381 logical,
intent(inout) :: flag(ixi^s,1:nw)
1382 integer,
intent(in) :: ixi^
l, ixo^
l
1383 double precision,
intent(in) :: w(ixi^s,nw)
1385 double precision :: tmp,b2,
b(ixo^s,1:
ndir)
1386 double precision :: v(ixo^s,1:
ndir),gamma2,inv_rho
1397 {
do ix^db=ixomin^db,ixomax^db \}
1398 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1401 {
do ix^db=ixomin^db,ixomax^db \}
1402 b2=(^
c&w(ix^d,
b^
c_)**2+)
1403 if(b2>smalldouble)
then
1408 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1409 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1410 inv_rho = 1d0/w(ix^d,
rho_)
1412 b2=b2*inv_rho*inv_squared_c
1414 gamma2=1.d0/(1.d0+b2)
1416 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1419 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1420 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1421 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1426 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1432 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1433 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1434 if(tmp<small_e) flag(ix^d,
e_)=.true.
1440 end subroutine mhd_check_w_semirelati
1442 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1445 logical,
intent(in) :: primitive
1446 integer,
intent(in) :: ixi^
l, ixo^
l
1447 double precision,
intent(in) :: w(ixi^s,nw)
1448 logical,
intent(inout) :: flag(ixi^s,1:nw)
1453 {
do ix^db=ixomin^db,ixomax^db\}
1459 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1463 end subroutine mhd_check_w_origin
1465 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1468 logical,
intent(in) :: primitive
1469 integer,
intent(in) :: ixi^
l, ixo^
l
1470 double precision,
intent(in) :: w(ixi^s,nw)
1471 logical,
intent(inout) :: flag(ixi^s,1:nw)
1473 double precision :: tmp
1477 {
do ix^db=ixomin^db,ixomax^db\}
1483 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1488 end subroutine mhd_check_w_split
1490 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1493 logical,
intent(in) :: primitive
1494 integer,
intent(in) :: ixi^
l, ixo^
l
1495 double precision,
intent(in) :: w(ixi^s,nw)
1496 logical,
intent(inout) :: flag(ixi^s,1:nw)
1501 {
do ix^db=ixomin^db,ixomax^db\}
1505 end subroutine mhd_check_w_noe
1507 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1510 logical,
intent(in) :: primitive
1511 integer,
intent(in) :: ixi^
l, ixo^
l
1512 double precision,
intent(in) :: w(ixi^s,nw)
1513 logical,
intent(inout) :: flag(ixi^s,1:nw)
1518 {
do ix^db=ixomin^db,ixomax^db\}
1523 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1527 end subroutine mhd_check_w_inte
1529 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1532 logical,
intent(in) :: primitive
1533 integer,
intent(in) :: ixi^
l, ixo^
l
1534 double precision,
intent(in) :: w(ixi^s,nw)
1535 logical,
intent(inout) :: flag(ixi^s,1:nw)
1540 {
do ix^db=ixomin^db,ixomax^db\}
1545 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1549 end subroutine mhd_check_w_hde
1552 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1554 integer,
intent(in) :: ixi^
l, ixo^
l
1555 double precision,
intent(inout) :: w(ixi^s, nw)
1556 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1560 {
do ix^db=ixomin^db,ixomax^db\}
1562 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1564 +(^
c&w(ix^
d,
b^
c_)**2+))
1569 end subroutine mhd_to_conserved_origin
1572 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1574 integer,
intent(in) :: ixi^
l, ixo^
l
1575 double precision,
intent(inout) :: w(ixi^s, nw)
1576 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1580 {
do ix^db=ixomin^db,ixomax^db\}
1585 end subroutine mhd_to_conserved_origin_noe
1588 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1590 integer,
intent(in) :: ixi^
l, ixo^
l
1591 double precision,
intent(inout) :: w(ixi^s, nw)
1592 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1596 {
do ix^db=ixomin^db,ixomax^db\}
1598 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1604 end subroutine mhd_to_conserved_hde
1607 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1609 integer,
intent(in) :: ixi^
l, ixo^
l
1610 double precision,
intent(inout) :: w(ixi^s, nw)
1611 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1615 {
do ix^db=ixomin^db,ixomax^db\}
1617 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1622 end subroutine mhd_to_conserved_inte
1625 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1627 integer,
intent(in) :: ixi^
l, ixo^
l
1628 double precision,
intent(inout) :: w(ixi^s, nw)
1629 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1631 double precision :: rho
1634 {
do ix^db=ixomin^db,ixomax^db\}
1637 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1638 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1639 +(^
c&w(ix^
d,
b^
c_)**2+))
1644 end subroutine mhd_to_conserved_split_rho
1647 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1649 integer,
intent(in) :: ixi^
l, ixo^
l
1650 double precision,
intent(inout) :: w(ixi^s, nw)
1651 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1653 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1656 {
do ix^db=ixomin^db,ixomax^db\}
1658 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1659 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1660 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1661 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1662 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1663 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1668 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1669 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1670 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1678 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1682 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1684 +(^
c&w(ix^
d,
b^
c_)**2+)&
1685 +(^
c&e(ix^
d,^
c)**2+)*inv_squared_c)
1693 end subroutine mhd_to_conserved_semirelati
1695 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1697 integer,
intent(in) :: ixi^
l, ixo^
l
1698 double precision,
intent(inout) :: w(ixi^s, nw)
1699 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1701 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1704 {
do ix^db=ixomin^db,ixomax^db\}
1706 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1707 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1708 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1709 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1710 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1711 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1716 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1717 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1718 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1728 end subroutine mhd_to_conserved_semirelati_noe
1731 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1733 integer,
intent(in) :: ixi^
l, ixo^
l
1734 double precision,
intent(inout) :: w(ixi^s, nw)
1735 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1737 double precision :: inv_rho
1742 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1745 {
do ix^db=ixomin^db,ixomax^db\}
1746 inv_rho = 1.d0/w(ix^
d,
rho_)
1750 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1752 +(^
c&w(ix^
d,
b^
c_)**2+)))
1755 end subroutine mhd_to_primitive_origin
1758 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1760 integer,
intent(in) :: ixi^
l, ixo^
l
1761 double precision,
intent(inout) :: w(ixi^s, nw)
1762 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1764 double precision :: inv_rho
1769 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1772 {
do ix^db=ixomin^db,ixomax^db\}
1773 inv_rho = 1.d0/w(ix^
d,
rho_)
1778 end subroutine mhd_to_primitive_origin_noe
1781 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1783 integer,
intent(in) :: ixi^
l, ixo^
l
1784 double precision,
intent(inout) :: w(ixi^s, nw)
1785 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1787 double precision :: inv_rho
1792 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1795 {
do ix^db=ixomin^db,ixomax^db\}
1796 inv_rho = 1d0/w(ix^
d,
rho_)
1803 end subroutine mhd_to_primitive_hde
1806 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1808 integer,
intent(in) :: ixi^
l, ixo^
l
1809 double precision,
intent(inout) :: w(ixi^s, nw)
1810 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1812 double precision :: inv_rho
1817 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1820 {
do ix^db=ixomin^db,ixomax^db\}
1822 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1824 inv_rho = 1.d0/w(ix^
d,
rho_)
1828 end subroutine mhd_to_primitive_inte
1831 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1833 integer,
intent(in) :: ixi^
l, ixo^
l
1834 double precision,
intent(inout) :: w(ixi^s, nw)
1835 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1837 double precision :: inv_rho
1842 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1845 {
do ix^db=ixomin^db,ixomax^db\}
1850 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1852 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1855 end subroutine mhd_to_primitive_split_rho
1858 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1860 integer,
intent(in) :: ixi^
l, ixo^
l
1861 double precision,
intent(inout) :: w(ixi^s, nw)
1862 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1864 double precision ::
b(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
1869 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1872 {
do ix^db=ixomin^db,ixomax^db\}
1873 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1874 if(b2>smalldouble)
then
1882 inv_rho=1.d0/w(ix^
d,
rho_)
1884 b2=b2*inv_rho*inv_squared_c
1886 gamma2=1.d0/(1.d0+b2)
1888 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1892 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1896 b(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1897 b(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1898 b(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1902 b(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1908 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1910 +(^
c&w(ix^
d,
b^
c_)**2+)&
1911 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
1915 end subroutine mhd_to_primitive_semirelati
1918 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1920 integer,
intent(in) :: ixi^
l, ixo^
l
1921 double precision,
intent(inout) :: w(ixi^s, nw)
1922 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1924 double precision ::
b(ixo^s,1:
ndir),tmp,b2,gamma2,inv_rho
1925 integer :: ix^
d, idir
1929 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1932 {
do ix^db=ixomin^db,ixomax^db\}
1933 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1934 if(b2>smalldouble)
then
1942 inv_rho=1.d0/w(ix^
d,
rho_)
1944 b2=b2*inv_rho*inv_squared_c
1946 gamma2=1.d0/(1.d0+b2)
1948 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1951 end subroutine mhd_to_primitive_semirelati_noe
1956 integer,
intent(in) :: ixi^
l, ixo^
l
1957 double precision,
intent(inout) :: w(ixi^s, nw)
1958 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1963 {
do ix^db=ixomin^db,ixomax^db\}
1966 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1968 +(^
c&w(ix^
d,
b^
c_)**2+))
1971 {
do ix^db=ixomin^db,ixomax^db\}
1973 w(ix^d,
e_)=w(ix^d,
e_)&
1974 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1975 +(^
c&w(ix^d,
b^
c_)**2+))
1982 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
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)
1990 {
do ix^db=ixomin^db,ixomax^db\}
1996 end subroutine mhd_ei_to_e_hde
1999 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
2001 integer,
intent(in) :: ixi^
l, ixo^
l
2002 double precision,
intent(inout) :: w(ixi^s, nw)
2003 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2005 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2006 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2008 end subroutine mhd_ei_to_e_semirelati
2013 integer,
intent(in) :: ixi^
l, ixo^
l
2014 double precision,
intent(inout) :: w(ixi^s, nw)
2015 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2020 {
do ix^db=ixomin^db,ixomax^db\}
2023 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2025 +(^
c&w(ix^
d,
b^
c_)**2+))
2028 {
do ix^db=ixomin^db,ixomax^db\}
2030 w(ix^d,
e_)=w(ix^d,
e_)&
2031 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2032 +(^
c&w(ix^d,
b^
c_)**2+))
2036 if(fix_small_values)
then
2037 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2043 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2045 integer,
intent(in) :: ixi^
l, ixo^
l
2046 double precision,
intent(inout) :: w(ixi^s, nw)
2047 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2051 {
do ix^db=ixomin^db,ixomax^db\}
2057 if(fix_small_values)
then
2058 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2061 end subroutine mhd_e_to_ei_hde
2064 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2066 integer,
intent(in) :: ixi^
l, ixo^
l
2067 double precision,
intent(inout) :: w(ixi^s, nw)
2068 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2070 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2071 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2073 end subroutine mhd_e_to_ei_semirelati
2075 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2078 logical,
intent(in) :: primitive
2079 integer,
intent(in) :: ixi^
l,ixo^
l
2080 double precision,
intent(inout) :: w(ixi^s,1:nw)
2081 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2082 character(len=*),
intent(in) :: subname
2084 double precision ::
b(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2085 double precision :: tmp, b2, gamma2, inv_rho
2087 logical :: flag(ixi^s,1:nw)
2096 {
do ix^db=ixomin^db,ixomax^db\}
2097 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2098 if(b2>smalldouble)
then
2105 inv_rho=1.d0/w(ix^
d,
rho_)
2107 b2=b2*inv_rho*inv_squared_c
2109 gamma2=1.d0/(1.d0+b2)
2111 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
2114 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2115 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2116 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2120 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2126 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2127 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2128 +(^
c&w(ix^
d,
b^
c_)**2+)&
2129 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
2136 select case (small_values_method)
2138 {
do ix^db=ixomin^db,ixomax^db\}
2139 if(flag(ix^d,
rho_))
then
2140 w(ix^d,
rho_) = small_density
2141 ^
c&w(ix^d,
m^
c_)=0.d0\
2145 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2147 if(flag(ix^d,
e_))
then
2148 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2149 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2156 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2159 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2161 w(ixo^s,
e_)=pressure(ixo^s)
2162 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2163 {
do ix^db=ixomin^db,ixomax^db\}
2164 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2165 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2170 if(.not.primitive)
then
2172 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2173 w(ixo^s,
e_)=pressure(ixo^s)
2175 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2179 end subroutine mhd_handle_small_values_semirelati
2181 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2184 logical,
intent(in) :: primitive
2185 integer,
intent(in) :: ixi^
l,ixo^
l
2186 double precision,
intent(inout) :: w(ixi^s,1:nw)
2187 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2188 character(len=*),
intent(in) :: subname
2191 logical :: flag(ixi^s,1:nw)
2193 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2198 {
do ix^db=ixomin^db,ixomax^db\}
2202 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2214 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2216 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2219 {
do ix^db=iximin^db,iximax^db\}
2220 w(ix^d,
e_)=w(ix^d,
e_)&
2221 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2223 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2225 {
do ix^db=iximin^db,iximax^db\}
2226 w(ix^d,
e_)=w(ix^d,
e_)&
2227 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2231 if(.not.primitive)
then
2233 {
do ix^db=ixomin^db,ixomax^db\}
2235 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2236 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2239 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2243 end subroutine mhd_handle_small_values_origin
2245 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2248 logical,
intent(in) :: primitive
2249 integer,
intent(in) :: ixi^
l,ixo^
l
2250 double precision,
intent(inout) :: w(ixi^s,1:nw)
2251 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2252 character(len=*),
intent(in) :: subname
2254 double precision :: rho
2256 logical :: flag(ixi^s,1:nw)
2258 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2263 {
do ix^db=ixomin^db,ixomax^db\}
2268 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2275 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2281 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2283 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2286 {
do ix^db=iximin^db,iximax^db\}
2288 w(ix^d,
e_)=w(ix^d,
e_)&
2289 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2291 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2293 {
do ix^db=iximin^db,iximax^db\}
2295 w(ix^d,
e_)=w(ix^d,
e_)&
2296 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2300 if(.not.primitive)
then
2302 {
do ix^db=ixomin^db,ixomax^db\}
2304 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2305 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2306 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2309 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2313 end subroutine mhd_handle_small_values_split
2315 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2318 logical,
intent(in) :: primitive
2319 integer,
intent(in) :: ixi^
l,ixo^
l
2320 double precision,
intent(inout) :: w(ixi^s,1:nw)
2321 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2322 character(len=*),
intent(in) :: subname
2325 logical :: flag(ixi^s,1:nw)
2327 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2332 {
do ix^db=ixomin^db,ixomax^db\}
2333 if(flag(ix^
d,
rho_))
then
2335 ^
c&w(ix^
d,
m^
c_)=0.d0\
2340 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2345 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2347 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2349 if(.not.primitive)
then
2351 {
do ix^db=ixomin^db,ixomax^db\}
2353 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2356 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2360 end subroutine mhd_handle_small_values_inte
2362 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2365 logical,
intent(in) :: primitive
2366 integer,
intent(in) :: ixi^
l,ixo^
l
2367 double precision,
intent(inout) :: w(ixi^s,1:nw)
2368 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2369 character(len=*),
intent(in) :: subname
2372 logical :: flag(ixi^s,1:nw)
2374 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2379 {
do ix^db=ixomin^db,ixomax^db\}
2383 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2389 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2391 if(.not.primitive)
then
2393 {
do ix^db=ixomin^db,ixomax^db\}
2397 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2401 end subroutine mhd_handle_small_values_noe
2403 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2406 logical,
intent(in) :: primitive
2407 integer,
intent(in) :: ixi^
l,ixo^
l
2408 double precision,
intent(inout) :: w(ixi^s,1:nw)
2409 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2410 character(len=*),
intent(in) :: subname
2413 logical :: flag(ixi^s,1:nw)
2415 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2420 {
do ix^db=ixomin^db,ixomax^db\}
2421 if(flag(ix^
d,
rho_))
then
2423 ^
c&w(ix^
d,
m^
c_)=0.d0\
2428 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2433 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2435 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2437 if(.not.primitive)
then
2439 {
do ix^db=ixomin^db,ixomax^db\}
2441 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2444 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2448 end subroutine mhd_handle_small_values_hde
2454 integer,
intent(in) :: ixi^
l, ixo^
l
2455 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2456 double precision,
intent(out) :: v(ixi^s,
ndir)
2458 double precision :: rho(ixi^s)
2463 rho(ixo^s)=1.d0/rho(ixo^s)
2466 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2472 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2475 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2476 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2477 double precision,
intent(inout) :: cmax(ixi^s)
2479 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2485 {
do ix^db=ixomin^db,ixomax^db \}
2496 cfast2=b2*inv_rho+cmax(ix^
d)
2497 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2498 if(avmincs2<zero) avmincs2=zero
2499 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2503 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2505 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2508 {
do ix^db=ixomin^db,ixomax^db \}
2518 b2=(^
c&w(ix^d,
b^
c_)**2+)
2519 cfast2=b2*inv_rho+cmax(ix^d)
2520 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2521 if(avmincs2<zero) avmincs2=zero
2522 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2526 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2528 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2532 end subroutine mhd_get_cmax_origin
2535 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2538 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2539 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2540 double precision,
intent(inout) :: cmax(ixi^s)
2542 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2548 {
do ix^db=ixomin^db,ixomax^db \}
2559 cfast2=b2*inv_rho+cmax(ix^
d)
2560 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2561 if(avmincs2<zero) avmincs2=zero
2562 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2566 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2568 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2571 {
do ix^db=ixomin^db,ixomax^db \}
2581 b2=(^
c&w(ix^d,
b^
c_)**2+)
2582 cfast2=b2*inv_rho+cmax(ix^d)
2583 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2584 if(avmincs2<zero) avmincs2=zero
2585 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2589 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2591 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2595 end subroutine mhd_get_cmax_origin_noe
2598 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2601 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2602 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2603 double precision,
intent(inout):: cmax(ixi^s)
2605 double precision :: csound, avmincs2, idim_alfven_speed2
2606 double precision :: inv_rho, alfven_speed2, gamma2
2609 {
do ix^db=ixomin^db,ixomax^db \}
2610 inv_rho=1.d0/w(ix^
d,
rho_)
2611 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2612 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2613 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2616 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2619 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2620 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2621 if(avmincs2<zero) avmincs2=zero
2623 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2624 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2627 end subroutine mhd_get_cmax_semirelati
2630 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2633 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2634 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2635 double precision,
intent(inout):: cmax(ixi^s)
2637 double precision :: csound, avmincs2, idim_alfven_speed2
2638 double precision :: inv_rho, alfven_speed2, gamma2
2641 {
do ix^db=ixomin^db,ixomax^db \}
2642 inv_rho=1.d0/w(ix^
d,
rho_)
2643 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2644 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2645 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2647 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2650 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2651 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2652 if(avmincs2<zero) avmincs2=zero
2654 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2655 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2658 end subroutine mhd_get_cmax_semirelati_noe
2660 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2663 integer,
intent(in) :: ixi^
l, ixo^
l
2664 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2665 double precision,
intent(inout) :: a2max(
ndim)
2666 double precision :: a2(ixi^s,
ndim,nw)
2667 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2672 hxo^
l=ixo^
l-
kr(i,^
d);
2673 gxo^
l=hxo^
l-
kr(i,^
d);
2674 jxo^
l=ixo^
l+
kr(i,^
d);
2675 kxo^
l=jxo^
l+
kr(i,^
d);
2676 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2677 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2678 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2680 end subroutine mhd_get_a2max
2683 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2686 integer,
intent(in) :: ixi^
l,ixo^
l
2687 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2689 double precision,
intent(inout) :: w(ixi^s,1:nw)
2690 double precision,
intent(out) :: tco_local,tmax_local
2692 double precision,
parameter :: trac_delta=0.25d0
2693 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
2694 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
2695 double precision,
dimension(ixI^S,1:ndim) :: gradt
2696 double precision :: bdir(
ndim)
2697 double precision :: ltrc,ltrp,altr(ixi^s)
2698 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d,ix^
d
2699 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2700 logical :: lrlt(ixi^s)
2703 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2705 call mhd_get_rfactor(w,x,ixi^
l,ixi^
l,te)
2706 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2709 tmax_local=maxval(te(ixo^s))
2719 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2721 where(lts(ixo^s) > trac_delta)
2724 if(any(lrlt(ixo^s)))
then
2725 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2736 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2737 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2738 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2739 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2741 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2753 gradt(ixo^s,idims)=tmp1(ixo^s)
2757 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
2759 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2765 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2768 if(sum(bdir(:)**2) .gt. zero)
then
2769 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2771 block%special_values(3:ndim+2)=bdir(1:ndim)
2773 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2774 where(tmp1(ixo^s)/=0.d0)
2775 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2777 tmp1(ixo^s)=bigdouble
2781 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2784 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2786 if(slab_uniform)
then
2787 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2789 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2792 where(lts(ixo^s) > trac_delta)
2795 if(any(lrlt(ixo^s)))
then
2796 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2798 block%special_values(1)=zero
2800 block%special_values(2)=tmax_local
2819 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2820 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2821 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2824 {
do ix^db=ixpmin^db,ixpmax^db\}
2826 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))+block%B0(ix^d,^
c,0)\
2828 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))\
2830 tmp1(ix^d)=1.d0/(dsqrt(^
c&bunitvec(ix^d,^
c)**2+)+smalldouble)
2832 ^d&bunitvec({ix^d},^d)=bunitvec({ix^d},^d)*tmp1({ix^d})\
2834 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec({ix^d},^d)+)/te(ix^d)
2836 if(slab_uniform)
then
2837 lts(ix^d)=min(^d&dxlevel(^d))*lts(ix^d)
2839 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2841 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2847 hxo^l=ixp^l-kr(idims,^d);
2848 jxo^l=ixp^l+kr(idims,^d);
2850 altr(ixp^s)=0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2852 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2855 block%wextra(ixp^s,
tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
2859 call mpistop(
"unknown mhd_trac_type")
2862 end subroutine mhd_get_tcutoff
2865 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2868 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2869 double precision,
intent(in) :: wprim(ixi^s, nw)
2870 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2871 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2873 double precision :: csound(ixi^s,
ndim)
2874 double precision,
allocatable :: tmp(:^
d&)
2875 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2879 allocate(tmp(ixa^s))
2881 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2882 csound(ixa^s,id)=tmp(ixa^s)
2885 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2886 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2887 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2888 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))
2892 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2893 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2894 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)))
2895 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2896 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2897 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)))
2902 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2903 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2904 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)))
2905 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2906 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2907 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)))
2911 end subroutine mhd_get_h_speed
2914 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2917 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2918 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2919 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2920 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2921 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
2922 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
2923 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2925 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
2926 double precision :: umean, dmean, tmp1, tmp2, tmp3
2933 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2934 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2935 if(
present(cmin))
then
2936 {
do ix^db=ixomin^db,ixomax^db\}
2937 tmp1=sqrt(wlp(ix^
d,
rho_))
2938 tmp2=sqrt(wrp(ix^
d,
rho_))
2939 tmp3=1.d0/(tmp1+tmp2)
2940 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
2941 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
2942 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
2943 cmin(ix^
d,1)=umean-dmean
2944 cmax(ix^
d,1)=umean+dmean
2946 if(h_correction)
then
2947 {
do ix^db=ixomin^db,ixomax^db\}
2948 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2949 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2953 {
do ix^db=ixomin^db,ixomax^db\}
2954 tmp1=sqrt(wlp(ix^d,
rho_))
2955 tmp2=sqrt(wrp(ix^d,
rho_))
2956 tmp3=1.d0/(tmp1+tmp2)
2957 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2958 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2959 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2960 cmax(ix^d,1)=abs(umean)+dmean
2964 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2965 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2966 if(
present(cmin))
then
2967 {
do ix^db=ixomin^db,ixomax^db\}
2968 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
2969 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
2971 if(h_correction)
then
2972 {
do ix^db=ixomin^db,ixomax^db\}
2973 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2974 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2978 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
2982 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2983 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2984 if(
present(cmin))
then
2985 {
do ix^db=ixomin^db,ixomax^db\}
2986 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2987 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2988 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2990 if(h_correction)
then
2991 {
do ix^db=ixomin^db,ixomax^db\}
2992 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2993 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2997 {
do ix^db=ixomin^db,ixomax^db\}
2998 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2999 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3004 end subroutine mhd_get_cbounds
3007 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3010 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3011 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3012 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3013 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3014 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3015 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3016 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3018 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3023 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3024 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3026 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3027 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3029 if(
present(cmin))
then
3030 {
do ix^db=ixomin^db,ixomax^db\}
3031 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3032 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3033 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3036 {
do ix^db=ixomin^db,ixomax^db\}
3037 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3038 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3042 end subroutine mhd_get_cbounds_semirelati
3045 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3048 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3049 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3050 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3051 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3052 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3053 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3054 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3056 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3057 double precision :: umean, dmean, tmp1, tmp2, tmp3
3064 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3065 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3066 if(
present(cmin))
then
3067 {
do ix^db=ixomin^db,ixomax^db\}
3070 tmp3=1.d0/(tmp1+tmp2)
3071 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3072 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3073 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3074 cmin(ix^
d,1)=umean-dmean
3075 cmax(ix^
d,1)=umean+dmean
3077 if(h_correction)
then
3078 {
do ix^db=ixomin^db,ixomax^db\}
3079 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3080 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3084 {
do ix^db=ixomin^db,ixomax^db\}
3087 tmp3=1.d0/(tmp1+tmp2)
3088 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3089 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3090 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3091 cmax(ix^d,1)=abs(umean)+dmean
3095 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3096 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3097 if(
present(cmin))
then
3098 {
do ix^db=ixomin^db,ixomax^db\}
3099 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3100 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3102 if(h_correction)
then
3103 {
do ix^db=ixomin^db,ixomax^db\}
3104 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3105 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3109 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3113 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3114 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3115 if(
present(cmin))
then
3116 {
do ix^db=ixomin^db,ixomax^db\}
3117 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3118 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3119 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3121 if(h_correction)
then
3122 {
do ix^db=ixomin^db,ixomax^db\}
3123 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3124 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3128 {
do ix^db=ixomin^db,ixomax^db\}
3129 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3130 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3135 end subroutine mhd_get_cbounds_split_rho
3138 subroutine mhd_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3141 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3142 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3143 double precision,
intent(in) :: cmax(ixi^s)
3144 double precision,
intent(in),
optional :: cmin(ixi^s)
3145 type(ct_velocity),
intent(inout):: vcts
3147 integer :: idime,idimn
3153 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3155 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3157 if(.not.
allocated(vcts%vbarC))
then
3158 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3159 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3162 if(
present(cmin))
then
3163 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3164 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3166 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3167 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3170 idimn=mod(idim,
ndir)+1
3171 idime=mod(idim+1,
ndir)+1
3173 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3174 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3175 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3176 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3177 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3179 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3180 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3181 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3182 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3183 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3185 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
3188 end subroutine mhd_get_ct_velocity
3191 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3194 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3195 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3196 double precision,
intent(out):: csound(ixo^s)
3198 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3205 {
do ix^db=ixomin^db,ixomax^db \}
3206 inv_rho=1.d0/w(ix^
d,
rho_)
3213 cfast2=b2*inv_rho+csound(ix^
d)
3214 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3216 if(avmincs2<zero) avmincs2=zero
3217 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3219 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3223 {
do ix^db=ixomin^db,ixomax^db \}
3224 inv_rho=1.d0/w(ix^d,
rho_)
3230 b2=(^
c&w(ix^d,
b^
c_)**2+)
3231 cfast2=b2*inv_rho+csound(ix^d)
3232 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3233 if(avmincs2<zero) avmincs2=zero
3234 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3236 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3241 end subroutine mhd_get_csound_prim
3244 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3247 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3248 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3249 double precision,
intent(out):: csound(ixo^s)
3251 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3258 {
do ix^db=ixomin^db,ixomax^db \}
3265 cfast2=b2*inv_rho+csound(ix^
d)
3266 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3268 if(avmincs2<zero) avmincs2=zero
3269 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3271 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3275 {
do ix^db=ixomin^db,ixomax^db \}
3281 b2=(^
c&w(ix^d,
b^
c_)**2+)
3282 cfast2=b2*inv_rho+csound(ix^d)
3283 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3284 if(avmincs2<zero) avmincs2=zero
3285 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3287 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3292 end subroutine mhd_get_csound_prim_split
3295 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3298 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3300 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3301 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3303 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3306 {
do ix^db=ixomin^db,ixomax^db\}
3307 inv_rho = 1.d0/w(ix^
d,
rho_)
3310 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3311 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3312 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3313 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3316 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3317 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3318 if(avmincs2<zero) avmincs2=zero
3320 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3323 end subroutine mhd_get_csound_semirelati
3326 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3329 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3331 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3332 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3334 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3337 {
do ix^db=ixomin^db,ixomax^db\}
3338 inv_rho = 1.d0/w(ix^
d,
rho_)
3341 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3342 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3343 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3344 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3347 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3348 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3349 if(avmincs2<zero) avmincs2=zero
3351 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3354 end subroutine mhd_get_csound_semirelati_noe
3357 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3360 integer,
intent(in) :: ixi^
l, ixo^
l
3361 double precision,
intent(in) :: w(ixi^s,nw)
3362 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3363 double precision,
intent(out):: pth(ixi^s)
3371 end subroutine mhd_get_pthermal_noe
3374 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3378 integer,
intent(in) :: ixi^
l, ixo^
l
3379 double precision,
intent(in) :: w(ixi^s,nw)
3380 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3381 double precision,
intent(out):: pth(ixi^s)
3385 {
do ix^db= ixomin^db,ixomax^db\}
3389 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3394 if(check_small_values.and..not.fix_small_values)
then
3395 {
do ix^db= ixomin^db,ixomax^db\}
3396 if(pth(ix^d)<small_pressure)
then
3397 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3398 " encountered when call mhd_get_pthermal_inte"
3399 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3400 write(*,*)
"Location: ", x(ix^d,:)
3401 write(*,*)
"Cell number: ", ix^d
3403 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3406 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3407 write(*,*)
"Saving status at the previous time step"
3413 end subroutine mhd_get_pthermal_inte
3416 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3420 integer,
intent(in) :: ixi^
l, ixo^
l
3421 double precision,
intent(in) :: w(ixi^s,nw)
3422 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3423 double precision,
intent(out):: pth(ixi^s)
3427 {
do ix^db=ixomin^db,ixomax^db\}
3432 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3433 +(^
c&w(ix^
d,
b^
c_)**2+)))
3438 if(check_small_values.and..not.fix_small_values)
then
3439 {
do ix^db=ixomin^db,ixomax^db\}
3440 if(pth(ix^d)<small_pressure)
then
3441 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3442 " encountered when call mhd_get_pthermal"
3443 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3444 write(*,*)
"Location: ", x(ix^d,:)
3445 write(*,*)
"Cell number: ", ix^d
3447 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3450 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3451 write(*,*)
"Saving status at the previous time step"
3457 end subroutine mhd_get_pthermal_origin
3460 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3464 integer,
intent(in) :: ixi^
l, ixo^
l
3465 double precision,
intent(in) :: w(ixi^s,nw)
3466 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3467 double precision,
intent(out):: pth(ixi^s)
3469 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3472 {
do ix^db=ixomin^db,ixomax^db\}
3473 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3474 if(b2>smalldouble)
then
3482 inv_rho=1.d0/w(ix^
d,
rho_)
3484 b2=b2*inv_rho*inv_squared_c
3486 gamma2=1.d0/(1.d0+b2)
3488 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3492 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3493 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3494 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3498 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3504 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3505 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3506 +(^
c&w(ix^
d,
b^
c_)**2+)&
3507 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3511 if(check_small_values.and..not.fix_small_values)
then
3512 {
do ix^db=ixomin^db,ixomax^db\}
3513 if(pth(ix^d)<small_pressure)
then
3514 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3515 " encountered when call mhd_get_pthermal_semirelati"
3516 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3517 write(*,*)
"Location: ", x(ix^d,:)
3518 write(*,*)
"Cell number: ", ix^d
3520 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3523 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3524 write(*,*)
"Saving status at the previous time step"
3530 end subroutine mhd_get_pthermal_semirelati
3533 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3537 integer,
intent(in) :: ixi^
l, ixo^
l
3538 double precision,
intent(in) :: w(ixi^s,nw)
3539 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3540 double precision,
intent(out):: pth(ixi^s)
3544 {
do ix^db= ixomin^db,ixomax^db\}
3545 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3548 if(check_small_values.and..not.fix_small_values)
then
3549 {
do ix^db= ixomin^db,ixomax^db\}
3550 if(pth(ix^d)<small_pressure)
then
3551 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3552 " encountered when call mhd_get_pthermal_hde"
3553 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3554 write(*,*)
"Location: ", x(ix^d,:)
3555 write(*,*)
"Cell number: ", ix^d
3557 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3560 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3561 write(*,*)
"Saving status at the previous time step"
3567 end subroutine mhd_get_pthermal_hde
3570 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3572 integer,
intent(in) :: ixi^
l, ixo^
l
3573 double precision,
intent(in) :: w(ixi^s, 1:nw)
3574 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3575 double precision,
intent(out):: res(ixi^s)
3576 res(ixo^s) = w(ixo^s,
te_)
3577 end subroutine mhd_get_temperature_from_te
3580 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3582 integer,
intent(in) :: ixi^
l, ixo^
l
3583 double precision,
intent(in) :: w(ixi^s, 1:nw)
3584 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3585 double precision,
intent(out):: res(ixi^s)
3587 double precision :: r(ixi^s)
3589 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3590 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3591 end subroutine mhd_get_temperature_from_eint
3594 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3596 integer,
intent(in) :: ixi^
l, ixo^
l
3597 double precision,
intent(in) :: w(ixi^s, 1:nw)
3598 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3599 double precision,
intent(out):: res(ixi^s)
3601 double precision :: r(ixi^s)
3603 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3605 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3607 end subroutine mhd_get_temperature_from_etot
3609 subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
3611 integer,
intent(in) :: ixi^
l, ixo^
l
3612 double precision,
intent(in) :: w(ixi^s, 1:nw)
3613 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3614 double precision,
intent(out):: res(ixi^s)
3616 double precision :: r(ixi^s)
3618 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3622 end subroutine mhd_get_temperature_from_etot_with_equi
3624 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3626 integer,
intent(in) :: ixi^
l, ixo^
l
3627 double precision,
intent(in) :: w(ixi^s, 1:nw)
3628 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3629 double precision,
intent(out):: res(ixi^s)
3631 double precision :: r(ixi^s)
3633 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3637 end subroutine mhd_get_temperature_from_eint_with_equi
3639 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3641 integer,
intent(in) :: ixi^
l, ixo^
l
3642 double precision,
intent(in) :: w(ixi^s, 1:nw)
3643 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3644 double precision,
intent(out):: res(ixi^s)
3646 double precision :: r(ixi^s)
3648 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3651 end subroutine mhd_get_temperature_equi
3653 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3655 integer,
intent(in) :: ixi^
l, ixo^
l
3656 double precision,
intent(in) :: w(ixi^s, 1:nw)
3657 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3658 double precision,
intent(out):: res(ixi^s)
3660 end subroutine mhd_get_rho_equi
3662 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3664 integer,
intent(in) :: ixi^
l, ixo^
l
3665 double precision,
intent(in) :: w(ixi^s, 1:nw)
3666 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3667 double precision,
intent(out):: res(ixi^s)
3669 end subroutine mhd_get_pe_equi
3672 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3676 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3678 double precision,
intent(in) :: wc(ixi^s,nw)
3680 double precision,
intent(in) :: w(ixi^s,nw)
3681 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3682 double precision,
intent(out) :: f(ixi^s,nwflux)
3684 double precision :: vhall(ixi^s,1:
ndir)
3685 double precision :: ptotal
3689 {
do ix^db=ixomin^db,ixomax^db\}
3702 {
do ix^db=ixomin^db,ixomax^db\}
3706 ^
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_)\
3707 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3709 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3712 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3713 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3715 ^
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_)\
3719 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3720 {
do ix^db=ixomin^db,ixomax^db\}
3721 if(total_energy)
then
3723 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3724 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3727 ^
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))\
3731 {
do ix^db=ixomin^db,ixomax^db\}
3732 f(ix^d,mag(idim))=w(ix^d,
psi_)
3734 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3739 {
do ix^db=ixomin^db,ixomax^db\}
3745 {
do ix^db=ixomin^db,ixomax^db\}
3746 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)
3751 end subroutine mhd_get_flux
3754 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3758 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3760 double precision,
intent(in) :: wc(ixi^s,nw)
3762 double precision,
intent(in) :: w(ixi^s,nw)
3763 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3764 double precision,
intent(out) :: f(ixi^s,nwflux)
3766 double precision :: vhall(ixi^s,1:
ndir)
3769 {
do ix^db=ixomin^db,ixomax^db\}
3780 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3781 {
do ix^db=ixomin^db,ixomax^db\}
3783 ^
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))\
3787 {
do ix^db=ixomin^db,ixomax^db\}
3788 f(ix^d,mag(idim))=w(ix^d,
psi_)
3790 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3795 {
do ix^db=ixomin^db,ixomax^db\}
3800 end subroutine mhd_get_flux_noe
3803 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3807 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3809 double precision,
intent(in) :: wc(ixi^s,nw)
3811 double precision,
intent(in) :: w(ixi^s,nw)
3812 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3813 double precision,
intent(out) :: f(ixi^s,nwflux)
3815 double precision :: vhall(ixi^s,1:
ndir)
3818 {
do ix^db=ixomin^db,ixomax^db\}
3831 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3832 {
do ix^db=ixomin^db,ixomax^db\}
3834 ^
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))\
3838 {
do ix^db=ixomin^db,ixomax^db\}
3839 f(ix^d,mag(idim))=w(ix^d,
psi_)
3841 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3846 {
do ix^db=ixomin^db,ixomax^db\}
3852 {
do ix^db=ixomin^db,ixomax^db\}
3853 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)
3858 end subroutine mhd_get_flux_hde
3861 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
3865 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3867 double precision,
intent(in) :: wc(ixi^s,nw)
3869 double precision,
intent(in) :: w(ixi^s,nw)
3870 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3871 double precision,
intent(out) :: f(ixi^s,nwflux)
3873 double precision :: vhall(ixi^s,1:
ndir)
3874 double precision :: ptotal, btotal(ixo^s,1:
ndir)
3877 {
do ix^db=ixomin^db,ixomax^db\}
3885 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3889 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
3893 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
3894 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3896 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
3900 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3903 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
3910 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
3911 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
3916 {
do ix^db=ixomin^db,ixomax^db\}
3917 f(ix^d,mag(idim))=w(ix^d,
psi_)
3919 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3924 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3925 {
do ix^db=ixomin^db,ixomax^db\}
3927 ^
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))\
3928 if(total_energy)
then
3930 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
3931 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3937 {
do ix^db=ixomin^db,ixomax^db\}
3942 {
do ix^db=ixomin^db,ixomax^db\}
3943 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
3948 end subroutine mhd_get_flux_split
3951 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
3955 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3957 double precision,
intent(in) :: wc(ixi^s,nw)
3959 double precision,
intent(in) :: w(ixi^s,nw)
3960 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3961 double precision,
intent(out) :: f(ixi^s,nwflux)
3963 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
3966 {
do ix^db=ixomin^db,ixomax^db\}
3971 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
3972 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
3973 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
3978 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
3983 e2=(^
c&e(ix^
d,^
c)**2+)
3990 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
3991 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
3992 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
3995 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
3996 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4009 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4011 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)
4018 {
do ix^db=ixomin^db,ixomax^db\}
4019 f(ix^d,mag(idim))=w(ix^d,
psi_)
4021 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4026 {
do ix^db=ixomin^db,ixomax^db\}
4031 end subroutine mhd_get_flux_semirelati
4033 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4037 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4039 double precision,
intent(in) :: wc(ixi^s,nw)
4041 double precision,
intent(in) :: w(ixi^s,nw)
4042 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4043 double precision,
intent(out) :: f(ixi^s,nwflux)
4045 double precision :: e(ixo^s,1:
ndir),e2
4048 {
do ix^db=ixomin^db,ixomax^db\}
4053 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4054 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4055 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4056 e2=(^
c&e(ix^
d,^
c)**2+)
4061 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4070 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4079 {
do ix^db=ixomin^db,ixomax^db\}
4080 f(ix^d,mag(idim))=w(ix^d,
psi_)
4082 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4087 {
do ix^db=ixomin^db,ixomax^db\}
4092 end subroutine mhd_get_flux_semirelati_noe
4098 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4100 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4101 double precision,
intent(in) :: qdt
4102 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4103 double precision,
intent(inout) :: w(ixi^s,1:nw)
4104 double precision :: tmp(ixi^s)
4105 double precision :: jxbxb(ixi^s,1:3)
4107 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4110 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4112 end subroutine add_source_ambipolar_internal_energy
4114 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4117 integer,
intent(in) :: ixi^
l, ixo^
l
4118 double precision,
intent(in) :: w(ixi^s,nw)
4119 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4120 double precision,
intent(out) :: res(:^
d&,:)
4122 double precision :: btot(ixi^s,1:3)
4123 double precision :: current(ixi^s,7-2*
ndir:3)
4124 double precision :: tmp(ixi^s),b2(ixi^s)
4125 integer :: idir, idirmin
4134 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4137 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4140 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4141 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4143 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4146 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4148 end subroutine mhd_get_jxbxb
4155 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4159 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4160 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4161 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4162 double precision,
intent(in) :: my_dt
4163 logical,
intent(in) :: fix_conserve_at_step
4165 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4166 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4167 double precision :: fe(ixi^s,
sdim:3)
4168 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4169 integer :: i, ixa^
l, ie_
4175 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4185 btot(ixa^s,1:3)=0.d0
4191 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4194 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4195 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4197 wres(ixo^s,
e_)=-tmp2(ixo^s)
4203 ff(ixa^s,1) = tmp(ixa^s,2)
4204 ff(ixa^s,2) = -tmp(ixa^s,1)
4206 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4207 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4208 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4211 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4213 ixamin^
d=ixomin^
d-1;
4214 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4223 ff(ixa^s,2) = tmp(ixa^s,3)
4224 ff(ixa^s,3) = -tmp(ixa^s,2)
4225 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4226 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4228 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4230 ff(ixa^s,1) = -tmp(ixa^s,3)
4232 ff(ixa^s,3) = tmp(ixa^s,1)
4233 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4234 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4235 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4239 ff(ixa^s,1) = tmp(ixa^s,2)
4240 ff(ixa^s,2) = -tmp(ixa^s,1)
4242 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4243 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4244 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4249 if(fix_conserve_at_step)
then
4250 fluxall=my_dt*fluxall
4257 end subroutine sts_set_source_ambipolar
4260 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4263 integer,
intent(in) :: ixi^
l, ixo^
l
4264 double precision,
intent(in) :: w(ixi^s,1:nw)
4265 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4267 double precision,
intent(in) :: ecc(ixi^s,1:3)
4268 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4269 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4271 integer :: hxc^
l,ixc^
l,ixa^
l
4272 integer :: idim1,idim2,idir,ix^
d
4278 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4280 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4281 ixamin^
d=ixcmin^
d+ix^
d;
4282 ixamax^
d=ixcmax^
d+ix^
d;
4283 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4285 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4291 ixcmin^d=ixomin^d-1;
4299 hxc^l=ixc^l-kr(idim2,^d);
4301 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4302 +lvc(idim1,idim2,idir)&
4307 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4310 end subroutine update_faces_ambipolar
4314 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4317 integer,
intent(in) :: ixi^
l, ixo^
l
4318 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4319 double precision,
intent(out) :: src(ixi^s)
4321 double precision :: ffc(ixi^s,1:
ndim)
4322 double precision :: dxinv(
ndim)
4323 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4329 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4331 ixbmin^
d=ixcmin^
d+ix^
d;
4332 ixbmax^
d=ixcmax^
d+ix^
d;
4335 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4337 ff(ixi^s,1:ndim)=0.d0
4339 ixb^l=ixo^l-kr(idims,^d);
4340 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4342 if({ ix^d==0 .and. ^d==idims | .or.})
then
4343 ixbmin^d=ixcmin^d-ix^d;
4344 ixbmax^d=ixcmax^d-ix^d;
4345 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4348 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4351 if(slab_uniform)
then
4353 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4354 ixb^l=ixo^l-kr(idims,^d);
4355 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4359 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4360 ixb^l=ixo^l-kr(idims,^d);
4361 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4363 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4365 end subroutine get_flux_on_cell_face
4369 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4372 integer,
intent(in) :: ixi^
l, ixo^
l
4373 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4374 double precision,
intent(in) :: w(ixi^s,1:nw)
4375 double precision :: dtnew
4377 double precision :: coef
4378 double precision :: dxarr(
ndim)
4379 double precision :: tmp(ixi^s)
4384 coef = maxval(abs(tmp(ixo^s)))
4391 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4393 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4396 end function get_ambipolar_dt
4404 integer,
intent(in) :: ixi^
l, ixo^
l
4405 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4406 double precision,
intent(inout) :: res(ixi^s)
4407 double precision :: tmp(ixi^s)
4408 double precision :: rho(ixi^s)
4417 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4421 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4428 integer,
intent(in) :: ixi^
l, ixo^
l
4429 double precision,
intent(in) :: qdt,dtfactor
4430 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4431 double precision,
intent(inout) :: w(ixi^s,1:nw)
4432 logical,
intent(in) :: qsourcesplit
4433 logical,
intent(inout) :: active
4440 if (.not. qsourcesplit)
then
4444 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4448 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4453 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4459 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4463 if (abs(
mhd_eta)>smalldouble)
then
4465 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4470 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4476 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4480 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4487 select case (type_divb)
4492 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4495 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4498 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4499 case (divb_janhunen)
4501 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4502 case (divb_lindejanhunen)
4504 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4505 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4506 case (divb_lindepowel)
4508 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4509 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4510 case (divb_lindeglm)
4512 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4513 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4514 case (divb_multigrid)
4519 call mpistop(
'Unknown divB fix')
4526 w,x,qsourcesplit,active,
rc_fl)
4536 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4545 if(.not.qsourcesplit)
then
4547 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4551 end subroutine mhd_add_source
4553 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4557 integer,
intent(in) :: ixi^
l, ixo^
l
4558 double precision,
intent(in) :: qdt,dtfactor
4559 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4560 double precision,
intent(inout) :: w(ixi^s,1:nw)
4561 double precision :: divv(ixi^s)
4577 end subroutine add_pe0_divv
4579 subroutine get_tau(ixI^L,ixO^L,w,Te,tau,sigT5)
4581 integer,
intent(in) :: ixi^
l, ixo^
l
4582 double precision,
dimension(ixI^S,1:nw),
intent(in) :: w
4583 double precision,
dimension(ixI^S),
intent(in) :: te
4584 double precision,
dimension(ixI^S),
intent(out) :: tau,sigt5
4586 double precision :: dxmin,taumin
4587 double precision,
dimension(ixI^S) :: sigt7,eint
4595 sigt7(ixo^s)=sigt5(ixo^s)*
block%wextra(ixo^s,
tcoff_)
4598 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4602 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4605 tau(ixo^s)=max(taumin*
dt,sigt7(ixo^s)/eint(ixo^s)/
cmax_global**2)
4606 end subroutine get_tau
4608 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4610 integer,
intent(in) :: ixi^
l,ixo^
l
4611 double precision,
intent(in) :: qdt
4612 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4613 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4614 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4616 double precision :: invdx
4617 double precision,
dimension(ixI^S) :: te,tau,sigt,htc_qsrc,tface,r
4618 double precision,
dimension(ixI^S) :: htc_esrc,bsum,bunit
4619 double precision,
dimension(ixI^S,1:ndim) :: btot
4621 integer :: hxc^
l,hxo^
l,ixc^
l,jxc^
l,jxo^
l,kxc^
l
4623 call mhd_get_rfactor(wctprim,x,ixi^
l,ixi^
l,r)
4625 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*w(ixi^s,
rho_))
4626 call get_tau(ixi^
l,ixo^
l,wctprim,te,tau,sigt)
4630 btot(ixo^s,idims)=wct(ixo^s,mag(idims))+
block%B0(ixo^s,idims,0)
4632 btot(ixo^s,idims)=wct(ixo^s,mag(idims))
4635 bsum(ixo^s)=sqrt(sum(btot(ixo^s,:)**2,dim=
ndim+1))+smalldouble
4639 ixcmin^
d=ixomin^
d-
kr(idims,^
d);ixcmax^
d=ixomax^
d;
4640 jxc^
l=ixc^
l+
kr(idims,^
d);
4641 kxc^
l=jxc^
l+
kr(idims,^
d);
4642 hxc^
l=ixc^
l-
kr(idims,^
d);
4643 hxo^
l=ixo^
l-
kr(idims,^
d);
4644 tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
4645 bunit(ixo^s)=btot(ixo^s,idims)/bsum(ixo^s)
4646 htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*bunit(ixo^s)*(tface(ixo^s)-tface(hxo^s))*invdx
4648 htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,
q_))/tau(ixo^s)
4649 w(ixo^s,
q_)=w(ixo^s,
q_)-qdt*htc_qsrc(ixo^s)
4650 end subroutine add_hypertc_source
4653 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4655 integer,
intent(in) :: ixi^
l, ixo^
l
4656 double precision,
intent(in) :: w(ixi^s,1:nw)
4657 double precision,
intent(inout) :: jxb(ixi^s,3)
4658 double precision :: a(ixi^s,3),
b(ixi^s,3)
4660 double precision :: current(ixi^s,7-2*
ndir:3)
4661 integer :: idir, idirmin
4666 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4670 b(ixo^s, idir) = w(ixo^s,mag(idir))
4679 a(ixo^s,idir)=current(ixo^s,idir)
4683 end subroutine get_lorentz_force
4687 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4689 integer,
intent(in) :: ixi^
l, ixo^
l
4690 double precision,
intent(in) :: w(ixi^s, nw)
4691 double precision,
intent(out) :: gamma_a2(ixo^s)
4692 double precision :: rho(ixi^s)
4698 rho(ixo^s) = w(ixo^s,
rho_)
4701 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4702 end subroutine mhd_gamma2_alfven
4706 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4708 integer,
intent(in) :: ixi^
l, ixo^
l
4709 double precision,
intent(in) :: w(ixi^s, nw)
4710 double precision :: gamma_a(ixo^s)
4712 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4713 gamma_a = sqrt(gamma_a)
4714 end function mhd_gamma_alfven
4718 integer,
intent(in) :: ixi^
l, ixo^
l
4719 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4720 double precision,
intent(out) :: rho(ixi^s)
4725 rho(ixo^s) = w(ixo^s,
rho_)
4731 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4734 integer,
intent(in) :: ixi^
l,ixo^
l, ie
4735 double precision,
intent(inout) :: w(ixi^s,1:nw)
4736 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4737 character(len=*),
intent(in) :: subname
4739 double precision :: rho(ixi^s)
4741 logical :: flag(ixi^s,1:nw)
4745 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4746 flag(ixo^s,ie)=.true.
4748 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4750 if(any(flag(ixo^s,ie)))
then
4754 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4757 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4763 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4766 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4772 end subroutine mhd_handle_small_ei
4774 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
4778 integer,
intent(in) :: ixi^
l, ixo^
l
4779 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4780 double precision,
intent(inout) :: w(ixi^s,1:nw)
4782 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
4791 end subroutine mhd_update_temperature
4794 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4797 integer,
intent(in) :: ixi^
l, ixo^
l
4798 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4799 double precision,
intent(inout) :: w(ixi^s,1:nw)
4801 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4813 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4818 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4821 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4827 if(total_energy)
then
4830 b(ixo^s,:)=wct(ixo^s,mag(:))
4839 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4842 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4846 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4850 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
4855 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4860 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
4862 end subroutine add_source_b0split
4865 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4869 integer,
intent(in) :: ixi^
l, ixo^
l
4870 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4871 double precision,
intent(inout) :: w(ixi^s,1:nw)
4872 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4874 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
4875 integer :: idir, idirmin, ix^
d
4879 {
do ix^db=iximin^db,iximax^db\}
4881 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
4882 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
4883 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
4885 call divvector(e,ixi^l,ixo^l,dive)
4887 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
4890 {
do ix^db=ixomin^db,ixomax^db\}
4891 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
4892 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
4893 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
4894 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
4895 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
4896 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
4900 end subroutine add_source_semirelativistic
4903 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4907 integer,
intent(in) :: ixi^
l, ixo^
l
4908 double precision,
intent(in) :: qdt
4909 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4910 double precision,
intent(inout) :: w(ixi^s,1:nw)
4911 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
4913 double precision :: divv(ixi^s), tmp
4925 {
do ix^db=ixomin^db,ixomax^db\}
4927 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
4928 if(w(ix^
d,
e_)<small_e)
then
4933 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
4936 if(fix_small_values)
then
4937 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
4939 end subroutine add_source_internal_e
4942 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4947 integer,
intent(in) :: ixi^
l, ixo^
l
4948 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4949 double precision,
intent(inout) :: w(ixi^s,1:nw)
4950 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4952 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
4953 double precision :: current(ixi^s,7-2*
ndir:3)
4954 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
4955 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
4956 integer :: idir, idirmin, idims, ix^
d
4961 b(ixo^s, idir) = wct(ixo^s,mag(idir))
4969 j(ixo^s,idir)=current(ixo^s,idir)
4987 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
4993 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5000 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)
5004 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5008 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5009 tmp(ixo^s)=sqrt(b2(ixo^s))
5010 where(tmp(ixo^s)>smalldouble)
5011 tmp(ixo^s)=1.d0/tmp(ixo^s)
5017 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5022 {
do ix^db=ixomin^db,ixomax^db\}
5024 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5026 b2(ix^
d)=vaoc/(1.d0+vaoc)
5029 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5032 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5036 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5039 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5040 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5043 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5046 end subroutine add_source_hydrodynamic_e
5052 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5057 integer,
intent(in) :: ixi^
l, ixo^
l
5058 double precision,
intent(in) :: qdt
5059 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5060 double precision,
intent(inout) :: w(ixi^s,1:nw)
5061 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5062 integer :: lxo^
l, kxo^
l
5064 double precision :: tmp(ixi^s),tmp2(ixi^s)
5067 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5068 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5077 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5078 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5085 gradeta(ixo^s,1:
ndim)=zero
5091 gradeta(ixo^s,idim)=tmp(ixo^s)
5098 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5105 tmp2(ixi^s)=bf(ixi^s,idir)
5107 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5108 jxo^
l=ixo^
l+
kr(idim,^
d);
5109 hxo^
l=ixo^
l-
kr(idim,^
d);
5110 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5111 tmp(ixo^s)=tmp(ixo^s)+&
5112 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5117 tmp2(ixi^s)=bf(ixi^s,idir)
5119 jxo^
l=ixo^
l+
kr(idim,^
d);
5120 hxo^
l=ixo^
l-
kr(idim,^
d);
5121 tmp(ixo^s)=tmp(ixo^s)+&
5122 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5127 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5131 do jdir=1,
ndim;
do kdir=idirmin,3
5132 if (
lvc(idir,jdir,kdir)/=0)
then
5133 if (
lvc(idir,jdir,kdir)==1)
then
5134 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5136 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5143 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5144 if(total_energy)
then
5145 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5151 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5154 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5156 end subroutine add_source_res1
5160 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5165 integer,
intent(in) :: ixi^
l, ixo^
l
5166 double precision,
intent(in) :: qdt
5167 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5168 double precision,
intent(inout) :: w(ixi^s,1:nw)
5171 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5172 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5173 integer :: ixa^
l,idir,idirmin,idirmin1
5177 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5178 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5188 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5193 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5202 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5205 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5210 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5212 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5214 if(total_energy)
then
5217 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5218 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5221 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5225 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5226 end subroutine add_source_res2
5230 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5234 integer,
intent(in) :: ixi^
l, ixo^
l
5235 double precision,
intent(in) :: qdt
5236 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5237 double precision,
intent(inout) :: w(ixi^s,1:nw)
5239 double precision :: current(ixi^s,7-2*
ndir:3)
5240 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5241 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5244 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5245 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5248 tmpvec(ixa^s,1:
ndir)=zero
5250 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5254 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5257 tmpvec(ixa^s,1:
ndir)=zero
5258 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5262 tmpvec2(ixa^s,1:
ndir)=zero
5263 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5266 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5269 if(total_energy)
then
5272 tmpvec2(ixa^s,1:
ndir)=zero
5273 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5274 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5275 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5276 end do;
end do;
end do
5278 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5279 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5282 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5284 end subroutine add_source_hyperres
5286 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5293 integer,
intent(in) :: ixi^
l, ixo^
l
5294 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5295 double precision,
intent(inout) :: w(ixi^s,1:nw)
5297 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5318 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5321 if(total_energy)
then
5330 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5339 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5343 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5345 end subroutine add_source_glm
5348 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5351 integer,
intent(in) :: ixi^
l, ixo^
l
5352 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5353 double precision,
intent(inout) :: w(ixi^s,1:nw)
5355 double precision :: divb(ixi^s), ba(1:
ndir)
5356 integer :: idir, ix^
d
5362 {
do ix^db=ixomin^db,ixomax^db\}
5367 if (total_energy)
then
5373 {
do ix^db=ixomin^db,ixomax^db\}
5375 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5377 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5378 if (total_energy)
then
5380 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5385 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5387 end subroutine add_source_powel
5389 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5394 integer,
intent(in) :: ixi^
l, ixo^
l
5395 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5396 double precision,
intent(inout) :: w(ixi^s,1:nw)
5398 double precision :: divb(ixi^s)
5399 integer :: idir, ix^
d
5404 {
do ix^db=ixomin^db,ixomax^db\}
5409 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5411 end subroutine add_source_janhunen
5413 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5418 integer,
intent(in) :: ixi^
l, ixo^
l
5419 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5420 double precision,
intent(inout) :: w(ixi^s,1:nw)
5422 double precision :: divb(ixi^s),graddivb(ixi^s)
5423 integer :: idim, idir, ixp^
l, i^
d, iside
5424 logical,
dimension(-1:1^D&) :: leveljump
5432 if(i^
d==0|.and.) cycle
5433 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5434 leveljump(i^
d)=.true.
5436 leveljump(i^
d)=.false.
5445 i^dd=kr(^dd,^d)*(2*iside-3);
5446 if (leveljump(i^dd))
then
5448 ixpmin^d=ixomin^d-i^d
5450 ixpmax^d=ixomax^d-i^d
5461 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5463 {
do i^db=ixpmin^db,ixpmax^db\}
5465 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5467 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5469 if (typedivbdiff==
'all' .and. total_energy)
then
5471 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5476 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5478 end subroutine add_source_linde
5485 integer,
intent(in) :: ixi^
l, ixo^
l
5486 double precision,
intent(in) :: w(ixi^s,1:nw)
5487 double precision :: divb(ixi^s), dsurface(ixi^s)
5489 double precision :: invb(ixo^s)
5490 integer :: ixa^
l,idims
5492 call get_divb(w,ixi^
l,ixo^
l,divb)
5494 where(invb(ixo^s)/=0.d0)
5495 invb(ixo^s)=1.d0/invb(ixo^s)
5498 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5500 ixamin^
d=ixomin^
d-1;
5501 ixamax^
d=ixomax^
d-1;
5502 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5504 ixa^
l=ixo^
l-
kr(idims,^
d);
5505 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5507 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5508 block%dvolume(ixo^s)/dsurface(ixo^s)
5519 integer,
intent(in) :: ixo^
l, ixi^
l
5520 double precision,
intent(in) :: w(ixi^s,1:nw)
5521 integer,
intent(out) :: idirmin
5524 double precision :: current(ixi^s,7-2*
ndir:3)
5525 integer :: idir, idirmin0
5531 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5532 block%J0(ixo^s,idirmin0:3)
5536 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5544 integer,
intent(in) :: ixi^
l, ixo^
l
5545 double precision,
intent(inout) :: dtnew
5546 double precision,
intent(in) ::
dx^
d
5547 double precision,
intent(in) :: w(ixi^s,1:nw)
5548 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5550 double precision :: dxarr(
ndim)
5551 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5552 integer :: idirmin,idim
5566 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5569 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5595 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5602 end subroutine mhd_get_dt
5605 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5610 integer,
intent(in) :: ixi^
l, ixo^
l
5611 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5612 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5614 double precision :: tmp,tmp1,invr,cot
5616 integer :: mr_,mphi_
5617 integer :: br_,bphi_
5620 br_=mag(1); bphi_=mag(1)-1+
phi_
5625 {
do ix^db=ixomin^db,ixomax^db\}
5628 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5633 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5638 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5639 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5640 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5641 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5642 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5644 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5645 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5646 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5649 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5654 {
do ix^db=ixomin^db,ixomax^db\}
5656 if(local_timestep)
then
5657 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5662 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5668 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5671 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5672 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5676 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5682 cot=1.d0/tan(x(ix^d,2))
5686 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5687 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5689 if(.not.stagger_grid)
then
5690 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5692 tmp=tmp+wprim(ix^d,
psi_)*cot
5694 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5699 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5700 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5701 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5703 if(.not.stagger_grid)
then
5704 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5706 tmp=tmp+wprim(ix^d,
psi_)*cot
5708 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5711 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5712 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5713 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5714 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5715 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5717 if(.not.stagger_grid)
then
5718 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5719 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5720 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5721 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5722 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5729 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5732 end subroutine mhd_add_source_geom
5735 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5740 integer,
intent(in) :: ixi^
l, ixo^
l
5741 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5742 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5744 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
5746 integer :: mr_,mphi_
5747 integer :: br_,bphi_
5750 br_=mag(1); bphi_=mag(1)-1+
phi_
5755 {
do ix^db=ixomin^db,ixomax^db\}
5758 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5769 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
5770 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
5771 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5776 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5782 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
5783 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
5784 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
5785 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5786 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5787 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
5789 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5790 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5791 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5794 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
5795 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
5800 {
do ix^db=ixomin^db,ixomax^db\}
5802 if(local_timestep)
then
5803 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
5809 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5810 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5811 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5815 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5822 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
5828 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
5831 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
5832 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
5833 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
5837 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
5843 cot=1.d0/tan(x(ix^d,2))
5847 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
5848 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
5850 if(.not.stagger_grid)
then
5851 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5853 tmp=tmp+wprim(ix^d,
psi_)*cot
5855 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5861 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
5862 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
5863 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
5864 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
5866 if(.not.stagger_grid)
then
5867 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5869 tmp=tmp+wprim(ix^d,
psi_)*cot
5871 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5874 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
5875 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
5876 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5877 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
5878 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
5879 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
5880 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
5882 if(.not.stagger_grid)
then
5883 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
5884 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5885 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5886 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5887 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5894 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5897 end subroutine mhd_add_source_geom_semirelati
5900 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5905 integer,
intent(in) :: ixi^
l, ixo^
l
5906 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5907 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5909 double precision :: tmp,tmp1,tmp2,invr,cot
5911 integer :: mr_,mphi_
5912 integer :: br_,bphi_
5915 br_=mag(1); bphi_=mag(1)-1+
phi_
5920 {
do ix^db=ixomin^db,ixomax^db\}
5923 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5928 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5933 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5934 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5935 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5936 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5937 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5939 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5940 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5941 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5944 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5949 {
do ix^db=ixomin^db,ixomax^db\}
5951 if(local_timestep)
then
5952 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5956 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5957 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
5960 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5964 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5965 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
5966 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
5968 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5969 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5974 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5980 cot=1.d0/tan(x(ix^d,2))
5985 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5986 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
5987 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
5989 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5990 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5993 if(.not.stagger_grid)
then
5995 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
5996 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
5998 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6001 tmp=tmp+wprim(ix^d,
psi_)*cot
6003 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6009 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6010 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6011 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6012 +(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)
6014 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6015 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6016 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6019 if(.not.stagger_grid)
then
6021 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6022 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6024 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6027 tmp=tmp+wprim(ix^d,
psi_)*cot
6029 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6033 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6034 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6035 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6036 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6037 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6038 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6039 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6040 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6041 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6043 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6044 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6045 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6046 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6047 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6050 if(.not.stagger_grid)
then
6052 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6053 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6054 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6055 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6056 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6057 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6058 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6059 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6060 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6062 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6063 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6064 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6065 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6066 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6074 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6077 end subroutine mhd_add_source_geom_split
6082 integer,
intent(in) :: ixi^
l, ixo^
l
6083 double precision,
intent(in) :: w(ixi^s, nw)
6084 double precision :: mge(ixo^s)
6087 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6089 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6093 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6096 integer,
intent(in) :: ixi^
l, ixo^
l
6097 double precision,
intent(in) :: w(ixi^s,nw)
6098 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6099 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6101 double precision :: current(ixi^s,7-2*
ndir:3)
6102 double precision :: rho(ixi^s)
6103 integer :: idir, idirmin, ix^
d
6108 do idir = idirmin,
ndir
6109 {
do ix^db=ixomin^db,ixomax^db\}
6110 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6114 end subroutine mhd_getv_hall
6116 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6119 integer,
intent(in) :: ixi^
l, ixo^
l
6120 double precision,
intent(in) :: w(ixi^s,nw)
6121 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6122 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6125 double precision :: current(ixi^s,7-2*
ndir:3)
6126 integer :: idir, idirmin
6133 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6134 do idir = idirmin, 3
6138 end subroutine mhd_get_jambi
6140 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6143 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6144 double precision,
intent(in) :: qt
6145 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6146 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6149 double precision :: db(ixo^s), dpsi(ixo^s)
6153 {
do ix^db=ixomin^db,ixomax^db\}
6154 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6155 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6156 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6157 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6166 {
do ix^db=ixomin^db,ixomax^db\}
6167 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6168 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6169 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6170 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6171 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6173 if(total_energy)
then
6174 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6175 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6177 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6179 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6182 if(total_energy)
then
6183 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6184 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6189 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6191 end subroutine mhd_modify_wlr
6193 subroutine mhd_boundary_adjust(igrid,psb)
6195 integer,
intent(in) :: igrid
6198 integer :: ib, idims, iside, ixo^
l, i^
d
6207 i^
d=
kr(^
d,idims)*(2*iside-3);
6208 if (neighbor_type(i^
d,igrid)/=1) cycle
6209 ib=(idims-1)*2+iside
6227 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6232 end subroutine mhd_boundary_adjust
6234 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6237 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6238 double precision,
intent(inout) :: w(ixg^s,1:nw)
6239 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6241 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6242 integer :: ix^
d,ixf^
l
6255 do ix1=ixfmax1,ixfmin1,-1
6256 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6257 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6258 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6261 do ix1=ixfmax1,ixfmin1,-1
6262 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6263 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6264 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6265 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6266 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6267 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6268 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6282 do ix1=ixfmax1,ixfmin1,-1
6283 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6284 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6285 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6286 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6287 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6288 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6291 do ix1=ixfmax1,ixfmin1,-1
6292 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6293 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6294 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6295 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6296 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6297 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6298 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6299 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6300 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6301 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6302 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6303 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6304 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6305 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6306 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6307 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6308 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6309 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6323 do ix1=ixfmin1,ixfmax1
6324 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6325 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6326 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6329 do ix1=ixfmin1,ixfmax1
6330 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6331 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6332 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6333 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6334 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6335 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6336 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6350 do ix1=ixfmin1,ixfmax1
6351 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6352 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6353 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6354 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6355 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6356 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6359 do ix1=ixfmin1,ixfmax1
6360 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6361 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6362 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6363 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6364 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6365 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6366 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6367 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6368 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6369 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6370 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6371 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6372 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6373 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6374 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6375 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6376 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6377 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6391 do ix2=ixfmax2,ixfmin2,-1
6392 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6393 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6394 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6397 do ix2=ixfmax2,ixfmin2,-1
6398 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6399 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6400 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6401 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6402 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6403 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6404 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6418 do ix2=ixfmax2,ixfmin2,-1
6419 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6420 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6421 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6422 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6423 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6424 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6427 do ix2=ixfmax2,ixfmin2,-1
6428 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6429 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6430 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6431 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6432 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6433 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6434 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6435 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6436 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6437 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6438 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6439 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6440 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6441 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6442 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6443 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6444 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6445 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6459 do ix2=ixfmin2,ixfmax2
6460 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6461 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6462 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6465 do ix2=ixfmin2,ixfmax2
6466 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6467 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6468 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6469 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6470 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6471 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6472 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6486 do ix2=ixfmin2,ixfmax2
6487 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6488 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6489 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6490 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6491 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6492 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6495 do ix2=ixfmin2,ixfmax2
6496 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6497 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6498 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6499 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6500 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6501 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6502 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6503 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6504 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6505 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6506 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6507 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6508 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6509 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6510 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6511 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6512 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6513 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6530 do ix3=ixfmax3,ixfmin3,-1
6531 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6532 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6533 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6534 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6535 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6536 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6539 do ix3=ixfmax3,ixfmin3,-1
6540 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6541 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6542 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6543 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6544 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6545 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6546 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6547 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6548 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6549 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6550 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6551 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6552 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6553 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6554 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6555 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6556 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6557 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6572 do ix3=ixfmin3,ixfmax3
6573 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6574 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6575 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6576 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6577 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6578 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6581 do ix3=ixfmin3,ixfmax3
6582 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6583 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6584 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6585 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6586 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6587 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6588 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6589 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6590 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6591 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6592 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6593 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6594 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6595 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6596 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6597 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6598 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6599 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6605 call mpistop(
"Special boundary is not defined for this region")
6608 end subroutine fixdivb_boundary
6617 double precision,
intent(in) :: qdt
6618 double precision,
intent(in) :: qt
6619 logical,
intent(inout) :: active
6622 integer,
parameter :: max_its = 50
6623 double precision :: residual_it(max_its), max_divb
6624 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6625 double precision :: res
6626 double precision,
parameter :: max_residual = 1
d-3
6627 double precision,
parameter :: residual_reduction = 1
d-10
6628 integer :: iigrid, igrid
6629 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6632 mg%operator_type = mg_laplacian
6640 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6641 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6644 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6645 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6647 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6648 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6651 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6652 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6656 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6657 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6658 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6666 do iigrid = 1, igridstail
6667 igrid = igrids(iigrid);
6670 lvl =
mg%boxes(id)%lvl
6671 nc =
mg%box_size_lvl(lvl)
6677 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6679 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6680 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6685 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6688 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6689 if (
mype == 0) print *,
"iteration vs residual"
6692 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6693 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6694 if (residual_it(n) < residual_reduction * max_divb)
exit
6696 if (
mype == 0 .and. n > max_its)
then
6697 print *,
"divb_multigrid warning: not fully converged"
6698 print *,
"current amplitude of divb: ", residual_it(max_its)
6699 print *,
"multigrid smallest grid: ", &
6700 mg%domain_size_lvl(:,
mg%lowest_lvl)
6701 print *,
"note: smallest grid ideally has <= 8 cells"
6702 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6703 print *,
"note: dx/dy/dz should be similar"
6707 call mg_fas_vcycle(
mg, max_res=res)
6708 if (res < max_residual)
exit
6710 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6715 do iigrid = 1, igridstail
6716 igrid = igrids(iigrid);
6725 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6729 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6731 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
6733 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6746 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6747 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6750 if(total_energy)
then
6752 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6755 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6764 subroutine mhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
6767 integer,
intent(in) :: ixi^
l, ixo^
l
6768 double precision,
intent(in) :: qt,qdt
6770 double precision,
intent(in) :: wprim(ixi^s,1:nw)
6771 type(state) :: sct, s
6772 type(ct_velocity) :: vcts
6773 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6774 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6778 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
6780 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
6782 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
6784 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6787 end subroutine mhd_update_faces
6790 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
6794 integer,
intent(in) :: ixi^
l, ixo^
l
6795 double precision,
intent(in) :: qt, qdt
6796 type(state) :: sct, s
6797 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6798 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6800 double precision :: circ(ixi^s,1:
ndim)
6802 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6803 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
6804 integer :: idim1,idim2,idir,iwdim1,iwdim2
6806 associate(bfaces=>s%ws,x=>s%x)
6813 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6820 i1kr^
d=
kr(idim1,^
d);
6823 i2kr^
d=
kr(idim2,^
d);
6826 if (
lvc(idim1,idim2,idir)==1)
then
6828 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6830 {
do ix^db=ixcmin^db,ixcmax^db\}
6831 fe(ix^
d,idir)=quarter*&
6832 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
6833 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
6835 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
6840 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
6848 if(
associated(usr_set_electric_field)) &
6849 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6851 circ(ixi^s,1:ndim)=zero
6856 ixcmin^d=ixomin^d-kr(idim1,^d);
6858 ixa^l=ixc^l-kr(idim2,^d);
6861 if(lvc(idim1,idim2,idir)==1)
then
6863 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6866 else if(lvc(idim1,idim2,idir)==-1)
then
6868 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6874 {
do ix^db=ixcmin^db,ixcmax^db\}
6876 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
6878 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
6885 end subroutine update_faces_average
6888 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
6893 integer,
intent(in) :: ixi^
l, ixo^
l
6894 double precision,
intent(in) :: qt, qdt
6896 double precision,
intent(in) :: wp(ixi^s,1:nw)
6897 type(state) :: sct, s
6898 type(ct_velocity) :: vcts
6899 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6900 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6902 double precision :: circ(ixi^s,1:
ndim)
6904 double precision :: ecc(ixi^s,
sdim:3)
6905 double precision :: ein(ixi^s,
sdim:3)
6907 double precision :: el(ixi^s),er(ixi^s)
6909 double precision :: elc,erc
6911 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6913 double precision :: jce(ixi^s,
sdim:3)
6915 double precision :: xs(ixgs^t,1:
ndim)
6916 double precision :: gradi(ixgs^t)
6917 integer :: ixc^
l,ixa^
l
6918 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
6920 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
6923 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6929 {
do ix^db=iximin^db,iximax^db\}
6932 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_)
6933 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_)
6934 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_)
6937 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
6944 {
do ix^db=iximin^db,iximax^db\}
6947 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
6948 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
6949 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6952 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6966 i1kr^d=kr(idim1,^d);
6969 i2kr^d=kr(idim2,^d);
6972 if (lvc(idim1,idim2,idir)==1)
then
6974 ixcmin^d=ixomin^d+kr(idir,^d)-1;
6977 {
do ix^db=ixcmin^db,ixcmax^db\}
6978 fe(ix^d,idir)=quarter*&
6979 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
6980 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
6985 ixamax^d=ixcmax^d+i1kr^d;
6986 {
do ix^db=ixamin^db,ixamax^db\}
6987 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
6988 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
6991 do ix^db=ixcmin^db,ixcmax^db\}
6992 if(vnorm(ix^d,idim1)>0.d0)
then
6994 else if(vnorm(ix^d,idim1)<0.d0)
then
6995 elc=el({ix^d+i1kr^d})
6997 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
6999 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7001 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7002 erc=er({ix^d+i1kr^d})
7004 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7006 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7011 ixamax^d=ixcmax^d+i2kr^d;
7012 {
do ix^db=ixamin^db,ixamax^db\}
7013 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7014 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7017 do ix^db=ixcmin^db,ixcmax^db\}
7018 if(vnorm(ix^d,idim2)>0.d0)
then
7020 else if(vnorm(ix^d,idim2)<0.d0)
then
7021 elc=el({ix^d+i2kr^d})
7023 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7025 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7027 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7028 erc=er({ix^d+i2kr^d})
7030 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7032 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7036 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7041 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7055 if (lvc(idim1,idim2,idir)==0) cycle
7057 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7058 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7061 xs(ixa^s,:)=x(ixa^s,:)
7062 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7063 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7064 if (lvc(idim1,idim2,idir)==1)
then
7065 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7067 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7074 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7076 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7080 {
do ix^db=ixomin^db,ixomax^db\}
7081 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7082 +ein(ix1,ix2-1,ix3-1,idir))
7083 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7084 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7086 else if(idir==2)
then
7087 {
do ix^db=ixomin^db,ixomax^db\}
7088 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7089 +ein(ix1-1,ix2,ix3-1,idir))
7090 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7091 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7094 {
do ix^db=ixomin^db,ixomax^db\}
7095 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7096 +ein(ix1-1,ix2-1,ix3,idir))
7097 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7098 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7104 {
do ix^db=ixomin^db,ixomax^db\}
7105 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7106 +ein(ix1-1,ix2-1,idir))
7107 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7108 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7113 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7119 if(
associated(usr_set_electric_field)) &
7120 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7122 circ(ixi^s,1:ndim)=zero
7127 ixcmin^d=ixomin^d-kr(idim1,^d);
7129 ixa^l=ixc^l-kr(idim2,^d);
7132 if(lvc(idim1,idim2,idir)==1)
then
7134 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7137 else if(lvc(idim1,idim2,idir)==-1)
then
7139 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7145 {
do ix^db=ixcmin^db,ixcmax^db\}
7147 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7149 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7156 end subroutine update_faces_contact
7159 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
7164 integer,
intent(in) :: ixi^
l, ixo^
l
7165 double precision,
intent(in) :: qt, qdt
7166 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7167 type(state) :: sct, s
7168 type(ct_velocity) :: vcts
7170 double precision :: vtill(ixi^s,2)
7171 double precision :: vtilr(ixi^s,2)
7172 double precision :: bfacetot(ixi^s,
ndim)
7173 double precision :: btill(ixi^s,
ndim)
7174 double precision :: btilr(ixi^s,
ndim)
7175 double precision :: cp(ixi^s,2)
7176 double precision :: cm(ixi^s,2)
7177 double precision :: circ(ixi^s,1:
ndim)
7179 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7180 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7181 integer :: idim1,idim2,idir,ix^
d
7183 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7184 cbarmax=>vcts%cbarmax)
7197 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
7213 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7217 idim2=mod(idir+1,3)+1
7219 jxc^
l=ixc^
l+
kr(idim1,^
d);
7220 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7224 vtill(ixi^s,2),vtilr(ixi^s,2))
7227 vtill(ixi^s,1),vtilr(ixi^s,1))
7233 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7234 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7236 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7237 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7240 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7243 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7247 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7248 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7250 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7251 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7255 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7256 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7257 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7258 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7259 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7260 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7261 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7262 /(cp(ixc^s,2)+cm(ixc^s,2))
7265 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7269 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7283 circ(ixi^s,1:
ndim)=zero
7288 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7292 if(
lvc(idim1,idim2,idir)/=0)
then
7293 hxc^
l=ixc^
l-
kr(idim2,^
d);
7295 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7296 +
lvc(idim1,idim2,idir)&
7302 {
do ix^db=ixcmin^db,ixcmax^db\}
7304 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7306 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7312 end subroutine update_faces_hll
7315 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
7320 integer,
intent(in) :: ixi^
l, ixo^
l
7321 type(state),
intent(in) :: sct, s
7323 double precision :: jce(ixi^s,
sdim:3)
7326 double precision :: jcc(ixi^s,7-2*
ndir:3)
7328 double precision :: xs(ixgs^t,1:
ndim)
7330 double precision :: eta(ixi^s)
7331 double precision :: gradi(ixgs^t)
7332 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7334 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7340 if (
lvc(idim1,idim2,idir)==0) cycle
7342 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7343 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7346 xs(ixb^s,:)=x(ixb^s,:)
7347 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7348 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7349 if (
lvc(idim1,idim2,idir)==1)
then
7350 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7352 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7359 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7367 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7368 jcc(ixc^s,idir)=0.d0
7370 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7371 ixamin^
d=ixcmin^
d+ix^
d;
7372 ixamax^
d=ixcmax^
d+ix^
d;
7373 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7375 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7376 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7381 end subroutine get_resistive_electric_field
7384 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7387 integer,
intent(in) :: ixi^
l, ixo^
l
7388 double precision,
intent(in) :: w(ixi^s,1:nw)
7389 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7390 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7392 double precision :: jxbxb(ixi^s,1:3)
7393 integer :: idir,ixa^
l,ixc^
l,ix^
d
7396 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7403 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7406 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7407 ixamin^
d=ixcmin^
d+ix^
d;
7408 ixamax^
d=ixcmax^
d+ix^
d;
7409 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7411 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7414 end subroutine get_ambipolar_electric_field
7420 integer,
intent(in) :: ixo^
l
7430 do ix^db=ixomin^db,ixomax^db\}
7432 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7433 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7434 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7435 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7436 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7437 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7440 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7441 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7442 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7443 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7486 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7487 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7488 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7490 double precision :: adummy(ixis^s,1:3)
7496 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7499 integer,
intent(in) :: ixi^
l, ixo^
l
7500 double precision,
intent(in) :: w(ixi^s,1:nw)
7501 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7502 double precision,
intent(out):: rfactor(ixi^s)
7504 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7508 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)
7510 end subroutine rfactor_from_temperature_ionization
7512 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7514 integer,
intent(in) :: ixi^
l, ixo^
l
7515 double precision,
intent(in) :: w(ixi^s,1:nw)
7516 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7517 double precision,
intent(out):: rfactor(ixi^s)
7521 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
double precision dt
global time step
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
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, rhov, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
subroutine ionization_degree_init()
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
integer, public, protected c_
logical, public, protected mhd_gravity
Whether gravity is added.
logical, public has_equi_rho0
whether split off equilibrium density
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
logical, public, protected mhd_glm_extended
Whether extended GLM-MHD is used with additional sources.
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
logical, public, protected mhd_hyperbolic_thermal_conduction
Whether thermal conduction is used.
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
subroutine, public mhd_e_to_ei(ixil, ixol, w, x)
Transform total energy to internal energy.
double precision, public mhd_adiab
The adiabatic constant.
logical, public, protected mhd_partial_ionization
Whether plasma is partially ionized.
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
double precision, public, protected rr
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
double precision, public mhd_gamma
The adiabatic index.
integer, public, protected mhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
logical, public, protected mhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected mhd_semirelativistic
Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
integer, public, protected mhd_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public, protected q_
Index of the heat flux q.
integer, public, protected mhd_n_tracer
Number of tracer species.
integer, public, protected te_
Indices of temperature.
integer, public, protected m
integer, public equi_rho0_
equi vars indices in the stateequi_vars array
integer, public, protected mhd_trac_type
Which TRAC method is used.
logical, public, protected mhd_cak_force
Whether CAK radiation line force is activated.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
logical, public, protected mhd_hall
Whether Hall-MHD is used.
type(te_fluid), allocatable, public te_fl_mhd
type of fluid for thermal emission synthesis
logical, public, protected mhd_ambipolar
Whether Ambipolar term is used.
double precision, public hypertc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixil, ixol)
Compute 2 times total magnetic energy.
subroutine, public multiplyambicoef(ixil, ixol, res, w, x)
multiply res by the ambipolar coefficient The ambipolar coefficient is calculated as -mhd_eta_ambi/rh...
logical, public partial_energy
Whether an internal or hydrodynamic energy equation is used.
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.
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_glm
Whether GLM-MHD is used to control div B.
logical, public clean_initial_divb
clean initial divB
procedure(sub_convert), pointer, public mhd_to_conserved
double precision, public mhd_eta
The MHD resistivity.
logical, public divbwave
Add divB wave in Roe solver.
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
double precision, public, protected mhd_trac_mask
Height of the mask used in the TRAC method.
procedure(mask_subroutine), pointer, public usr_mask_ambipolar
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
procedure(sub_get_pthermal), pointer, public mhd_get_temperature
integer, public equi_pe0_
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
procedure(sub_convert), pointer, public mhd_to_primitive
logical, public has_equi_pe0
whether split off equilibrium thermal pressure
logical, public, protected mhd_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
logical, public, protected mhd_particles
Whether particles module is added.
integer, public, protected b
subroutine, public mhd_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
subroutine, public get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
double precision, public mhd_etah
Hall resistivity.
subroutine, public mhd_get_v(w, x, ixil, ixol, v)
Calculate v vector.
double precision, public mhd_eta_ambi
The MHD ambipolar coefficient.
logical, public, protected mhd_hydrodynamic_e
Whether hydrodynamic energy is solved instead of total energy.
subroutine, public mhd_phys_init()
logical, public, protected mhd_trac
Whether TRAC method is used.
logical, public, protected eq_state_units
type(rc_fluid), allocatable, public rc_fl
type of fluid for radiative cooling
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
integer, public, protected rho_
Index of the density (in the w array)
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
integer, public, protected tweight_
logical, public, protected mhd_ambipolar_sts
Whether Ambipolar term is implemented using supertimestepping.
procedure(sub_get_pthermal), pointer, public mhd_get_pthermal
subroutine, public mhd_ei_to_e(ixil, ixol, w, x)
Transform internal energy to total energy.
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected mhd_4th_order
MHD fourth order.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine, public mhd_get_rho(w, x, ixil, ixol, rho)
integer, public, protected psi_
Indices of the GLM psi.
logical, public mhd_equi_thermal
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_add_source(qdt, dtfactor, ixil, ixol, wct, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rotating_frame_init()
Initialize the module.
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicut timestep for the TC (mhd implementation)
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine, public tc_get_mhd_params(fl, read_mhd_params)
Init TC coefficients: MHD case.
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
procedure(special_resistivity), pointer usr_special_resistivity
procedure(phys_gravity), pointer usr_gravity
procedure(set_equi_vars), pointer usr_set_equi_vars
procedure(set_electric_field), pointer usr_set_electric_field
The module add viscous source terms and check time step.
subroutine viscosity_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
subroutine viscosity_get_dt(w, ixil, ixol, dtnew, dxd, x)
The data structure that contains information about a tree node/grid block.