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
568 allocate(iw_vector(nvector))
569 iw_vector(1) =
mom(1) - 1
570 iw_vector(2) = mag(1) - 1
573 if (.not.
allocated(flux_type))
then
574 allocate(flux_type(
ndir, nwflux))
575 flux_type = flux_default
576 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
577 call mpistop(
"phys_check error: flux_type has wrong shape")
580 if(nwflux>mag(
ndir))
then
582 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
587 flux_type(:,
psi_)=flux_special
589 flux_type(idir,mag(idir))=flux_special
593 flux_type(idir,mag(idir))=flux_tvdlf
599 phys_get_dt => mhd_get_dt
602 phys_get_cmax => mhd_get_cmax_semirelati
604 phys_get_cmax => mhd_get_cmax_semirelati_noe
608 phys_get_cmax => mhd_get_cmax_origin
610 phys_get_cmax => mhd_get_cmax_origin_noe
613 phys_get_a2max => mhd_get_a2max
614 phys_get_tcutoff => mhd_get_tcutoff
615 phys_get_h_speed => mhd_get_h_speed
617 phys_get_cbounds => mhd_get_cbounds_split_rho
619 phys_get_cbounds => mhd_get_cbounds_semirelati
621 phys_get_cbounds => mhd_get_cbounds
624 phys_to_primitive => mhd_to_primitive_hde
626 phys_to_conserved => mhd_to_conserved_hde
630 phys_to_primitive => mhd_to_primitive_semirelati
632 phys_to_conserved => mhd_to_conserved_semirelati
635 phys_to_primitive => mhd_to_primitive_semirelati_noe
637 phys_to_conserved => mhd_to_conserved_semirelati_noe
642 phys_to_primitive => mhd_to_primitive_split_rho
644 phys_to_conserved => mhd_to_conserved_split_rho
647 phys_to_primitive => mhd_to_primitive_inte
649 phys_to_conserved => mhd_to_conserved_inte
652 phys_to_primitive => mhd_to_primitive_origin
654 phys_to_conserved => mhd_to_conserved_origin
657 phys_to_primitive => mhd_to_primitive_origin_noe
659 phys_to_conserved => mhd_to_conserved_origin_noe
664 phys_get_flux => mhd_get_flux_hde
667 phys_get_flux => mhd_get_flux_semirelati
669 phys_get_flux => mhd_get_flux_semirelati_noe
673 phys_get_flux => mhd_get_flux_split
675 phys_get_flux => mhd_get_flux
677 phys_get_flux => mhd_get_flux_noe
682 phys_add_source_geom => mhd_add_source_geom_semirelati
684 phys_add_source_geom => mhd_add_source_geom_split
686 phys_add_source_geom => mhd_add_source_geom
688 phys_add_source => mhd_add_source
689 phys_check_params => mhd_check_params
690 phys_write_info => mhd_write_info
693 phys_handle_small_values => mhd_handle_small_values_inte
694 mhd_handle_small_values => mhd_handle_small_values_inte
695 phys_check_w => mhd_check_w_inte
697 phys_handle_small_values => mhd_handle_small_values_hde
698 mhd_handle_small_values => mhd_handle_small_values_hde
699 phys_check_w => mhd_check_w_hde
701 phys_handle_small_values => mhd_handle_small_values_semirelati
702 mhd_handle_small_values => mhd_handle_small_values_semirelati
703 phys_check_w => mhd_check_w_semirelati
705 phys_handle_small_values => mhd_handle_small_values_origin
706 mhd_handle_small_values => mhd_handle_small_values_origin
707 phys_check_w => mhd_check_w_origin
709 phys_handle_small_values => mhd_handle_small_values_noe
710 mhd_handle_small_values => mhd_handle_small_values_noe
711 phys_check_w => mhd_check_w_noe
715 phys_get_pthermal => mhd_get_pthermal_inte
718 phys_get_pthermal => mhd_get_pthermal_hde
721 phys_get_pthermal => mhd_get_pthermal_semirelati
724 phys_get_pthermal => mhd_get_pthermal_origin
727 phys_get_pthermal => mhd_get_pthermal_noe
732 phys_set_equi_vars => set_equi_vars_grid
735 if(type_divb==divb_glm)
then
736 phys_modify_wlr => mhd_modify_wlr
741 mhd_get_rfactor=>rfactor_from_temperature_ionization
742 phys_update_temperature => mhd_update_temperature
746 mhd_get_rfactor=>rfactor_from_constant_ionization
769 phys_get_ct_velocity => mhd_get_ct_velocity
770 phys_update_faces => mhd_update_faces
772 phys_modify_wlr => mhd_modify_wlr
774 phys_boundary_adjust => mhd_boundary_adjust
783 call mhd_physical_units()
789 call mpistop(
"thermal conduction needs mhd_energy=T")
792 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
795 call mpistop(
"radiative cooling needs mhd_energy=T")
806 if(phys_internal_e)
then
808 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
810 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
814 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot_with_equi
816 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
820 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
822 tc_fl%has_equi = .true.
823 tc_fl%get_temperature_equi => mhd_get_temperature_equi
824 tc_fl%get_rho_equi => mhd_get_rho_equi
826 tc_fl%has_equi = .false.
829 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
853 rc_fl%get_var_Rfactor => mhd_get_rfactor
857 rc_fl%has_equi = .true.
858 rc_fl%get_rho_equi => mhd_get_rho_equi
859 rc_fl%get_pthermal_equi => mhd_get_pe_equi
861 rc_fl%has_equi = .false.
867 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
869 phys_te_images => mhd_te_images
885 if (particles_eta < zero) particles_eta =
mhd_eta
886 if (particles_etah < zero) particles_eta =
mhd_etah
888 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
889 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
903 phys_wider_stencil = 2
905 phys_wider_stencil = 1
913 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
925 phys_wider_stencil = 2
927 phys_wider_stencil = 1
941 subroutine mhd_te_images
946 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
948 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
950 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
952 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
955 call mpistop(
"Error in synthesize emission: Unknown convert_type")
957 end subroutine mhd_te_images
963 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
967 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
968 double precision,
intent(in) :: x(ixi^s,1:
ndim)
969 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
970 double precision,
intent(in) :: my_dt
971 logical,
intent(in) :: fix_conserve_at_step
973 end subroutine mhd_sts_set_source_tc_mhd
975 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
982 integer,
intent(in) :: ixi^
l, ixo^
l
983 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
984 double precision,
intent(in) :: w(ixi^s,1:nw)
985 double precision :: dtnew
988 end function mhd_get_tc_dt_mhd
990 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
993 integer,
intent(in) :: ixi^
l,ixo^
l
994 double precision,
intent(inout) :: w(ixi^s,1:nw)
995 double precision,
intent(in) :: x(ixi^s,1:
ndim)
996 integer,
intent(in) :: step
997 character(len=140) :: error_msg
999 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1000 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1001 end subroutine mhd_tc_handle_small_e
1004 subroutine tc_params_read_mhd(fl)
1006 type(tc_fluid),
intent(inout) :: fl
1008 double precision :: tc_k_para=0d0
1009 double precision :: tc_k_perp=0d0
1012 logical :: tc_perpendicular=.false.
1013 logical :: tc_saturate=.false.
1014 character(len=std_len) :: tc_slope_limiter=
"MC"
1016 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1020 read(
unitpar, tc_list,
end=111)
1024 fl%tc_perpendicular = tc_perpendicular
1025 fl%tc_saturate = tc_saturate
1026 fl%tc_k_para = tc_k_para
1027 fl%tc_k_perp = tc_k_perp
1028 select case(tc_slope_limiter)
1030 fl%tc_slope_limiter = 0
1033 fl%tc_slope_limiter = 1
1036 fl%tc_slope_limiter = 2
1039 fl%tc_slope_limiter = 3
1042 fl%tc_slope_limiter = 4
1044 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1046 end subroutine tc_params_read_mhd
1050 subroutine rc_params_read(fl)
1053 type(rc_fluid),
intent(inout) :: fl
1055 double precision :: cfrac=0.1d0
1058 double precision :: rad_cut_hgt=0.5d0
1059 double precision :: rad_cut_dey=0.15d0
1062 integer :: ncool = 4000
1064 logical :: tfix=.false.
1066 logical :: rc_split=.false.
1067 logical :: rad_cut=.false.
1069 character(len=std_len) :: coolcurve=
'JCcorona'
1071 character(len=std_len) :: coolmethod=
'exact'
1073 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1077 read(
unitpar, rc_list,
end=111)
1082 fl%coolcurve=coolcurve
1083 fl%coolmethod=coolmethod
1086 fl%rc_split=rc_split
1089 fl%rad_cut_hgt=rad_cut_hgt
1090 fl%rad_cut_dey=rad_cut_dey
1091 end subroutine rc_params_read
1095 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1098 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1099 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1101 double precision :: delx(ixi^s,1:
ndim)
1102 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1103 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1109 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1113 hxo^
l=ixo^
l-
kr(idims,^
d);
1119 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1122 xshift^
d=half*(one-
kr(^
d,idims));
1129 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1133 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1136 end subroutine set_equi_vars_grid_faces
1139 subroutine set_equi_vars_grid(igrid)
1143 integer,
intent(in) :: igrid
1149 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1151 end subroutine set_equi_vars_grid
1154 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1156 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1157 double precision,
intent(in) :: w(ixi^s, 1:nw)
1158 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1159 double precision :: wnew(ixo^s, 1:nwc)
1166 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1172 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1176 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1180 if(
b0field .and. total_energy)
then
1181 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1182 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1186 end function convert_vars_splitting
1188 subroutine mhd_check_params
1196 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1197 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1201 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1202 inv_gamma_1=1.d0/gamma_1
1207 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1212 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1217 end subroutine mhd_check_params
1219 subroutine mhd_physical_units()
1221 double precision :: mp,kb,miu0,c_lightspeed
1222 double precision :: a,
b
1233 c_lightspeed=const_c
1294 end subroutine mhd_physical_units
1296 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1299 logical,
intent(in) :: primitive
1300 logical,
intent(inout) :: flag(ixi^s,1:nw)
1301 integer,
intent(in) :: ixi^
l, ixo^
l
1302 double precision,
intent(in) :: w(ixi^s,nw)
1304 double precision :: tmp,b2,
b(ixo^s,1:
ndir)
1305 double precision :: v(ixo^s,1:
ndir),gamma2,inv_rho
1316 {
do ix^db=ixomin^db,ixomax^db \}
1317 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1320 {
do ix^db=ixomin^db,ixomax^db \}
1321 b2=(^
c&w(ix^d,
b^
c_)**2+)
1322 if(b2>smalldouble)
then
1327 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1328 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1329 inv_rho = 1d0/w(ix^d,
rho_)
1331 b2=b2*inv_rho*inv_squared_c
1333 gamma2=1.d0/(1.d0+b2)
1335 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1338 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1339 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1340 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1345 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1351 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1352 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1353 if(tmp<small_e) flag(ix^d,
e_)=.true.
1359 end subroutine mhd_check_w_semirelati
1361 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1364 logical,
intent(in) :: primitive
1365 integer,
intent(in) :: ixi^
l, ixo^
l
1366 double precision,
intent(in) :: w(ixi^s,nw)
1367 logical,
intent(inout) :: flag(ixi^s,1:nw)
1369 double precision :: tmp
1373 {
do ix^db=ixomin^db,ixomax^db\}
1387 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1391 if(tmp<small_e) flag(ix^
d,
e_) = .true.
1396 end subroutine mhd_check_w_origin
1398 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1401 logical,
intent(in) :: primitive
1402 integer,
intent(in) :: ixi^
l, ixo^
l
1403 double precision,
intent(in) :: w(ixi^s,nw)
1404 logical,
intent(inout) :: flag(ixi^s,1:nw)
1409 {
do ix^db=ixomin^db,ixomax^db\}
1417 end subroutine mhd_check_w_noe
1419 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1422 logical,
intent(in) :: primitive
1423 integer,
intent(in) :: ixi^
l, ixo^
l
1424 double precision,
intent(in) :: w(ixi^s,nw)
1425 logical,
intent(inout) :: flag(ixi^s,1:nw)
1430 {
do ix^db=ixomin^db,ixomax^db\}
1446 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1451 end subroutine mhd_check_w_inte
1453 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1456 logical,
intent(in) :: primitive
1457 integer,
intent(in) :: ixi^
l, ixo^
l
1458 double precision,
intent(in) :: w(ixi^s,nw)
1459 logical,
intent(inout) :: flag(ixi^s,1:nw)
1464 {
do ix^db=ixomin^db,ixomax^db\}
1469 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1473 end subroutine mhd_check_w_hde
1476 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1478 integer,
intent(in) :: ixi^
l, ixo^
l
1479 double precision,
intent(inout) :: w(ixi^s, nw)
1480 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1484 {
do ix^db=ixomin^db,ixomax^db\}
1486 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1488 +(^
c&w(ix^
d,
b^
c_)**2+))
1493 end subroutine mhd_to_conserved_origin
1496 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1498 integer,
intent(in) :: ixi^
l, ixo^
l
1499 double precision,
intent(inout) :: w(ixi^s, nw)
1500 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1504 {
do ix^db=ixomin^db,ixomax^db\}
1509 end subroutine mhd_to_conserved_origin_noe
1512 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1514 integer,
intent(in) :: ixi^
l, ixo^
l
1515 double precision,
intent(inout) :: w(ixi^s, nw)
1516 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1520 {
do ix^db=ixomin^db,ixomax^db\}
1522 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1528 end subroutine mhd_to_conserved_hde
1531 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1533 integer,
intent(in) :: ixi^
l, ixo^
l
1534 double precision,
intent(inout) :: w(ixi^s, nw)
1535 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1539 {
do ix^db=ixomin^db,ixomax^db\}
1541 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1546 end subroutine mhd_to_conserved_inte
1549 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1551 integer,
intent(in) :: ixi^
l, ixo^
l
1552 double precision,
intent(inout) :: w(ixi^s, nw)
1553 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1555 double precision :: rho
1558 {
do ix^db=ixomin^db,ixomax^db\}
1561 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1562 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1563 +(^
c&w(ix^
d,
b^
c_)**2+))
1568 end subroutine mhd_to_conserved_split_rho
1571 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1573 integer,
intent(in) :: ixi^
l, ixo^
l
1574 double precision,
intent(inout) :: w(ixi^s, nw)
1575 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1577 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1580 {
do ix^db=ixomin^db,ixomax^db\}
1582 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1583 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1584 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1585 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1586 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1587 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1592 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1593 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1594 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1602 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1606 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1608 +(^
c&w(ix^
d,
b^
c_)**2+)&
1609 +(^
c&e(ix^
d,^
c)**2+)*inv_squared_c)
1617 end subroutine mhd_to_conserved_semirelati
1619 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1621 integer,
intent(in) :: ixi^
l, ixo^
l
1622 double precision,
intent(inout) :: w(ixi^s, nw)
1623 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1625 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1628 {
do ix^db=ixomin^db,ixomax^db\}
1630 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1631 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1632 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1633 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1634 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1635 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1640 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1641 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1642 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1652 end subroutine mhd_to_conserved_semirelati_noe
1655 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1657 integer,
intent(in) :: ixi^
l, ixo^
l
1658 double precision,
intent(inout) :: w(ixi^s, nw)
1659 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1661 double precision :: inv_rho
1666 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1669 {
do ix^db=ixomin^db,ixomax^db\}
1670 inv_rho = 1.d0/w(ix^
d,
rho_)
1674 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1676 +(^
c&w(ix^
d,
b^
c_)**2+)))
1679 end subroutine mhd_to_primitive_origin
1682 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1684 integer,
intent(in) :: ixi^
l, ixo^
l
1685 double precision,
intent(inout) :: w(ixi^s, nw)
1686 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1688 double precision :: inv_rho
1693 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1696 {
do ix^db=ixomin^db,ixomax^db\}
1697 inv_rho = 1.d0/w(ix^
d,
rho_)
1702 end subroutine mhd_to_primitive_origin_noe
1705 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1707 integer,
intent(in) :: ixi^
l, ixo^
l
1708 double precision,
intent(inout) :: w(ixi^s, nw)
1709 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1711 double precision :: inv_rho
1716 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1719 {
do ix^db=ixomin^db,ixomax^db\}
1720 inv_rho = 1d0/w(ix^
d,
rho_)
1727 end subroutine mhd_to_primitive_hde
1730 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1732 integer,
intent(in) :: ixi^
l, ixo^
l
1733 double precision,
intent(inout) :: w(ixi^s, nw)
1734 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1736 double precision :: inv_rho
1741 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1744 {
do ix^db=ixomin^db,ixomax^db\}
1746 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1748 inv_rho = 1.d0/w(ix^
d,
rho_)
1752 end subroutine mhd_to_primitive_inte
1755 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1757 integer,
intent(in) :: ixi^
l, ixo^
l
1758 double precision,
intent(inout) :: w(ixi^s, nw)
1759 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1761 double precision :: inv_rho
1766 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1769 {
do ix^db=ixomin^db,ixomax^db\}
1774 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1776 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1779 end subroutine mhd_to_primitive_split_rho
1782 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1784 integer,
intent(in) :: ixi^
l, ixo^
l
1785 double precision,
intent(inout) :: w(ixi^s, nw)
1786 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1788 double precision ::
b(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
1793 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1796 {
do ix^db=ixomin^db,ixomax^db\}
1797 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1798 if(b2>smalldouble)
then
1806 inv_rho=1.d0/w(ix^
d,
rho_)
1808 b2=b2*inv_rho*inv_squared_c
1810 gamma2=1.d0/(1.d0+b2)
1812 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1816 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1820 b(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1821 b(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1822 b(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1826 b(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1832 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1834 +(^
c&w(ix^
d,
b^
c_)**2+)&
1835 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
1839 end subroutine mhd_to_primitive_semirelati
1842 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1844 integer,
intent(in) :: ixi^
l, ixo^
l
1845 double precision,
intent(inout) :: w(ixi^s, nw)
1846 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1848 double precision ::
b(ixo^s,1:
ndir),tmp,b2,gamma2,inv_rho
1849 integer :: ix^
d, idir
1853 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1856 {
do ix^db=ixomin^db,ixomax^db\}
1857 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1858 if(b2>smalldouble)
then
1866 inv_rho=1.d0/w(ix^
d,
rho_)
1868 b2=b2*inv_rho*inv_squared_c
1870 gamma2=1.d0/(1.d0+b2)
1872 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1875 end subroutine mhd_to_primitive_semirelati_noe
1880 integer,
intent(in) :: ixi^
l, ixo^
l
1881 double precision,
intent(inout) :: w(ixi^s, nw)
1882 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1887 {
do ix^db=ixomin^db,ixomax^db\}
1890 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1892 +(^
c&w(ix^
d,
b^
c_)**2+))
1895 {
do ix^db=ixomin^db,ixomax^db\}
1897 w(ix^d,
e_)=w(ix^d,
e_)&
1898 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1899 +(^
c&w(ix^d,
b^
c_)**2+))
1906 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
1908 integer,
intent(in) :: ixi^
l, ixo^
l
1909 double precision,
intent(inout) :: w(ixi^s, nw)
1910 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1915 {
do ix^db=ixomin^db,ixomax^db\}
1918 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1922 {
do ix^db=ixomin^db,ixomax^db\}
1924 w(ix^d,
e_)=w(ix^d,
e_)&
1925 +half*(^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)
1929 end subroutine mhd_ei_to_e_hde
1932 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
1934 integer,
intent(in) :: ixi^
l, ixo^
l
1935 double precision,
intent(inout) :: w(ixi^s, nw)
1936 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1938 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
1939 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
1941 end subroutine mhd_ei_to_e_semirelati
1946 integer,
intent(in) :: ixi^
l, ixo^
l
1947 double precision,
intent(inout) :: w(ixi^s, nw)
1948 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1953 {
do ix^db=ixomin^db,ixomax^db\}
1956 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1958 +(^
c&w(ix^
d,
b^
c_)**2+))
1961 {
do ix^db=ixomin^db,ixomax^db\}
1963 w(ix^d,
e_)=w(ix^d,
e_)&
1964 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1965 +(^
c&w(ix^d,
b^
c_)**2+))
1969 if(fix_small_values)
then
1970 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
1976 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
1978 integer,
intent(in) :: ixi^
l, ixo^
l
1979 double precision,
intent(inout) :: w(ixi^s, nw)
1980 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1985 {
do ix^db=ixomin^db,ixomax^db\}
1988 -half*(^
c&w(ix^
d,
m^
c_)**2+)/&
1992 {
do ix^db=ixomin^db,ixomax^db\}
1994 w(ix^d,
e_)=w(ix^d,
e_)&
1995 -half*(^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)
1999 if(fix_small_values)
then
2000 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2003 end subroutine mhd_e_to_ei_hde
2006 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2008 integer,
intent(in) :: ixi^
l, ixo^
l
2009 double precision,
intent(inout) :: w(ixi^s, nw)
2010 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2012 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2013 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2015 end subroutine mhd_e_to_ei_semirelati
2017 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2020 logical,
intent(in) :: primitive
2021 integer,
intent(in) :: ixi^
l,ixo^
l
2022 double precision,
intent(inout) :: w(ixi^s,1:nw)
2023 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2024 character(len=*),
intent(in) :: subname
2026 double precision ::
b(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2027 double precision :: tmp, b2, gamma2, inv_rho
2029 logical :: flag(ixi^s,1:nw)
2038 {
do ix^db=iximin^db,iximax^db\}
2039 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2040 if(b2>smalldouble)
then
2047 inv_rho=1.d0/w(ix^
d,
rho_)
2049 b2=b2*inv_rho*inv_squared_c
2051 gamma2=1.d0/(1.d0+b2)
2053 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
2056 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2057 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2058 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2062 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2068 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2069 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2070 +(^
c&w(ix^
d,
b^
c_)**2+)&
2071 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
2078 select case (small_values_method)
2080 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = small_density
2082 if(small_values_fix_iw(
m^
c_))
then
2083 if(flag({ix^d},
rho_)) w({ix^d},
m^
c_)=0.0d0
2088 where(flag(ixo^s,
e_)) w(ixo^s,
p_) = small_pressure
2090 {
do ix^db=ixomin^db,ixomax^db\}
2091 if(flag(ix^d,
e_))
then
2092 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2093 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2100 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2103 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2105 w(ixi^s,
e_)=pressure(ixi^s)
2106 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2107 {
do ix^db=iximin^db,iximax^db\}
2108 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2109 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2114 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2118 end subroutine mhd_handle_small_values_semirelati
2120 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2123 logical,
intent(in) :: primitive
2124 integer,
intent(in) :: ixi^
l,ixo^
l
2125 double precision,
intent(inout) :: w(ixi^s,1:nw)
2126 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2127 character(len=*),
intent(in) :: subname
2129 double precision :: rho
2130 integer :: idir, ix^
d
2131 logical :: flag(ixi^s,1:nw)
2133 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2138 {
do ix^db=ixomin^db,ixomax^db\}
2148 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2160 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2164 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))
2170 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2172 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2175 {
do ix^db=iximin^db,iximax^db\}
2181 w(ix^d,
e_)=w(ix^d,
e_)&
2182 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2184 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2186 {
do ix^db=iximin^db,iximax^db\}
2192 w(ix^d,
e_)=w(ix^d,
e_)&
2193 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2197 if(.not.primitive)
then
2199 {
do ix^db=iximin^db,iximax^db\}
2205 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2206 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2207 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2210 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2214 end subroutine mhd_handle_small_values_origin
2216 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2219 logical,
intent(in) :: primitive
2220 integer,
intent(in) :: ixi^
l,ixo^
l
2221 double precision,
intent(inout) :: w(ixi^s,1:nw)
2222 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2223 character(len=*),
intent(in) :: subname
2225 double precision :: rho
2227 logical :: flag(ixi^s,1:nw)
2229 call phys_check_w(primitive, ixi^
l, ixi^
l, w, flag)
2234 {
do ix^db=ixomin^db,ixomax^db\}
2242 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2263 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2265 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2267 if(.not.primitive)
then
2269 {
do ix^db=iximin^db,iximax^db\}
2275 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2276 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2279 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2283 end subroutine mhd_handle_small_values_inte
2285 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2288 logical,
intent(in) :: primitive
2289 integer,
intent(in) :: ixi^
l,ixo^
l
2290 double precision,
intent(inout) :: w(ixi^s,1:nw)
2291 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2292 character(len=*),
intent(in) :: subname
2295 logical :: flag(ixi^s,1:nw)
2297 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2303 where(flag(ixo^s,
rho_)) w(ixo^s,
rho_) = &
2310 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2317 if(.not.primitive)
then
2321 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))/(w(ixo^s,
rho_)+&
2326 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
2334 end subroutine mhd_handle_small_values_noe
2336 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2339 logical,
intent(in) :: primitive
2340 integer,
intent(in) :: ixi^
l,ixo^
l
2341 double precision,
intent(inout) :: w(ixi^s,1:nw)
2342 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2343 character(len=*),
intent(in) :: subname
2346 logical :: flag(ixi^s,1:nw)
2348 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2356 where(flag(ixo^s,
rho_)) w(ixo^s,
mom(idir)) = 0.0d0
2362 where(flag(ixo^s,
e_))
2363 w(ixo^s,
e_) = small_e+half*sum(w(ixo^s,
mom(:))**2,dim=
ndim+1)/w(ixo^s,
rho_)
2372 if(.not.primitive)
then
2375 w(ixo^s,
p_)=gamma_1*(w(ixo^s,
e_)-half*sum(w(ixo^s,
mom(:))**2,dim=
ndim+1)/w(ixo^s,
rho_))
2378 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))/w(ixo^s,
rho_)
2385 end subroutine mhd_handle_small_values_hde
2391 integer,
intent(in) :: ixi^
l, ixo^
l
2392 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2393 double precision,
intent(out) :: v(ixi^s,
ndir)
2395 double precision :: rho(ixi^s)
2400 rho(ixo^s)=1.d0/rho(ixo^s)
2403 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2409 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2412 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2413 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2414 double precision,
intent(inout) :: cmax(ixi^s)
2416 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2422 {
do ix^db=ixomin^db,ixomax^db \}
2433 cfast2=b2*inv_rho+cmax(ix^
d)
2434 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2435 if(avmincs2<zero) avmincs2=zero
2436 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2440 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2442 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2445 {
do ix^db=ixomin^db,ixomax^db \}
2455 b2=(^
c&w(ix^d,
b^
c_)**2+)
2456 cfast2=b2*inv_rho+cmax(ix^d)
2457 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2458 if(avmincs2<zero) avmincs2=zero
2459 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2463 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2465 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2469 end subroutine mhd_get_cmax_origin
2472 subroutine mhd_get_cmax_origin_noe(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_noe
2535 subroutine mhd_get_cmax_semirelati(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 :: csound, avmincs2, idim_alfven_speed2
2543 double precision :: inv_rho, alfven_speed2, gamma2
2546 {
do ix^db=ixomin^db,ixomax^db \}
2547 inv_rho=1.d0/w(ix^
d,
rho_)
2548 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2549 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2550 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2553 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2556 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2557 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2558 if(avmincs2<zero) avmincs2=zero
2560 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2561 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2564 end subroutine mhd_get_cmax_semirelati
2567 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2570 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2571 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2572 double precision,
intent(inout):: cmax(ixi^s)
2574 double precision :: csound, avmincs2, idim_alfven_speed2
2575 double precision :: inv_rho, alfven_speed2, gamma2
2578 {
do ix^db=ixomin^db,ixomax^db \}
2579 inv_rho=1.d0/w(ix^
d,
rho_)
2580 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2581 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2582 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2584 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2587 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2588 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2589 if(avmincs2<zero) avmincs2=zero
2591 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2592 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2595 end subroutine mhd_get_cmax_semirelati_noe
2597 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2600 integer,
intent(in) :: ixi^
l, ixo^
l
2601 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2602 double precision,
intent(inout) :: a2max(
ndim)
2603 double precision :: a2(ixi^s,
ndim,nw)
2604 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2609 hxo^
l=ixo^
l-
kr(i,^
d);
2610 gxo^
l=hxo^
l-
kr(i,^
d);
2611 jxo^
l=ixo^
l+
kr(i,^
d);
2612 kxo^
l=jxo^
l+
kr(i,^
d);
2613 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2614 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2615 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2617 end subroutine mhd_get_a2max
2620 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2623 integer,
intent(in) :: ixi^
l,ixo^
l
2624 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2626 double precision,
intent(inout) :: w(ixi^s,1:nw)
2627 double precision,
intent(out) :: tco_local,tmax_local
2629 double precision,
parameter :: trac_delta=0.25d0
2630 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
2631 double precision,
dimension(ixI^S,1:ndir) :: bunitvec
2632 double precision,
dimension(ixI^S,1:ndim) :: gradt
2633 double precision :: bdir(
ndim)
2634 double precision :: ltrc,ltrp,altr(ixi^s)
2635 integer :: idims,jxo^
l,hxo^
l,ixa^
d,ixb^
d,ix^
d
2636 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2637 logical :: lrlt(ixi^s)
2640 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2642 call mhd_get_rfactor(w,x,ixi^
l,ixi^
l,te)
2643 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2646 tmax_local=maxval(te(ixo^s))
2656 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2658 where(lts(ixo^s) > trac_delta)
2661 if(any(lrlt(ixo^s)))
then
2662 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2673 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2674 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2675 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2676 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2678 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2690 gradt(ixo^s,idims)=tmp1(ixo^s)
2694 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+
block%B0(ixo^s,:,0)
2696 bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2702 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2705 if(sum(bdir(:)**2) .gt. zero)
then
2706 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2708 block%special_values(3:ndim+2)=bdir(1:ndim)
2710 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2711 where(tmp1(ixo^s)/=0.d0)
2712 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2714 tmp1(ixo^s)=bigdouble
2718 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2721 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2723 if(slab_uniform)
then
2724 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2726 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2729 where(lts(ixo^s) > trac_delta)
2732 if(any(lrlt(ixo^s)))
then
2733 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2735 block%special_values(1)=zero
2737 block%special_values(2)=tmax_local
2756 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2757 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2758 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2761 {
do ix^db=ixpmin^db,ixpmax^db\}
2763 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))+block%B0(ix^d,^
c,0)\
2765 ^
c&bunitvec(ix^d,^
c)=w(ix^d,iw_mag(^
c))\
2767 tmp1(ix^d)=1.d0/(dsqrt(^
c&bunitvec(ix^d,^
c)**2+)+smalldouble)
2769 ^d&bunitvec({ix^d},^d)=bunitvec({ix^d},^d)*tmp1({ix^d})\
2771 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec({ix^d},^d)+)/te(ix^d)
2773 if(slab_uniform)
then
2774 lts(ix^d)=min(^d&dxlevel(^d))*lts(ix^d)
2776 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2778 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2784 hxo^l=ixp^l-kr(idims,^d);
2785 jxo^l=ixp^l+kr(idims,^d);
2787 altr(ixp^s)=0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2789 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2792 block%wextra(ixp^s,
tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
2796 call mpistop(
"unknown mhd_trac_type")
2799 end subroutine mhd_get_tcutoff
2802 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2805 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2806 double precision,
intent(in) :: wprim(ixi^s, nw)
2807 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2808 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2810 double precision :: csound(ixi^s,
ndim)
2811 double precision,
allocatable :: tmp(:^
d&)
2812 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2816 allocate(tmp(ixa^s))
2818 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2819 csound(ixa^s,id)=tmp(ixa^s)
2822 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2823 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2824 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2825 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))
2829 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2830 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2831 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)))
2832 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2833 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2834 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)))
2839 ixamax^
d=jxcmax^
d+
kr(id,^
d);
2840 ixamin^
d=jxcmin^
d+
kr(id,^
d);
2841 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)))
2842 ixamax^
d=jxcmax^
d-
kr(id,^
d);
2843 ixamin^
d=jxcmin^
d-
kr(id,^
d);
2844 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)))
2848 end subroutine mhd_get_h_speed
2851 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2854 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2855 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2856 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2857 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2858 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
2859 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
2860 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2862 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
2863 double precision :: umean, dmean, tmp1, tmp2, tmp3
2870 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
2871 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
2872 if(
present(cmin))
then
2873 {
do ix^db=ixomin^db,ixomax^db\}
2874 tmp1=sqrt(wlp(ix^
d,
rho_))
2875 tmp2=sqrt(wrp(ix^
d,
rho_))
2876 tmp3=1.d0/(tmp1+tmp2)
2877 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
2878 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
2879 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
2880 cmin(ix^
d,1)=umean-dmean
2881 cmax(ix^
d,1)=umean+dmean
2883 if(h_correction)
then
2884 {
do ix^db=ixomin^db,ixomax^db\}
2885 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2886 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2890 {
do ix^db=ixomin^db,ixomax^db\}
2891 tmp1=sqrt(wlp(ix^d,
rho_))
2892 tmp2=sqrt(wrp(ix^d,
rho_))
2893 tmp3=1.d0/(tmp1+tmp2)
2894 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
2895 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
2896 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
2897 cmax(ix^d,1)=abs(umean)+dmean
2901 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
2902 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2903 if(
present(cmin))
then
2904 {
do ix^db=ixomin^db,ixomax^db\}
2905 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
2906 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
2908 if(h_correction)
then
2909 {
do ix^db=ixomin^db,ixomax^db\}
2910 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2911 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2915 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
2919 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2920 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2921 if(
present(cmin))
then
2922 {
do ix^db=ixomin^db,ixomax^db\}
2923 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2924 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
2925 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2927 if(h_correction)
then
2928 {
do ix^db=ixomin^db,ixomax^db\}
2929 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2930 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2934 {
do ix^db=ixomin^db,ixomax^db\}
2935 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2936 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2941 end subroutine mhd_get_cbounds
2944 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2947 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2948 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2949 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2950 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2951 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
2952 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
2953 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2955 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
2960 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
2961 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
2963 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
2964 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
2966 if(
present(cmin))
then
2967 {
do ix^db=ixomin^db,ixomax^db\}
2968 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
2969 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
2970 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
2973 {
do ix^db=ixomin^db,ixomax^db\}
2974 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
2975 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
2979 end subroutine mhd_get_cbounds_semirelati
2982 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2985 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2986 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
2987 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
2988 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2989 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
2990 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
2991 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
2993 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
2994 double precision :: umean, dmean, tmp1, tmp2, tmp3
3001 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3002 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3003 if(
present(cmin))
then
3004 {
do ix^db=ixomin^db,ixomax^db\}
3007 tmp3=1.d0/(tmp1+tmp2)
3008 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3009 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3010 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3011 cmin(ix^
d,1)=umean-dmean
3012 cmax(ix^
d,1)=umean+dmean
3014 if(h_correction)
then
3015 {
do ix^db=ixomin^db,ixomax^db\}
3016 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3017 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3021 {
do ix^db=ixomin^db,ixomax^db\}
3024 tmp3=1.d0/(tmp1+tmp2)
3025 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3026 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3027 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3028 cmax(ix^d,1)=abs(umean)+dmean
3032 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3033 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3034 if(
present(cmin))
then
3035 {
do ix^db=ixomin^db,ixomax^db\}
3036 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3037 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3039 if(h_correction)
then
3040 {
do ix^db=ixomin^db,ixomax^db\}
3041 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3042 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3046 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3050 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3051 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3052 if(
present(cmin))
then
3053 {
do ix^db=ixomin^db,ixomax^db\}
3054 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3055 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3056 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3058 if(h_correction)
then
3059 {
do ix^db=ixomin^db,ixomax^db\}
3060 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3061 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3065 {
do ix^db=ixomin^db,ixomax^db\}
3066 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3067 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3072 end subroutine mhd_get_cbounds_split_rho
3075 subroutine mhd_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3078 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3079 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3080 double precision,
intent(in) :: cmax(ixi^s)
3081 double precision,
intent(in),
optional :: cmin(ixi^s)
3082 type(ct_velocity),
intent(inout):: vcts
3084 integer :: idime,idimn
3090 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3092 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3094 if(.not.
allocated(vcts%vbarC))
then
3095 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3096 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3099 if(
present(cmin))
then
3100 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3101 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3103 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3104 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3107 idimn=mod(idim,
ndir)+1
3108 idime=mod(idim+1,
ndir)+1
3110 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3111 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3112 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3113 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3114 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3116 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3117 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3118 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3119 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3120 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3122 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
3125 end subroutine mhd_get_ct_velocity
3128 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3131 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3132 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3133 double precision,
intent(out):: csound(ixo^s)
3135 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3142 {
do ix^db=ixomin^db,ixomax^db \}
3143 inv_rho=1.d0/w(ix^
d,
rho_)
3150 cfast2=b2*inv_rho+csound(ix^
d)
3151 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3153 if(avmincs2<zero) avmincs2=zero
3154 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3156 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3160 {
do ix^db=ixomin^db,ixomax^db \}
3161 inv_rho=1.d0/w(ix^d,
rho_)
3167 b2=(^
c&w(ix^d,
b^
c_)**2+)
3168 cfast2=b2*inv_rho+csound(ix^d)
3169 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3170 if(avmincs2<zero) avmincs2=zero
3171 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3173 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3178 end subroutine mhd_get_csound_prim
3181 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3184 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3185 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3186 double precision,
intent(out):: csound(ixo^s)
3188 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3195 {
do ix^db=ixomin^db,ixomax^db \}
3202 cfast2=b2*inv_rho+csound(ix^
d)
3203 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3205 if(avmincs2<zero) avmincs2=zero
3206 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3208 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3212 {
do ix^db=ixomin^db,ixomax^db \}
3218 b2=(^
c&w(ix^d,
b^
c_)**2+)
3219 cfast2=b2*inv_rho+csound(ix^d)
3220 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3221 if(avmincs2<zero) avmincs2=zero
3222 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3224 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3229 end subroutine mhd_get_csound_prim_split
3232 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3235 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3237 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3238 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3240 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3243 {
do ix^db=ixomin^db,ixomax^db\}
3244 inv_rho = 1.d0/w(ix^
d,
rho_)
3247 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3248 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3249 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3250 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3253 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3254 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3255 if(avmincs2<zero) avmincs2=zero
3257 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3260 end subroutine mhd_get_csound_semirelati
3263 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3266 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3268 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3269 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3271 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3274 {
do ix^db=ixomin^db,ixomax^db\}
3275 inv_rho = 1.d0/w(ix^
d,
rho_)
3278 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3279 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3280 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3281 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3284 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3285 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3286 if(avmincs2<zero) avmincs2=zero
3288 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3291 end subroutine mhd_get_csound_semirelati_noe
3294 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3297 integer,
intent(in) :: ixi^
l, ixo^
l
3298 double precision,
intent(in) :: w(ixi^s,nw)
3299 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3300 double precision,
intent(out):: pth(ixi^s)
3308 end subroutine mhd_get_pthermal_noe
3311 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3315 integer,
intent(in) :: ixi^
l, ixo^
l
3316 double precision,
intent(in) :: w(ixi^s,nw)
3317 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3318 double precision,
intent(out):: pth(ixi^s)
3322 {
do ix^db= ixomin^db,ixomax^db\}
3326 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3331 if(check_small_values.and..not.fix_small_values)
then
3332 {
do ix^db= ixomin^db,ixomax^db\}
3333 if(pth(ix^d)<small_pressure)
then
3334 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3335 " encountered when call mhd_get_pthermal_inte"
3336 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3337 write(*,*)
"Location: ", x(ix^d,:)
3338 write(*,*)
"Cell number: ", ix^d
3340 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3343 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3344 write(*,*)
"Saving status at the previous time step"
3350 end subroutine mhd_get_pthermal_inte
3353 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3357 integer,
intent(in) :: ixi^
l, ixo^
l
3358 double precision,
intent(in) :: w(ixi^s,nw)
3359 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3360 double precision,
intent(out):: pth(ixi^s)
3364 {
do ix^db=ixomin^db,ixomax^db\}
3369 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3370 +(^
c&w(ix^
d,
b^
c_)**2+)))
3375 if(check_small_values.and..not.fix_small_values)
then
3376 {
do ix^db=ixomin^db,ixomax^db\}
3377 if(pth(ix^d)<small_pressure)
then
3378 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3379 " encountered when call mhd_get_pthermal"
3380 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3381 write(*,*)
"Location: ", x(ix^d,:)
3382 write(*,*)
"Cell number: ", ix^d
3384 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3387 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3388 write(*,*)
"Saving status at the previous time step"
3394 end subroutine mhd_get_pthermal_origin
3397 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3401 integer,
intent(in) :: ixi^
l, ixo^
l
3402 double precision,
intent(in) :: w(ixi^s,nw)
3403 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3404 double precision,
intent(out):: pth(ixi^s)
3406 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3409 {
do ix^db=ixomin^db,ixomax^db\}
3410 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3411 if(b2>smalldouble)
then
3419 inv_rho=1.d0/w(ix^
d,
rho_)
3421 b2=b2*inv_rho*inv_squared_c
3423 gamma2=1.d0/(1.d0+b2)
3425 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3429 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3430 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3431 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3435 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3441 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3442 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3443 +(^
c&w(ix^
d,
b^
c_)**2+)&
3444 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3448 if(check_small_values.and..not.fix_small_values)
then
3449 {
do ix^db=ixomin^db,ixomax^db\}
3450 if(pth(ix^d)<small_pressure)
then
3451 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3452 " encountered when call mhd_get_pthermal_semirelati"
3453 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3454 write(*,*)
"Location: ", x(ix^d,:)
3455 write(*,*)
"Cell number: ", ix^d
3457 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3460 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3461 write(*,*)
"Saving status at the previous time step"
3467 end subroutine mhd_get_pthermal_semirelati
3470 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3474 integer,
intent(in) :: ixi^
l, ixo^
l
3475 double precision,
intent(in) :: w(ixi^s,nw)
3476 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3477 double precision,
intent(out):: pth(ixi^s)
3481 {
do ix^db= ixomin^db,ixomax^db\}
3482 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3485 if(check_small_values.and..not.fix_small_values)
then
3486 {
do ix^db= ixomin^db,ixomax^db\}
3487 if(pth(ix^d)<small_pressure)
then
3488 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3489 " encountered when call mhd_get_pthermal_hde"
3490 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3491 write(*,*)
"Location: ", x(ix^d,:)
3492 write(*,*)
"Cell number: ", ix^d
3494 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3497 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3498 write(*,*)
"Saving status at the previous time step"
3504 end subroutine mhd_get_pthermal_hde
3507 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3509 integer,
intent(in) :: ixi^
l, ixo^
l
3510 double precision,
intent(in) :: w(ixi^s, 1:nw)
3511 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3512 double precision,
intent(out):: res(ixi^s)
3513 res(ixo^s) = w(ixo^s,
te_)
3514 end subroutine mhd_get_temperature_from_te
3517 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3519 integer,
intent(in) :: ixi^
l, ixo^
l
3520 double precision,
intent(in) :: w(ixi^s, 1:nw)
3521 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3522 double precision,
intent(out):: res(ixi^s)
3524 double precision :: r(ixi^s)
3526 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3527 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3528 end subroutine mhd_get_temperature_from_eint
3531 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3533 integer,
intent(in) :: ixi^
l, ixo^
l
3534 double precision,
intent(in) :: w(ixi^s, 1:nw)
3535 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3536 double precision,
intent(out):: res(ixi^s)
3538 double precision :: r(ixi^s)
3540 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3542 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3544 end subroutine mhd_get_temperature_from_etot
3546 subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
3548 integer,
intent(in) :: ixi^
l, ixo^
l
3549 double precision,
intent(in) :: w(ixi^s, 1:nw)
3550 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3551 double precision,
intent(out):: res(ixi^s)
3553 double precision :: r(ixi^s)
3555 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3559 end subroutine mhd_get_temperature_from_etot_with_equi
3561 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3563 integer,
intent(in) :: ixi^
l, ixo^
l
3564 double precision,
intent(in) :: w(ixi^s, 1:nw)
3565 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3566 double precision,
intent(out):: res(ixi^s)
3568 double precision :: r(ixi^s)
3570 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3574 end subroutine mhd_get_temperature_from_eint_with_equi
3576 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3578 integer,
intent(in) :: ixi^
l, ixo^
l
3579 double precision,
intent(in) :: w(ixi^s, 1:nw)
3580 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3581 double precision,
intent(out):: res(ixi^s)
3583 double precision :: r(ixi^s)
3585 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3588 end subroutine mhd_get_temperature_equi
3590 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3592 integer,
intent(in) :: ixi^
l, ixo^
l
3593 double precision,
intent(in) :: w(ixi^s, 1:nw)
3594 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3595 double precision,
intent(out):: res(ixi^s)
3597 end subroutine mhd_get_rho_equi
3599 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3601 integer,
intent(in) :: ixi^
l, ixo^
l
3602 double precision,
intent(in) :: w(ixi^s, 1:nw)
3603 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3604 double precision,
intent(out):: res(ixi^s)
3606 end subroutine mhd_get_pe_equi
3609 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3613 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3615 double precision,
intent(in) :: wc(ixi^s,nw)
3617 double precision,
intent(in) :: w(ixi^s,nw)
3618 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3619 double precision,
intent(out) :: f(ixi^s,nwflux)
3621 double precision :: vhall(ixi^s,1:
ndir)
3622 double precision :: ptotal
3626 {
do ix^db=ixomin^db,ixomax^db\}
3639 {
do ix^db=ixomin^db,ixomax^db\}
3643 ^
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_)\
3644 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3646 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3649 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3650 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3652 ^
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_)\
3656 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3657 {
do ix^db=ixomin^db,ixomax^db\}
3658 if(total_energy)
then
3660 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3661 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3664 ^
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))\
3668 {
do ix^db=ixomin^db,ixomax^db\}
3669 f(ix^d,mag(idim))=w(ix^d,
psi_)
3671 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3676 {
do ix^db=ixomin^db,ixomax^db\}
3682 {
do ix^db=ixomin^db,ixomax^db\}
3683 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)
3688 end subroutine mhd_get_flux
3691 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3695 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3697 double precision,
intent(in) :: wc(ixi^s,nw)
3699 double precision,
intent(in) :: w(ixi^s,nw)
3700 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3701 double precision,
intent(out) :: f(ixi^s,nwflux)
3703 double precision :: vhall(ixi^s,1:
ndir)
3706 {
do ix^db=ixomin^db,ixomax^db\}
3717 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3718 {
do ix^db=ixomin^db,ixomax^db\}
3720 ^
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))\
3724 {
do ix^db=ixomin^db,ixomax^db\}
3725 f(ix^d,mag(idim))=w(ix^d,
psi_)
3727 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3732 {
do ix^db=ixomin^db,ixomax^db\}
3737 end subroutine mhd_get_flux_noe
3740 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3744 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3746 double precision,
intent(in) :: wc(ixi^s,nw)
3748 double precision,
intent(in) :: w(ixi^s,nw)
3749 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3750 double precision,
intent(out) :: f(ixi^s,nwflux)
3752 double precision :: vhall(ixi^s,1:
ndir)
3755 {
do ix^db=ixomin^db,ixomax^db\}
3768 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3769 {
do ix^db=ixomin^db,ixomax^db\}
3771 ^
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))\
3775 {
do ix^db=ixomin^db,ixomax^db\}
3776 f(ix^d,mag(idim))=w(ix^d,
psi_)
3778 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3783 {
do ix^db=ixomin^db,ixomax^db\}
3789 {
do ix^db=ixomin^db,ixomax^db\}
3790 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)
3795 end subroutine mhd_get_flux_hde
3798 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
3802 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3804 double precision,
intent(in) :: wc(ixi^s,nw)
3806 double precision,
intent(in) :: w(ixi^s,nw)
3807 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3808 double precision,
intent(out) :: f(ixi^s,nwflux)
3810 double precision :: vhall(ixi^s,1:
ndir)
3811 double precision :: ptotal, btotal(ixo^s,1:
ndir)
3814 {
do ix^db=ixomin^db,ixomax^db\}
3822 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3826 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
3830 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
3831 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3833 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
3837 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3840 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
3847 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
3848 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
3853 {
do ix^db=ixomin^db,ixomax^db\}
3854 f(ix^d,mag(idim))=w(ix^d,
psi_)
3856 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3861 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3862 {
do ix^db=ixomin^db,ixomax^db\}
3864 ^
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))\
3865 if(total_energy)
then
3867 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
3868 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3874 {
do ix^db=ixomin^db,ixomax^db\}
3879 {
do ix^db=ixomin^db,ixomax^db\}
3880 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
3885 end subroutine mhd_get_flux_split
3888 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
3892 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3894 double precision,
intent(in) :: wc(ixi^s,nw)
3896 double precision,
intent(in) :: w(ixi^s,nw)
3897 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3898 double precision,
intent(out) :: f(ixi^s,nwflux)
3900 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
3903 {
do ix^db=ixomin^db,ixomax^db\}
3908 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
3909 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
3910 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
3915 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
3920 e2=(^
c&e(ix^
d,^
c)**2+)
3927 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
3928 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
3929 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
3932 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
3933 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
3946 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
3948 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)
3955 {
do ix^db=ixomin^db,ixomax^db\}
3956 f(ix^d,mag(idim))=w(ix^d,
psi_)
3958 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
3963 {
do ix^db=ixomin^db,ixomax^db\}
3968 end subroutine mhd_get_flux_semirelati
3970 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3974 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3976 double precision,
intent(in) :: wc(ixi^s,nw)
3978 double precision,
intent(in) :: w(ixi^s,nw)
3979 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3980 double precision,
intent(out) :: f(ixi^s,nwflux)
3982 double precision :: e(ixo^s,1:
ndir),e2
3985 {
do ix^db=ixomin^db,ixomax^db\}
3990 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
3991 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
3992 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
3993 e2=(^
c&e(ix^
d,^
c)**2+)
3998 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4007 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4009 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)
4016 {
do ix^db=ixomin^db,ixomax^db\}
4017 f(ix^d,mag(idim))=w(ix^d,
psi_)
4019 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4024 {
do ix^db=ixomin^db,ixomax^db\}
4029 end subroutine mhd_get_flux_semirelati_noe
4035 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4037 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4038 double precision,
intent(in) :: qdt
4039 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4040 double precision,
intent(inout) :: w(ixi^s,1:nw)
4041 double precision :: tmp(ixi^s)
4042 double precision :: jxbxb(ixi^s,1:3)
4044 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4047 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4049 end subroutine add_source_ambipolar_internal_energy
4051 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4054 integer,
intent(in) :: ixi^
l, ixo^
l
4055 double precision,
intent(in) :: w(ixi^s,nw)
4056 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4057 double precision,
intent(out) :: res(:^
d&,:)
4059 double precision :: btot(ixi^s,1:3)
4060 double precision :: current(ixi^s,7-2*
ndir:3)
4061 double precision :: tmp(ixi^s),b2(ixi^s)
4062 integer :: idir, idirmin
4071 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4074 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4077 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4078 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4080 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4083 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4085 end subroutine mhd_get_jxbxb
4092 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4096 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4097 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4098 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4099 double precision,
intent(in) :: my_dt
4100 logical,
intent(in) :: fix_conserve_at_step
4102 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4103 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4104 double precision :: fe(ixi^s,
sdim:3)
4105 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4106 integer :: i, ixa^
l, ie_
4112 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4122 btot(ixa^s,1:3)=0.d0
4128 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4131 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4132 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4134 wres(ixo^s,
e_)=-tmp2(ixo^s)
4140 ff(ixa^s,1) = tmp(ixa^s,2)
4141 ff(ixa^s,2) = -tmp(ixa^s,1)
4143 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4144 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4145 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4148 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4150 ixamin^
d=ixomin^
d-1;
4151 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4160 ff(ixa^s,2) = tmp(ixa^s,3)
4161 ff(ixa^s,3) = -tmp(ixa^s,2)
4162 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4163 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4165 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4167 ff(ixa^s,1) = -tmp(ixa^s,3)
4169 ff(ixa^s,3) = tmp(ixa^s,1)
4170 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4171 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4172 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4176 ff(ixa^s,1) = tmp(ixa^s,2)
4177 ff(ixa^s,2) = -tmp(ixa^s,1)
4179 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4180 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4181 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4186 if(fix_conserve_at_step)
then
4187 fluxall=my_dt*fluxall
4194 end subroutine sts_set_source_ambipolar
4197 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4200 integer,
intent(in) :: ixi^
l, ixo^
l
4201 double precision,
intent(in) :: w(ixi^s,1:nw)
4202 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4204 double precision,
intent(in) :: ecc(ixi^s,1:3)
4205 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4206 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4208 integer :: hxc^
l,ixc^
l,ixa^
l
4209 integer :: idim1,idim2,idir,ix^
d
4215 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4217 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4218 ixamin^
d=ixcmin^
d+ix^
d;
4219 ixamax^
d=ixcmax^
d+ix^
d;
4220 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4222 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4228 ixcmin^d=ixomin^d-1;
4236 hxc^l=ixc^l-kr(idim2,^d);
4238 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4239 +lvc(idim1,idim2,idir)&
4244 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4247 end subroutine update_faces_ambipolar
4251 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4254 integer,
intent(in) :: ixi^
l, ixo^
l
4255 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4256 double precision,
intent(out) :: src(ixi^s)
4258 double precision :: ffc(ixi^s,1:
ndim)
4259 double precision :: dxinv(
ndim)
4260 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4266 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4268 ixbmin^
d=ixcmin^
d+ix^
d;
4269 ixbmax^
d=ixcmax^
d+ix^
d;
4272 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4274 ff(ixi^s,1:ndim)=0.d0
4276 ixb^l=ixo^l-kr(idims,^d);
4277 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4279 if({ ix^d==0 .and. ^d==idims | .or.})
then
4280 ixbmin^d=ixcmin^d-ix^d;
4281 ixbmax^d=ixcmax^d-ix^d;
4282 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4285 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4288 if(slab_uniform)
then
4290 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4291 ixb^l=ixo^l-kr(idims,^d);
4292 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4296 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4297 ixb^l=ixo^l-kr(idims,^d);
4298 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4300 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4302 end subroutine get_flux_on_cell_face
4306 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4309 integer,
intent(in) :: ixi^
l, ixo^
l
4310 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4311 double precision,
intent(in) :: w(ixi^s,1:nw)
4312 double precision :: dtnew
4314 double precision :: coef
4315 double precision :: dxarr(
ndim)
4316 double precision :: tmp(ixi^s)
4321 coef = maxval(abs(tmp(ixo^s)))
4328 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4330 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4333 end function get_ambipolar_dt
4341 integer,
intent(in) :: ixi^
l, ixo^
l
4342 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4343 double precision,
intent(inout) :: res(ixi^s)
4344 double precision :: tmp(ixi^s)
4345 double precision :: rho(ixi^s)
4354 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4358 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4365 integer,
intent(in) :: ixi^
l, ixo^
l
4366 double precision,
intent(in) :: qdt,dtfactor
4367 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4368 double precision,
intent(inout) :: w(ixi^s,1:nw)
4369 logical,
intent(in) :: qsourcesplit
4370 logical,
intent(inout) :: active
4377 if (.not. qsourcesplit)
then
4381 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4385 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4390 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4396 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4400 if (abs(
mhd_eta)>smalldouble)
then
4402 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4407 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4413 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4417 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4424 select case (type_divb)
4429 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4432 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4435 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4436 case (divb_janhunen)
4438 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4439 case (divb_lindejanhunen)
4441 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4442 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4443 case (divb_lindepowel)
4445 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4446 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4447 case (divb_lindeglm)
4449 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4450 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4451 case (divb_multigrid)
4456 call mpistop(
'Unknown divB fix')
4463 w,x,qsourcesplit,active,
rc_fl)
4473 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4482 if(.not.qsourcesplit)
then
4484 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4488 end subroutine mhd_add_source
4490 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4494 integer,
intent(in) :: ixi^
l, ixo^
l
4495 double precision,
intent(in) :: qdt,dtfactor
4496 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4497 double precision,
intent(inout) :: w(ixi^s,1:nw)
4498 double precision :: divv(ixi^s)
4514 end subroutine add_pe0_divv
4516 subroutine get_tau(ixI^L,ixO^L,w,Te,tau,sigT5)
4518 integer,
intent(in) :: ixi^
l, ixo^
l
4519 double precision,
dimension(ixI^S,1:nw),
intent(in) :: w
4520 double precision,
dimension(ixI^S),
intent(in) :: te
4521 double precision,
dimension(ixI^S),
intent(out) :: tau,sigt5
4523 double precision :: dxmin,taumin
4524 double precision,
dimension(ixI^S) :: sigt7,eint
4532 sigt7(ixo^s)=sigt5(ixo^s)*
block%wextra(ixo^s,
tcoff_)
4535 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4539 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4542 tau(ixo^s)=max(taumin*
dt,sigt7(ixo^s)/eint(ixo^s)/
cmax_global**2)
4543 end subroutine get_tau
4545 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4547 integer,
intent(in) :: ixi^
l,ixo^
l
4548 double precision,
intent(in) :: qdt
4549 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4550 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4551 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4553 double precision :: invdx
4554 double precision,
dimension(ixI^S) :: te,tau,sigt,htc_qsrc,tface,r
4555 double precision,
dimension(ixI^S) :: htc_esrc,bsum,bunit
4556 double precision,
dimension(ixI^S,1:ndim) :: btot
4558 integer :: hxc^
l,hxo^
l,ixc^
l,jxc^
l,jxo^
l,kxc^
l
4560 call mhd_get_rfactor(wctprim,x,ixi^
l,ixi^
l,r)
4562 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*w(ixi^s,
rho_))
4563 call get_tau(ixi^
l,ixo^
l,wctprim,te,tau,sigt)
4567 btot(ixo^s,idims)=wct(ixo^s,mag(idims))+
block%B0(ixo^s,idims,0)
4569 btot(ixo^s,idims)=wct(ixo^s,mag(idims))
4572 bsum(ixo^s)=sqrt(sum(btot(ixo^s,:)**2,dim=
ndim+1))+smalldouble
4576 ixcmin^
d=ixomin^
d-
kr(idims,^
d);ixcmax^
d=ixomax^
d;
4577 jxc^
l=ixc^
l+
kr(idims,^
d);
4578 kxc^
l=jxc^
l+
kr(idims,^
d);
4579 hxc^
l=ixc^
l-
kr(idims,^
d);
4580 hxo^
l=ixo^
l-
kr(idims,^
d);
4581 tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
4582 bunit(ixo^s)=btot(ixo^s,idims)/bsum(ixo^s)
4583 htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*bunit(ixo^s)*(tface(ixo^s)-tface(hxo^s))*invdx
4585 htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,
q_))/tau(ixo^s)
4586 w(ixo^s,
q_)=w(ixo^s,
q_)-qdt*htc_qsrc(ixo^s)
4587 end subroutine add_hypertc_source
4590 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4592 integer,
intent(in) :: ixi^
l, ixo^
l
4593 double precision,
intent(in) :: w(ixi^s,1:nw)
4594 double precision,
intent(inout) :: jxb(ixi^s,3)
4595 double precision :: a(ixi^s,3),
b(ixi^s,3)
4597 double precision :: current(ixi^s,7-2*
ndir:3)
4598 integer :: idir, idirmin
4603 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4607 b(ixo^s, idir) = w(ixo^s,mag(idir))
4616 a(ixo^s,idir)=current(ixo^s,idir)
4620 end subroutine get_lorentz_force
4624 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4626 integer,
intent(in) :: ixi^
l, ixo^
l
4627 double precision,
intent(in) :: w(ixi^s, nw)
4628 double precision,
intent(out) :: gamma_a2(ixo^s)
4629 double precision :: rho(ixi^s)
4635 rho(ixo^s) = w(ixo^s,
rho_)
4638 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4639 end subroutine mhd_gamma2_alfven
4643 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4645 integer,
intent(in) :: ixi^
l, ixo^
l
4646 double precision,
intent(in) :: w(ixi^s, nw)
4647 double precision :: gamma_a(ixo^s)
4649 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4650 gamma_a = sqrt(gamma_a)
4651 end function mhd_gamma_alfven
4655 integer,
intent(in) :: ixi^
l, ixo^
l
4656 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4657 double precision,
intent(out) :: rho(ixi^s)
4662 rho(ixo^s) = w(ixo^s,
rho_)
4668 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4671 integer,
intent(in) :: ixi^
l,ixo^
l, ie
4672 double precision,
intent(inout) :: w(ixi^s,1:nw)
4673 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4674 character(len=*),
intent(in) :: subname
4676 double precision :: rho(ixi^s)
4678 logical :: flag(ixi^s,1:nw)
4682 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4683 flag(ixo^s,ie)=.true.
4685 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4687 if(any(flag(ixo^s,ie)))
then
4691 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4694 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4700 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4703 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4709 end subroutine mhd_handle_small_ei
4711 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
4715 integer,
intent(in) :: ixi^
l, ixo^
l
4716 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4717 double precision,
intent(inout) :: w(ixi^s,1:nw)
4719 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
4728 end subroutine mhd_update_temperature
4731 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4734 integer,
intent(in) :: ixi^
l, ixo^
l
4735 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4736 double precision,
intent(inout) :: w(ixi^s,1:nw)
4738 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4750 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4755 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4758 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4764 if(total_energy)
then
4767 b(ixo^s,:)=wct(ixo^s,mag(:))
4776 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4779 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4783 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4787 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
4792 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4797 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
4799 end subroutine add_source_b0split
4802 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4806 integer,
intent(in) :: ixi^
l, ixo^
l
4807 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4808 double precision,
intent(inout) :: w(ixi^s,1:nw)
4809 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4811 double precision :: v(ixi^s,1:3),e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
4812 integer :: idir, idirmin, ix^
d
4814 {
do ix^db=iximin^db,iximax^db\}
4817 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
4818 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
4819 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
4824 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
4832 call divvector(e,ixi^l,ixo^l,dive)
4834 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
4836 call cross_product(ixi^l,ixo^l,e,curle,v)
4839 {
do ix^db=ixomin^db,ixomax^db\}
4840 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)+qdt*(inv_squared_c0-inv_squared_c)*&
4841 (e(ix^d,^
c)*dive(ix^d)-v(ix^d,^
c))\
4844 end subroutine add_source_semirelativistic
4847 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4851 integer,
intent(in) :: ixi^
l, ixo^
l
4852 double precision,
intent(in) :: qdt
4853 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4854 double precision,
intent(inout) :: w(ixi^s,1:nw)
4855 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
4857 double precision :: divv(ixi^s), tmp
4869 {
do ix^db=ixomin^db,ixomax^db\}
4871 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
4872 if(w(ix^
d,
e_)<small_e)
then
4877 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
4880 if(fix_small_values)
then
4881 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
4883 end subroutine add_source_internal_e
4886 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4891 integer,
intent(in) :: ixi^
l, ixo^
l
4892 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4893 double precision,
intent(inout) :: w(ixi^s,1:nw)
4894 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4896 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
4897 double precision :: current(ixi^s,7-2*
ndir:3)
4898 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
4899 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
4900 integer :: idir, idirmin, idims, ix^
d
4905 b(ixo^s, idir) = wct(ixo^s,mag(idir))
4913 j(ixo^s,idir)=current(ixo^s,idir)
4931 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
4937 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
4944 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)
4948 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4952 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
4953 tmp(ixo^s)=sqrt(b2(ixo^s))
4954 where(tmp(ixo^s)>smalldouble)
4955 tmp(ixo^s)=1.d0/tmp(ixo^s)
4961 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
4966 {
do ix^db=ixomin^db,ixomax^db\}
4968 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
4970 b2(ix^
d)=vaoc/(1.d0+vaoc)
4973 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
4976 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
4980 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
4983 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
4984 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
4987 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
4990 end subroutine add_source_hydrodynamic_e
4996 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5001 integer,
intent(in) :: ixi^
l, ixo^
l
5002 double precision,
intent(in) :: qdt
5003 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5004 double precision,
intent(inout) :: w(ixi^s,1:nw)
5005 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5006 integer :: lxo^
l, kxo^
l
5008 double precision :: tmp(ixi^s),tmp2(ixi^s)
5011 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5012 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5021 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5022 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5029 gradeta(ixo^s,1:
ndim)=zero
5035 gradeta(ixo^s,idim)=tmp(ixo^s)
5042 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5049 tmp2(ixi^s)=bf(ixi^s,idir)
5051 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5052 jxo^
l=ixo^
l+
kr(idim,^
d);
5053 hxo^
l=ixo^
l-
kr(idim,^
d);
5054 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5055 tmp(ixo^s)=tmp(ixo^s)+&
5056 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5061 tmp2(ixi^s)=bf(ixi^s,idir)
5063 jxo^
l=ixo^
l+
kr(idim,^
d);
5064 hxo^
l=ixo^
l-
kr(idim,^
d);
5065 tmp(ixo^s)=tmp(ixo^s)+&
5066 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5071 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5075 do jdir=1,
ndim;
do kdir=idirmin,3
5076 if (
lvc(idir,jdir,kdir)/=0)
then
5077 if (
lvc(idir,jdir,kdir)==1)
then
5078 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5080 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5087 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5088 if(total_energy)
then
5089 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5095 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5098 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5100 end subroutine add_source_res1
5104 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5109 integer,
intent(in) :: ixi^
l, ixo^
l
5110 double precision,
intent(in) :: qdt
5111 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5112 double precision,
intent(inout) :: w(ixi^s,1:nw)
5115 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5116 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5117 integer :: ixa^
l,idir,idirmin,idirmin1
5121 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5122 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5132 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5137 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5146 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5149 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5154 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5156 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5158 if(total_energy)
then
5161 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5162 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5165 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5169 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5170 end subroutine add_source_res2
5174 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5178 integer,
intent(in) :: ixi^
l, ixo^
l
5179 double precision,
intent(in) :: qdt
5180 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5181 double precision,
intent(inout) :: w(ixi^s,1:nw)
5183 double precision :: current(ixi^s,7-2*
ndir:3)
5184 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5185 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5188 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5189 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5192 tmpvec(ixa^s,1:
ndir)=zero
5194 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5198 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5201 tmpvec(ixa^s,1:
ndir)=zero
5202 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5206 tmpvec2(ixa^s,1:
ndir)=zero
5207 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5210 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5213 if(total_energy)
then
5216 tmpvec2(ixa^s,1:
ndir)=zero
5217 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5218 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5219 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5220 end do;
end do;
end do
5222 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5223 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5226 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5228 end subroutine add_source_hyperres
5230 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5237 integer,
intent(in) :: ixi^
l, ixo^
l
5238 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5239 double precision,
intent(inout) :: w(ixi^s,1:nw)
5241 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5262 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5265 if(total_energy)
then
5274 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5283 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5287 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5289 end subroutine add_source_glm
5292 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5295 integer,
intent(in) :: ixi^
l, ixo^
l
5296 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5297 double precision,
intent(inout) :: w(ixi^s,1:nw)
5299 double precision :: divb(ixi^s), ba(1:
ndir)
5300 integer :: idir, ix^
d
5306 {
do ix^db=ixomin^db,ixomax^db\}
5311 if (total_energy)
then
5317 {
do ix^db=ixomin^db,ixomax^db\}
5319 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5321 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5322 if (total_energy)
then
5324 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5329 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5331 end subroutine add_source_powel
5333 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5338 integer,
intent(in) :: ixi^
l, ixo^
l
5339 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5340 double precision,
intent(inout) :: w(ixi^s,1:nw)
5342 double precision :: divb(ixi^s)
5343 integer :: idir, ix^
d
5348 {
do ix^db=ixomin^db,ixomax^db\}
5353 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5355 end subroutine add_source_janhunen
5357 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5362 integer,
intent(in) :: ixi^
l, ixo^
l
5363 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5364 double precision,
intent(inout) :: w(ixi^s,1:nw)
5366 double precision :: divb(ixi^s),graddivb(ixi^s)
5367 integer :: idim, idir, ixp^
l, i^
d, iside
5368 logical,
dimension(-1:1^D&) :: leveljump
5376 if(i^
d==0|.and.) cycle
5377 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5378 leveljump(i^
d)=.true.
5380 leveljump(i^
d)=.false.
5389 i^dd=kr(^dd,^d)*(2*iside-3);
5390 if (leveljump(i^dd))
then
5392 ixpmin^d=ixomin^d-i^d
5394 ixpmax^d=ixomax^d-i^d
5405 select case(typegrad)
5407 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5409 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
5413 if (slab_uniform)
then
5414 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5416 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5417 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5420 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
5422 if (typedivbdiff==
'all' .and. total_energy)
then
5424 w(ixp^s,
e_)=w(ixp^s,
e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
5428 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5430 end subroutine add_source_linde
5437 integer,
intent(in) :: ixi^
l, ixo^
l
5438 double precision,
intent(in) :: w(ixi^s,1:nw)
5439 double precision :: divb(ixi^s), dsurface(ixi^s)
5441 double precision :: invb(ixo^s)
5442 integer :: ixa^
l,idims
5444 call get_divb(w,ixi^
l,ixo^
l,divb)
5446 where(invb(ixo^s)/=0.d0)
5447 invb(ixo^s)=1.d0/invb(ixo^s)
5450 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5452 ixamin^
d=ixomin^
d-1;
5453 ixamax^
d=ixomax^
d-1;
5454 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5456 ixa^
l=ixo^
l-
kr(idims,^
d);
5457 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5459 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5460 block%dvolume(ixo^s)/dsurface(ixo^s)
5471 integer,
intent(in) :: ixo^
l, ixi^
l
5472 double precision,
intent(in) :: w(ixi^s,1:nw)
5473 integer,
intent(out) :: idirmin
5476 double precision :: current(ixi^s,7-2*
ndir:3)
5477 integer :: idir, idirmin0
5483 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5484 block%J0(ixo^s,idirmin0:3)
5488 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5496 integer,
intent(in) :: ixi^
l, ixo^
l
5497 double precision,
intent(inout) :: dtnew
5498 double precision,
intent(in) ::
dx^
d
5499 double precision,
intent(in) :: w(ixi^s,1:nw)
5500 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5502 double precision :: dxarr(
ndim)
5503 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5504 integer :: idirmin,idim
5518 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5521 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5547 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5554 end subroutine mhd_get_dt
5557 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5562 integer,
intent(in) :: ixi^
l, ixo^
l
5563 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5564 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5566 double precision :: tmp,tmp1,invr,cot
5568 integer :: mr_,mphi_
5569 integer :: br_,bphi_
5572 br_=mag(1); bphi_=mag(1)-1+
phi_
5577 {
do ix^db=ixomin^db,ixomax^db\}
5580 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5585 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5590 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5591 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5592 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5593 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5594 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5596 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5597 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5598 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5601 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5606 {
do ix^db=ixomin^db,ixomax^db\}
5608 if(local_timestep)
then
5609 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5614 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5620 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5623 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5624 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5628 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5634 cot=1.d0/tan(x(ix^d,2))
5638 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5639 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5641 if(.not.stagger_grid)
then
5642 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5644 tmp=tmp+wprim(ix^d,
psi_)*cot
5646 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5651 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5652 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5653 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5655 if(.not.stagger_grid)
then
5656 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5658 tmp=tmp+wprim(ix^d,
psi_)*cot
5660 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5663 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5664 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5665 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5666 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5667 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5669 if(.not.stagger_grid)
then
5670 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5671 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5672 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5673 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5674 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5681 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5684 end subroutine mhd_add_source_geom
5687 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5692 integer,
intent(in) :: ixi^
l, ixo^
l
5693 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5694 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5696 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
5698 integer :: mr_,mphi_
5699 integer :: br_,bphi_
5702 br_=mag(1); bphi_=mag(1)-1+
phi_
5707 {
do ix^db=ixomin^db,ixomax^db\}
5710 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5721 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
5722 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
5723 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5728 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5734 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
5735 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
5736 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
5737 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5738 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5739 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
5741 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5742 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5743 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5746 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
5747 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
5752 {
do ix^db=ixomin^db,ixomax^db\}
5754 if(local_timestep)
then
5755 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
5761 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5762 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5763 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5767 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5774 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
5780 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
5783 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
5784 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
5785 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
5789 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
5795 cot=1.d0/tan(x(ix^d,2))
5799 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
5800 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
5802 if(.not.stagger_grid)
then
5803 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5805 tmp=tmp+wprim(ix^d,
psi_)*cot
5807 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5813 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
5814 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
5815 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
5816 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
5818 if(.not.stagger_grid)
then
5819 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5821 tmp=tmp+wprim(ix^d,
psi_)*cot
5823 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5826 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
5827 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
5828 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5829 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
5830 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
5831 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
5832 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
5834 if(.not.stagger_grid)
then
5835 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
5836 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5837 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5838 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5839 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5846 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5849 end subroutine mhd_add_source_geom_semirelati
5852 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5857 integer,
intent(in) :: ixi^
l, ixo^
l
5858 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5859 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5861 double precision :: tmp,tmp1,tmp2,invr,cot
5863 integer :: mr_,mphi_
5864 integer :: br_,bphi_
5867 br_=mag(1); bphi_=mag(1)-1+
phi_
5872 {
do ix^db=ixomin^db,ixomax^db\}
5875 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5880 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5885 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5886 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5887 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5888 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5889 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5891 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5892 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5893 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5896 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5901 {
do ix^db=ixomin^db,ixomax^db\}
5903 if(local_timestep)
then
5904 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5908 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5909 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
5912 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5916 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5917 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
5918 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
5920 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5921 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5926 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5932 cot=1.d0/tan(x(ix^d,2))
5937 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5938 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
5939 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
5941 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5942 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5945 if(.not.stagger_grid)
then
5947 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
5948 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
5950 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5953 tmp=tmp+wprim(ix^d,
psi_)*cot
5955 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5961 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5962 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
5963 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
5964 +(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)
5966 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5967 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5968 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5971 if(.not.stagger_grid)
then
5973 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
5974 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
5976 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5979 tmp=tmp+wprim(ix^d,
psi_)*cot
5981 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5985 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5986 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5987 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5988 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
5989 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
5990 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5991 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
5992 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
5993 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
5995 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5996 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5997 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5998 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5999 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6002 if(.not.stagger_grid)
then
6004 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6005 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6006 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6007 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6008 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6009 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6010 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6011 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6012 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6014 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6015 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6016 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6017 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6018 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6026 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6029 end subroutine mhd_add_source_geom_split
6034 integer,
intent(in) :: ixi^
l, ixo^
l
6035 double precision,
intent(in) :: w(ixi^s, nw)
6036 double precision :: mge(ixo^s)
6039 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6041 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6045 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6048 integer,
intent(in) :: ixi^
l, ixo^
l
6049 double precision,
intent(in) :: w(ixi^s,nw)
6050 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6051 double precision,
intent(inout) :: vhall(ixi^s,1:3)
6053 double precision :: current(ixi^s,7-2*
ndir:3)
6054 double precision :: rho(ixi^s)
6055 integer :: idir, idirmin, ix^
d
6060 do idir = idirmin, 3
6061 {
do ix^db=ixomin^db,ixomax^db\}
6062 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6066 end subroutine mhd_getv_hall
6068 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6071 integer,
intent(in) :: ixi^
l, ixo^
l
6072 double precision,
intent(in) :: w(ixi^s,nw)
6073 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6074 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6077 double precision :: current(ixi^s,7-2*
ndir:3)
6078 integer :: idir, idirmin
6085 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6086 do idir = idirmin, 3
6090 end subroutine mhd_get_jambi
6092 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6095 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6096 double precision,
intent(in) :: qt
6097 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6098 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6101 double precision :: db(ixo^s), dpsi(ixo^s)
6105 {
do ix^db=ixomin^db,ixomax^db\}
6106 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6107 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6108 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6109 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6118 {
do ix^db=ixomin^db,ixomax^db\}
6119 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6120 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6121 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6122 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6123 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6125 if(total_energy)
then
6126 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6127 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6129 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6131 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6134 if(total_energy)
then
6135 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6136 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6141 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6143 end subroutine mhd_modify_wlr
6145 subroutine mhd_boundary_adjust(igrid,psb)
6147 integer,
intent(in) :: igrid
6150 integer :: ib, idims, iside, ixo^
l, i^
d
6159 i^
d=
kr(^
d,idims)*(2*iside-3);
6160 if (neighbor_type(i^
d,igrid)/=1) cycle
6161 ib=(idims-1)*2+iside
6179 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6184 end subroutine mhd_boundary_adjust
6186 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6189 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6190 double precision,
intent(inout) :: w(ixg^s,1:nw)
6191 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6193 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6194 integer :: ix^
d,ixf^
l
6207 do ix1=ixfmax1,ixfmin1,-1
6208 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6209 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6210 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6213 do ix1=ixfmax1,ixfmin1,-1
6214 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6215 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6216 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6217 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6218 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6219 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6220 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6234 do ix1=ixfmax1,ixfmin1,-1
6235 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6236 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6237 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6238 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6239 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6240 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6243 do ix1=ixfmax1,ixfmin1,-1
6244 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6245 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6246 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6247 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6248 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6249 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6250 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6251 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6252 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6253 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6254 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6255 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6256 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6257 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6258 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6259 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6260 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6261 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6275 do ix1=ixfmin1,ixfmax1
6276 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6277 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6278 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6281 do ix1=ixfmin1,ixfmax1
6282 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6283 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6284 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6285 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6286 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6287 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6288 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6302 do ix1=ixfmin1,ixfmax1
6303 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6304 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6305 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6306 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6307 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6308 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6311 do ix1=ixfmin1,ixfmax1
6312 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6313 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6314 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6315 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6316 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6317 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6318 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6319 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6320 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6321 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6322 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6323 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6324 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6325 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6326 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6327 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6328 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6329 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6343 do ix2=ixfmax2,ixfmin2,-1
6344 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6345 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6346 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6349 do ix2=ixfmax2,ixfmin2,-1
6350 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6351 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6352 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6353 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6354 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6355 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6356 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6370 do ix2=ixfmax2,ixfmin2,-1
6371 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6372 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6373 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6374 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6375 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6376 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6379 do ix2=ixfmax2,ixfmin2,-1
6380 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6381 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6382 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6383 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6384 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6385 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6386 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6387 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6388 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6389 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6390 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6391 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6392 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6393 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6394 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6395 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6396 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6397 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6411 do ix2=ixfmin2,ixfmax2
6412 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6413 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6414 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6417 do ix2=ixfmin2,ixfmax2
6418 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6419 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6420 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6421 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6422 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6423 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6424 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6438 do ix2=ixfmin2,ixfmax2
6439 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6440 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6441 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6442 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6443 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6444 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6447 do ix2=ixfmin2,ixfmax2
6448 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6449 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6450 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6451 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6452 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6453 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6454 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6455 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6456 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6457 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6458 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6459 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6460 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6461 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6462 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6463 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6464 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6465 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6482 do ix3=ixfmax3,ixfmin3,-1
6483 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6484 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6485 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6486 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6487 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6488 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6491 do ix3=ixfmax3,ixfmin3,-1
6492 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6493 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6494 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6495 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6496 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6497 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6498 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6499 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6500 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6501 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6502 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6503 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6504 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6505 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6506 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6507 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6508 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6509 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6524 do ix3=ixfmin3,ixfmax3
6525 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6526 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6527 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6528 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6529 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6530 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6533 do ix3=ixfmin3,ixfmax3
6534 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6535 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6536 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6537 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6538 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6539 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6540 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6541 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6542 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6543 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6544 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6545 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6546 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6547 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6548 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6549 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6550 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6551 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6557 call mpistop(
"Special boundary is not defined for this region")
6560 end subroutine fixdivb_boundary
6569 double precision,
intent(in) :: qdt
6570 double precision,
intent(in) :: qt
6571 logical,
intent(inout) :: active
6574 integer,
parameter :: max_its = 50
6575 double precision :: residual_it(max_its), max_divb
6576 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6577 double precision :: res
6578 double precision,
parameter :: max_residual = 1
d-3
6579 double precision,
parameter :: residual_reduction = 1
d-10
6580 integer :: iigrid, igrid
6581 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6584 mg%operator_type = mg_laplacian
6592 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6593 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6596 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6597 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6599 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6600 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6603 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6604 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6608 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6609 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6610 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6618 do iigrid = 1, igridstail
6619 igrid = igrids(iigrid);
6622 lvl =
mg%boxes(id)%lvl
6623 nc =
mg%box_size_lvl(lvl)
6629 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6631 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6632 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6637 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6640 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6641 if (
mype == 0) print *,
"iteration vs residual"
6644 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6645 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6646 if (residual_it(n) < residual_reduction * max_divb)
exit
6648 if (
mype == 0 .and. n > max_its)
then
6649 print *,
"divb_multigrid warning: not fully converged"
6650 print *,
"current amplitude of divb: ", residual_it(max_its)
6651 print *,
"multigrid smallest grid: ", &
6652 mg%domain_size_lvl(:,
mg%lowest_lvl)
6653 print *,
"note: smallest grid ideally has <= 8 cells"
6654 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6655 print *,
"note: dx/dy/dz should be similar"
6659 call mg_fas_vcycle(
mg, max_res=res)
6660 if (res < max_residual)
exit
6662 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6667 do iigrid = 1, igridstail
6668 igrid = igrids(iigrid);
6677 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6681 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6683 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
6685 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6698 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6699 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6702 if(total_energy)
then
6704 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6707 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6716 subroutine mhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
6719 integer,
intent(in) :: ixi^
l, ixo^
l
6720 double precision,
intent(in) :: qt,qdt
6722 double precision,
intent(in) :: wprim(ixi^s,1:nw)
6723 type(state) :: sct, s
6724 type(ct_velocity) :: vcts
6725 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6726 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6730 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
6732 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
6734 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
6736 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6739 end subroutine mhd_update_faces
6742 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
6746 integer,
intent(in) :: ixi^
l, ixo^
l
6747 double precision,
intent(in) :: qt, qdt
6748 type(state) :: sct, s
6749 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6750 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6752 double precision :: circ(ixi^s,1:
ndim)
6754 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6755 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
6756 integer :: idim1,idim2,idir,iwdim1,iwdim2
6758 associate(bfaces=>s%ws,x=>s%x)
6765 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6772 i1kr^
d=
kr(idim1,^
d);
6775 i2kr^
d=
kr(idim2,^
d);
6778 if (
lvc(idim1,idim2,idir)==1)
then
6780 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6782 {
do ix^db=ixcmin^db,ixcmax^db\}
6783 fe(ix^
d,idir)=quarter*&
6784 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
6785 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
6787 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
6792 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
6800 if(
associated(usr_set_electric_field)) &
6801 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6803 circ(ixi^s,1:ndim)=zero
6808 ixcmin^d=ixomin^d-kr(idim1,^d);
6810 ixa^l=ixc^l-kr(idim2,^d);
6813 if(lvc(idim1,idim2,idir)==1)
then
6815 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6818 else if(lvc(idim1,idim2,idir)==-1)
then
6820 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6827 where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
6828 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6830 circ(ixc^s,idim1)=zero
6833 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6838 end subroutine update_faces_average
6841 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
6846 integer,
intent(in) :: ixi^
l, ixo^
l
6847 double precision,
intent(in) :: qt, qdt
6849 double precision,
intent(in) :: wp(ixi^s,1:nw)
6850 type(state) :: sct, s
6851 type(ct_velocity) :: vcts
6852 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6853 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6855 double precision :: circ(ixi^s,1:
ndim)
6857 double precision :: ecc(ixi^s,
sdim:3)
6858 double precision :: ein(ixi^s,
sdim:3)
6860 double precision :: el(ixi^s),er(ixi^s)
6862 double precision :: elc,erc
6864 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6866 double precision :: jce(ixi^s,
sdim:3)
6868 double precision :: xs(ixgs^t,1:
ndim)
6869 double precision :: gradi(ixgs^t)
6870 integer :: ixc^
l,ixa^
l
6871 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
6873 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
6876 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6882 {
do ix^db=iximin^db,iximax^db\}
6885 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_)
6886 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_)
6887 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_)
6890 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
6897 {
do ix^db=iximin^db,iximax^db\}
6900 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
6901 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
6902 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6905 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6919 i1kr^d=kr(idim1,^d);
6922 i2kr^d=kr(idim2,^d);
6925 if (lvc(idim1,idim2,idir)==1)
then
6927 ixcmin^d=ixomin^d+kr(idir,^d)-1;
6930 {
do ix^db=ixcmin^db,ixcmax^db\}
6931 fe(ix^d,idir)=quarter*&
6932 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
6933 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
6938 ixamax^d=ixcmax^d+i1kr^d;
6939 {
do ix^db=ixamin^db,ixamax^db\}
6940 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
6941 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
6944 do ix^db=ixcmin^db,ixcmax^db\}
6945 if(vnorm(ix^d,idim1)>0.d0)
then
6947 else if(vnorm(ix^d,idim1)<0.d0)
then
6948 elc=el({ix^d+i1kr^d})
6950 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
6952 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
6954 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
6955 erc=er({ix^d+i1kr^d})
6957 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
6959 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
6964 ixamax^d=ixcmax^d+i2kr^d;
6965 {
do ix^db=ixamin^db,ixamax^db\}
6966 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
6967 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
6970 do ix^db=ixcmin^db,ixcmax^db\}
6971 if(vnorm(ix^d,idim2)>0.d0)
then
6973 else if(vnorm(ix^d,idim2)<0.d0)
then
6974 elc=el({ix^d+i2kr^d})
6976 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
6978 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
6980 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
6981 erc=er({ix^d+i2kr^d})
6983 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
6985 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
6989 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
6994 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7008 if (lvc(idim1,idim2,idir)==0) cycle
7010 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7011 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7014 xs(ixa^s,:)=x(ixa^s,:)
7015 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7016 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7017 if (lvc(idim1,idim2,idir)==1)
then
7018 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7020 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7027 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7029 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7033 {
do ix^db=ixomin^db,ixomax^db\}
7034 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7035 +ein(ix1,ix2-1,ix3-1,idir))
7036 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7037 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7039 else if(idir==2)
then
7040 {
do ix^db=ixomin^db,ixomax^db\}
7041 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7042 +ein(ix1-1,ix2,ix3-1,idir))
7043 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7044 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7047 {
do ix^db=ixomin^db,ixomax^db\}
7048 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7049 +ein(ix1-1,ix2-1,ix3,idir))
7050 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7051 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7057 {
do ix^db=ixomin^db,ixomax^db\}
7058 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7059 +ein(ix1-1,ix2-1,idir))
7060 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7061 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7066 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7072 if(
associated(usr_set_electric_field)) &
7073 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7075 circ(ixi^s,1:ndim)=zero
7080 ixcmin^d=ixomin^d-kr(idim1,^d);
7082 ixa^l=ixc^l-kr(idim2,^d);
7085 if(lvc(idim1,idim2,idir)==1)
then
7087 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7090 else if(lvc(idim1,idim2,idir)==-1)
then
7092 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7099 where(s%surfaceC(ixc^s,idim1) > smalldouble)
7100 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
7102 circ(ixc^s,idim1)=zero
7105 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
7110 end subroutine update_faces_contact
7113 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
7118 integer,
intent(in) :: ixi^
l, ixo^
l
7119 double precision,
intent(in) :: qt, qdt
7120 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7121 type(state) :: sct, s
7122 type(ct_velocity) :: vcts
7124 double precision :: vtill(ixi^s,2)
7125 double precision :: vtilr(ixi^s,2)
7126 double precision :: bfacetot(ixi^s,
ndim)
7127 double precision :: btill(ixi^s,
ndim)
7128 double precision :: btilr(ixi^s,
ndim)
7129 double precision :: cp(ixi^s,2)
7130 double precision :: cm(ixi^s,2)
7131 double precision :: circ(ixi^s,1:
ndim)
7133 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7134 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7135 integer :: idim1,idim2,idir
7137 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7138 cbarmax=>vcts%cbarmax)
7151 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
7167 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7171 idim2=mod(idir+1,3)+1
7173 jxc^
l=ixc^
l+
kr(idim1,^
d);
7174 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7178 vtill(ixi^s,2),vtilr(ixi^s,2))
7181 vtill(ixi^s,1),vtilr(ixi^s,1))
7187 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7188 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7190 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7191 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7194 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7197 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7201 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7202 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7204 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7205 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7209 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7210 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7211 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7212 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7213 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7214 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7215 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7216 /(cp(ixc^s,2)+cm(ixc^s,2))
7219 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7223 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7237 circ(ixi^s,1:
ndim)=zero
7242 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7246 if(
lvc(idim1,idim2,idir)/=0)
then
7247 hxc^
l=ixc^
l-
kr(idim2,^
d);
7249 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7250 +
lvc(idim1,idim2,idir)&
7257 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
7258 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
7260 circ(ixc^s,idim1)=zero
7263 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
7267 end subroutine update_faces_hll
7270 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
7275 integer,
intent(in) :: ixi^
l, ixo^
l
7276 type(state),
intent(in) :: sct, s
7278 double precision :: jce(ixi^s,
sdim:3)
7281 double precision :: jcc(ixi^s,7-2*
ndir:3)
7283 double precision :: xs(ixgs^t,1:
ndim)
7285 double precision :: eta(ixi^s)
7286 double precision :: gradi(ixgs^t)
7287 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7289 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7295 if (
lvc(idim1,idim2,idir)==0) cycle
7297 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7298 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7301 xs(ixb^s,:)=x(ixb^s,:)
7302 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7303 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7304 if (
lvc(idim1,idim2,idir)==1)
then
7305 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7307 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7314 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7322 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7323 jcc(ixc^s,idir)=0.d0
7325 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7326 ixamin^
d=ixcmin^
d+ix^
d;
7327 ixamax^
d=ixcmax^
d+ix^
d;
7328 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7330 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7331 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7336 end subroutine get_resistive_electric_field
7339 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7342 integer,
intent(in) :: ixi^
l, ixo^
l
7343 double precision,
intent(in) :: w(ixi^s,1:nw)
7344 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7345 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7347 double precision :: jxbxb(ixi^s,1:3)
7348 integer :: idir,ixa^
l,ixc^
l,ix^
d
7351 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7358 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7361 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7362 ixamin^
d=ixcmin^
d+ix^
d;
7363 ixamax^
d=ixcmax^
d+ix^
d;
7364 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7366 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7369 end subroutine get_ambipolar_electric_field
7375 integer,
intent(in) :: ixo^
l
7385 do ix^db=ixomin^db,ixomax^db\}
7387 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7388 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7389 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7390 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7391 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7392 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7395 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7396 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7397 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7398 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7441 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7442 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7443 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7445 double precision :: adummy(ixis^s,1:3)
7451 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7454 integer,
intent(in) :: ixi^
l, ixo^
l
7455 double precision,
intent(in) :: w(ixi^s,1:nw)
7456 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7457 double precision,
intent(out):: rfactor(ixi^s)
7459 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7463 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)
7465 end subroutine rfactor_from_temperature_ionization
7467 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7469 integer,
intent(in) :: ixi^
l, ixo^
l
7470 double precision,
intent(in) :: w(ixi^s,1:nw)
7471 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7472 double precision,
intent(out):: rfactor(ixi^s)
7476 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.