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 gradientx(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),.false.)
2758 call gradientq(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims))
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)
4499 double precision :: divv(ixi^s)
4516 end subroutine add_pe0_divv
4518 subroutine get_tau(ixI^L,ixO^L,w,Te,tau,sigT5)
4520 integer,
intent(in) :: ixi^
l, ixo^
l
4521 double precision,
dimension(ixI^S,1:nw),
intent(in) :: w
4522 double precision,
dimension(ixI^S),
intent(in) :: te
4523 double precision,
dimension(ixI^S),
intent(out) :: tau,sigt5
4525 double precision :: dxmin,taumin
4526 double precision,
dimension(ixI^S) :: sigt7,eint
4534 sigt7(ixo^s)=sigt5(ixo^s)*
block%wextra(ixo^s,
tcoff_)
4537 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4541 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
4544 tau(ixo^s)=max(taumin*
dt,sigt7(ixo^s)/eint(ixo^s)/
cmax_global**2)
4545 end subroutine get_tau
4547 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4549 integer,
intent(in) :: ixi^
l,ixo^
l
4550 double precision,
intent(in) :: qdt
4551 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4552 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4553 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4555 double precision :: invdx
4556 double precision,
dimension(ixI^S) :: te,tau,sigt,htc_qsrc,tface,r
4557 double precision,
dimension(ixI^S) :: htc_esrc,bsum,bunit
4558 double precision,
dimension(ixI^S,1:ndim) :: btot
4560 integer :: hxc^
l,hxo^
l,ixc^
l,jxc^
l,jxo^
l,kxc^
l
4562 call mhd_get_rfactor(wctprim,x,ixi^
l,ixi^
l,r)
4564 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*w(ixi^s,
rho_))
4565 call get_tau(ixi^
l,ixo^
l,wctprim,te,tau,sigt)
4569 btot(ixo^s,idims)=wct(ixo^s,mag(idims))+
block%B0(ixo^s,idims,0)
4571 btot(ixo^s,idims)=wct(ixo^s,mag(idims))
4574 bsum(ixo^s)=sqrt(sum(btot(ixo^s,:)**2,dim=
ndim+1))+smalldouble
4578 ixcmin^
d=ixomin^
d-
kr(idims,^
d);ixcmax^
d=ixomax^
d;
4579 jxc^
l=ixc^
l+
kr(idims,^
d);
4580 kxc^
l=jxc^
l+
kr(idims,^
d);
4581 hxc^
l=ixc^
l-
kr(idims,^
d);
4582 hxo^
l=ixo^
l-
kr(idims,^
d);
4583 tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
4584 bunit(ixo^s)=btot(ixo^s,idims)/bsum(ixo^s)
4585 htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*bunit(ixo^s)*(tface(ixo^s)-tface(hxo^s))*invdx
4587 htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,
q_))/tau(ixo^s)
4588 w(ixo^s,
q_)=w(ixo^s,
q_)-qdt*htc_qsrc(ixo^s)
4589 end subroutine add_hypertc_source
4592 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4594 integer,
intent(in) :: ixi^
l, ixo^
l
4595 double precision,
intent(in) :: w(ixi^s,1:nw)
4596 double precision,
intent(inout) :: jxb(ixi^s,3)
4597 double precision :: a(ixi^s,3),
b(ixi^s,3)
4599 double precision :: current(ixi^s,7-2*
ndir:3)
4600 integer :: idir, idirmin
4605 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4609 b(ixo^s, idir) = w(ixo^s,mag(idir))
4618 a(ixo^s,idir)=current(ixo^s,idir)
4622 end subroutine get_lorentz_force
4626 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4628 integer,
intent(in) :: ixi^
l, ixo^
l
4629 double precision,
intent(in) :: w(ixi^s, nw)
4630 double precision,
intent(out) :: gamma_a2(ixo^s)
4631 double precision :: rho(ixi^s)
4637 rho(ixo^s) = w(ixo^s,
rho_)
4640 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4641 end subroutine mhd_gamma2_alfven
4645 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4647 integer,
intent(in) :: ixi^
l, ixo^
l
4648 double precision,
intent(in) :: w(ixi^s, nw)
4649 double precision :: gamma_a(ixo^s)
4651 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4652 gamma_a = sqrt(gamma_a)
4653 end function mhd_gamma_alfven
4657 integer,
intent(in) :: ixi^
l, ixo^
l
4658 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4659 double precision,
intent(out) :: rho(ixi^s)
4664 rho(ixo^s) = w(ixo^s,
rho_)
4670 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4673 integer,
intent(in) :: ixi^
l,ixo^
l, ie
4674 double precision,
intent(inout) :: w(ixi^s,1:nw)
4675 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4676 character(len=*),
intent(in) :: subname
4678 double precision :: rho(ixi^s)
4680 logical :: flag(ixi^s,1:nw)
4684 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4685 flag(ixo^s,ie)=.true.
4687 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4689 if(any(flag(ixo^s,ie)))
then
4693 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4696 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4702 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4705 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4711 end subroutine mhd_handle_small_ei
4713 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
4717 integer,
intent(in) :: ixi^
l, ixo^
l
4718 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4719 double precision,
intent(inout) :: w(ixi^s,1:nw)
4721 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
4730 end subroutine mhd_update_temperature
4733 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4736 integer,
intent(in) :: ixi^
l, ixo^
l
4737 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4738 double precision,
intent(inout) :: w(ixi^s,1:nw)
4740 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4752 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4757 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4760 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4766 if(total_energy)
then
4769 b(ixo^s,:)=wct(ixo^s,mag(:))
4778 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4781 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4785 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4789 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
4794 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4799 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
4801 end subroutine add_source_b0split
4804 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4808 integer,
intent(in) :: ixi^
l, ixo^
l
4809 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4810 double precision,
intent(inout) :: w(ixi^s,1:nw)
4811 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4813 double precision :: v(ixi^s,1:3),e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
4814 integer :: idir, idirmin, ix^
d
4816 {
do ix^db=iximin^db,iximax^db\}
4819 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
4820 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
4821 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
4826 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
4834 call divvector(e,ixi^l,ixo^l,dive)
4836 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
4838 call cross_product(ixi^l,ixo^l,e,curle,v)
4841 {
do ix^db=ixomin^db,ixomax^db\}
4842 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)+qdt*(inv_squared_c0-inv_squared_c)*&
4843 (e(ix^d,^
c)*dive(ix^d)-v(ix^d,^
c))\
4846 end subroutine add_source_semirelativistic
4849 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4853 integer,
intent(in) :: ixi^
l, ixo^
l
4854 double precision,
intent(in) :: qdt
4855 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4856 double precision,
intent(inout) :: w(ixi^s,1:nw)
4857 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
4859 double precision :: divv(ixi^s), tmp
4864 call divvector(wctprim(ixi^s,
mom(:)),ixi^
l,ixo^
l,divv,sixthorder=.true.)
4866 call divvector(wctprim(ixi^s,
mom(:)),ixi^
l,ixo^
l,divv,fourthorder=.true.)
4871 {
do ix^db=ixomin^db,ixomax^db\}
4873 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
4874 if(w(ix^
d,
e_)<small_e)
then
4879 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
4882 if(fix_small_values)
then
4883 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
4885 end subroutine add_source_internal_e
4888 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4893 integer,
intent(in) :: ixi^
l, ixo^
l
4894 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4895 double precision,
intent(inout) :: w(ixi^s,1:nw)
4896 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
4898 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
4899 double precision :: current(ixi^s,7-2*
ndir:3)
4900 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
4901 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
4902 integer :: idir, idirmin, idims, ix^
d
4907 b(ixo^s, idir) = wct(ixo^s,mag(idir))
4915 j(ixo^s,idir)=current(ixo^s,idir)
4933 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
4939 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
4946 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)
4950 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4954 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
4955 tmp(ixo^s)=sqrt(b2(ixo^s))
4956 where(tmp(ixo^s)>smalldouble)
4957 tmp(ixo^s)=1.d0/tmp(ixo^s)
4963 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
4968 {
do ix^db=ixomin^db,ixomax^db\}
4970 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
4972 b2(ix^
d)=vaoc/(1.d0+vaoc)
4975 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
4978 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
4982 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
4985 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
4986 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
4989 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
4992 end subroutine add_source_hydrodynamic_e
4998 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5003 integer,
intent(in) :: ixi^
l, ixo^
l
5004 double precision,
intent(in) :: qdt
5005 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5006 double precision,
intent(inout) :: w(ixi^s,1:nw)
5007 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5008 integer :: lxo^
l, kxo^
l
5010 double precision :: tmp(ixi^s),tmp2(ixi^s)
5013 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5014 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5023 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5024 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5031 gradeta(ixo^s,1:
ndim)=zero
5037 gradeta(ixo^s,idim)=tmp(ixo^s)
5044 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5051 tmp2(ixi^s)=bf(ixi^s,idir)
5053 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5054 jxo^
l=ixo^
l+
kr(idim,^
d);
5055 hxo^
l=ixo^
l-
kr(idim,^
d);
5056 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5057 tmp(ixo^s)=tmp(ixo^s)+&
5058 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5063 tmp2(ixi^s)=bf(ixi^s,idir)
5065 jxo^
l=ixo^
l+
kr(idim,^
d);
5066 hxo^
l=ixo^
l-
kr(idim,^
d);
5067 tmp(ixo^s)=tmp(ixo^s)+&
5068 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5073 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5077 do jdir=1,
ndim;
do kdir=idirmin,3
5078 if (
lvc(idir,jdir,kdir)/=0)
then
5079 if (
lvc(idir,jdir,kdir)==1)
then
5080 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5082 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5089 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5090 if(total_energy)
then
5091 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5097 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5100 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5102 end subroutine add_source_res1
5106 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5111 integer,
intent(in) :: ixi^
l, ixo^
l
5112 double precision,
intent(in) :: qdt
5113 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5114 double precision,
intent(inout) :: w(ixi^s,1:nw)
5117 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5118 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5119 integer :: ixa^
l,idir,idirmin,idirmin1
5123 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5124 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5134 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5139 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5148 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5151 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5156 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5158 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5160 if(total_energy)
then
5163 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5164 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5167 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5171 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5172 end subroutine add_source_res2
5176 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5180 integer,
intent(in) :: ixi^
l, ixo^
l
5181 double precision,
intent(in) :: qdt
5182 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5183 double precision,
intent(inout) :: w(ixi^s,1:nw)
5185 double precision :: current(ixi^s,7-2*
ndir:3)
5186 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5187 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5190 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5191 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5194 tmpvec(ixa^s,1:
ndir)=zero
5196 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5200 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5203 tmpvec(ixa^s,1:
ndir)=zero
5204 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5208 tmpvec2(ixa^s,1:
ndir)=zero
5209 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5212 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5215 if(total_energy)
then
5218 tmpvec2(ixa^s,1:
ndir)=zero
5219 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5220 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5221 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5222 end do;
end do;
end do
5224 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5225 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5228 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5230 end subroutine add_source_hyperres
5232 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5239 integer,
intent(in) :: ixi^
l, ixo^
l
5240 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5241 double precision,
intent(inout) :: w(ixi^s,1:nw)
5243 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5264 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5267 if(total_energy)
then
5276 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5285 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5289 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5291 end subroutine add_source_glm
5294 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5297 integer,
intent(in) :: ixi^
l, ixo^
l
5298 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5299 double precision,
intent(inout) :: w(ixi^s,1:nw)
5301 double precision :: divb(ixi^s), ba(1:
ndir)
5302 integer :: idir, ix^
d
5308 {
do ix^db=ixomin^db,ixomax^db\}
5313 if (total_energy)
then
5319 {
do ix^db=ixomin^db,ixomax^db\}
5321 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5323 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5324 if (total_energy)
then
5326 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5331 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5333 end subroutine add_source_powel
5335 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5340 integer,
intent(in) :: ixi^
l, ixo^
l
5341 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5342 double precision,
intent(inout) :: w(ixi^s,1:nw)
5344 double precision :: divb(ixi^s)
5345 integer :: idir, ix^
d
5350 {
do ix^db=ixomin^db,ixomax^db\}
5355 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5357 end subroutine add_source_janhunen
5359 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5364 integer,
intent(in) :: ixi^
l, ixo^
l
5365 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5366 double precision,
intent(inout) :: w(ixi^s,1:nw)
5368 double precision :: divb(ixi^s),graddivb(ixi^s)
5369 integer :: idim, idir, ixp^
l, i^
d, iside
5370 logical,
dimension(-1:1^D&) :: leveljump
5378 if(i^
d==0|.and.) cycle
5379 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5380 leveljump(i^
d)=.true.
5382 leveljump(i^
d)=.false.
5391 i^dd=kr(^dd,^d)*(2*iside-3);
5392 if (leveljump(i^dd))
then
5394 ixpmin^d=ixomin^d-i^d
5396 ixpmax^d=ixomax^d-i^d
5407 select case(typegrad)
5409 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5411 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
5415 if (slab_uniform)
then
5416 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5418 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5419 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5422 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
5424 if (typedivbdiff==
'all' .and. total_energy)
then
5426 w(ixp^s,
e_)=w(ixp^s,
e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
5430 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5432 end subroutine add_source_linde
5439 integer,
intent(in) :: ixi^
l, ixo^
l
5440 double precision,
intent(in) :: w(ixi^s,1:nw)
5441 double precision :: divb(ixi^s), dsurface(ixi^s)
5443 double precision :: invb(ixo^s)
5444 integer :: ixa^
l,idims
5446 call get_divb(w,ixi^
l,ixo^
l,divb)
5448 where(invb(ixo^s)/=0.d0)
5449 invb(ixo^s)=1.d0/invb(ixo^s)
5452 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5454 ixamin^
d=ixomin^
d-1;
5455 ixamax^
d=ixomax^
d-1;
5456 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5458 ixa^
l=ixo^
l-
kr(idims,^
d);
5459 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5461 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5462 block%dvolume(ixo^s)/dsurface(ixo^s)
5473 integer,
intent(in) :: ixo^
l, ixi^
l
5474 double precision,
intent(in) :: w(ixi^s,1:nw)
5475 integer,
intent(out) :: idirmin
5478 double precision :: current(ixi^s,7-2*
ndir:3)
5479 integer :: idir, idirmin0
5485 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5486 block%J0(ixo^s,idirmin0:3)
5490 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5498 integer,
intent(in) :: ixi^
l, ixo^
l
5499 double precision,
intent(inout) :: dtnew
5500 double precision,
intent(in) ::
dx^
d
5501 double precision,
intent(in) :: w(ixi^s,1:nw)
5502 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5504 double precision :: dxarr(
ndim)
5505 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5506 integer :: idirmin,idim
5520 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5523 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5549 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5556 end subroutine mhd_get_dt
5559 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5564 integer,
intent(in) :: ixi^
l, ixo^
l
5565 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5566 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5568 double precision :: tmp,tmp1,invr,cot
5570 integer :: mr_,mphi_
5571 integer :: br_,bphi_
5574 br_=mag(1); bphi_=mag(1)-1+
phi_
5579 {
do ix^db=ixomin^db,ixomax^db\}
5582 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5587 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5592 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5593 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5594 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5595 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5596 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5598 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5599 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5600 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5603 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5608 {
do ix^db=ixomin^db,ixomax^db\}
5610 if(local_timestep)
then
5611 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5616 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5622 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5625 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5626 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5630 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5636 cot=1.d0/tan(x(ix^d,2))
5640 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5641 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5643 if(.not.stagger_grid)
then
5644 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5646 tmp=tmp+wprim(ix^d,
psi_)*cot
5648 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5653 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5654 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5655 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5657 if(.not.stagger_grid)
then
5658 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5660 tmp=tmp+wprim(ix^d,
psi_)*cot
5662 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5665 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5666 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5667 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5668 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5669 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5671 if(.not.stagger_grid)
then
5672 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5673 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5674 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5675 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5676 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5683 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5686 end subroutine mhd_add_source_geom
5689 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5694 integer,
intent(in) :: ixi^
l, ixo^
l
5695 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5696 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5698 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
5700 integer :: mr_,mphi_
5701 integer :: br_,bphi_
5704 br_=mag(1); bphi_=mag(1)-1+
phi_
5709 {
do ix^db=ixomin^db,ixomax^db\}
5712 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5723 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
5724 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
5725 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5730 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5736 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
5737 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
5738 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
5739 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5740 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5741 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
5743 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5744 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5745 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5748 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
5749 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
5754 {
do ix^db=ixomin^db,ixomax^db\}
5756 if(local_timestep)
then
5757 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
5763 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5764 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5765 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5769 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5776 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
5782 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
5785 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
5786 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
5787 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
5791 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
5797 cot=1.d0/tan(x(ix^d,2))
5801 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
5802 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
5804 if(.not.stagger_grid)
then
5805 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5807 tmp=tmp+wprim(ix^d,
psi_)*cot
5809 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5815 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
5816 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
5817 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
5818 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
5820 if(.not.stagger_grid)
then
5821 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5823 tmp=tmp+wprim(ix^d,
psi_)*cot
5825 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
5828 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
5829 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
5830 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5831 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
5832 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
5833 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
5834 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
5836 if(.not.stagger_grid)
then
5837 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
5838 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5839 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5840 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5841 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5848 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5851 end subroutine mhd_add_source_geom_semirelati
5854 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5859 integer,
intent(in) :: ixi^
l, ixo^
l
5860 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5861 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5863 double precision :: tmp,tmp1,tmp2,invr,cot
5865 integer :: mr_,mphi_
5866 integer :: br_,bphi_
5869 br_=mag(1); bphi_=mag(1)-1+
phi_
5874 {
do ix^db=ixomin^db,ixomax^db\}
5877 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5882 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5887 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5888 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5889 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5890 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5891 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5893 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5894 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5895 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5898 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5903 {
do ix^db=ixomin^db,ixomax^db\}
5905 if(local_timestep)
then
5906 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5910 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5911 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
5914 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5918 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5919 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
5920 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
5922 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5923 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5928 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5934 cot=1.d0/tan(x(ix^d,2))
5939 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5940 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
5941 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
5943 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5944 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5947 if(.not.stagger_grid)
then
5949 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
5950 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
5952 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5955 tmp=tmp+wprim(ix^d,
psi_)*cot
5957 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5963 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5964 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
5965 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
5966 +(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)
5968 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5969 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5970 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5973 if(.not.stagger_grid)
then
5975 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
5976 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
5978 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5981 tmp=tmp+wprim(ix^d,
psi_)*cot
5983 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5987 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5988 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5989 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5990 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
5991 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
5992 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5993 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
5994 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
5995 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
5997 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5998 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5999 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6000 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6001 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6004 if(.not.stagger_grid)
then
6006 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6007 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6008 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6009 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6010 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6011 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6012 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6013 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6014 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6016 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6017 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6018 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6019 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6020 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6028 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6031 end subroutine mhd_add_source_geom_split
6036 integer,
intent(in) :: ixi^
l, ixo^
l
6037 double precision,
intent(in) :: w(ixi^s, nw)
6038 double precision :: mge(ixo^s)
6041 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6043 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6047 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6050 integer,
intent(in) :: ixi^
l, ixo^
l
6051 double precision,
intent(in) :: w(ixi^s,nw)
6052 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6053 double precision,
intent(inout) :: vhall(ixi^s,1:3)
6055 double precision :: current(ixi^s,7-2*
ndir:3)
6056 double precision :: rho(ixi^s)
6057 integer :: idir, idirmin, ix^
d
6062 do idir = idirmin, 3
6063 {
do ix^db=ixomin^db,ixomax^db\}
6064 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6068 end subroutine mhd_getv_hall
6070 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6073 integer,
intent(in) :: ixi^
l, ixo^
l
6074 double precision,
intent(in) :: w(ixi^s,nw)
6075 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6076 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6079 double precision :: current(ixi^s,7-2*
ndir:3)
6080 integer :: idir, idirmin
6087 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6088 do idir = idirmin, 3
6092 end subroutine mhd_get_jambi
6094 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6097 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6098 double precision,
intent(in) :: qt
6099 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6100 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6103 double precision :: db(ixo^s), dpsi(ixo^s)
6107 {
do ix^db=ixomin^db,ixomax^db\}
6108 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6109 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6110 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6111 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6120 {
do ix^db=ixomin^db,ixomax^db\}
6121 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6122 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6123 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6124 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6125 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6127 if(total_energy)
then
6128 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6129 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6131 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6133 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6136 if(total_energy)
then
6137 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6138 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6143 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6145 end subroutine mhd_modify_wlr
6147 subroutine mhd_boundary_adjust(igrid,psb)
6149 integer,
intent(in) :: igrid
6152 integer :: ib, idims, iside, ixo^
l, i^
d
6161 i^
d=
kr(^
d,idims)*(2*iside-3);
6162 if (neighbor_type(i^
d,igrid)/=1) cycle
6163 ib=(idims-1)*2+iside
6181 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6186 end subroutine mhd_boundary_adjust
6188 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6191 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6192 double precision,
intent(inout) :: w(ixg^s,1:nw)
6193 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6195 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6196 integer :: ix^
d,ixf^
l
6209 do ix1=ixfmax1,ixfmin1,-1
6210 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6211 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6212 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6215 do ix1=ixfmax1,ixfmin1,-1
6216 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6217 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6218 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6219 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6220 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6221 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6222 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6236 do ix1=ixfmax1,ixfmin1,-1
6237 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6238 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6239 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6240 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6241 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6242 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6245 do ix1=ixfmax1,ixfmin1,-1
6246 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6247 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6248 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6249 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6250 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6251 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6252 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6253 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6254 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6255 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6256 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6257 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6258 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6259 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6260 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6261 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6262 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6263 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6277 do ix1=ixfmin1,ixfmax1
6278 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6279 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6280 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6283 do ix1=ixfmin1,ixfmax1
6284 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6285 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6286 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6287 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6288 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6289 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6290 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6304 do ix1=ixfmin1,ixfmax1
6305 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6306 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6307 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6308 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6309 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6310 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6313 do ix1=ixfmin1,ixfmax1
6314 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6315 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6316 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6317 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6318 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6319 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6320 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6321 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6322 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6323 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6324 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6325 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6326 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6327 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6328 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6329 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6330 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6331 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6345 do ix2=ixfmax2,ixfmin2,-1
6346 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6347 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6348 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6351 do ix2=ixfmax2,ixfmin2,-1
6352 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6353 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6354 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6355 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6356 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6357 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6358 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6372 do ix2=ixfmax2,ixfmin2,-1
6373 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6374 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6375 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6376 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6377 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6378 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6381 do ix2=ixfmax2,ixfmin2,-1
6382 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6383 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6384 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6385 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6386 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6387 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6388 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6389 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6390 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6391 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6392 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6393 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6394 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6395 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6396 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6397 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6398 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6399 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6413 do ix2=ixfmin2,ixfmax2
6414 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6415 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6416 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6419 do ix2=ixfmin2,ixfmax2
6420 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6421 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6422 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6423 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6424 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6425 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6426 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6440 do ix2=ixfmin2,ixfmax2
6441 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6442 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6443 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6444 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6445 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6446 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6449 do ix2=ixfmin2,ixfmax2
6450 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6451 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6452 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6453 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6454 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6455 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6456 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6457 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6458 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6459 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6460 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6461 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6462 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6463 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6464 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6465 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6466 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6467 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6484 do ix3=ixfmax3,ixfmin3,-1
6485 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6486 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6487 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6488 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6489 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6490 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6493 do ix3=ixfmax3,ixfmin3,-1
6494 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6495 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6496 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6497 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6498 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6499 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6500 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6501 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6502 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6503 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6504 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6505 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6506 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6507 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6508 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6509 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6510 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6511 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6526 do ix3=ixfmin3,ixfmax3
6527 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6528 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6529 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6530 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6531 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6532 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6535 do ix3=ixfmin3,ixfmax3
6536 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6537 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6538 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6539 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6540 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6541 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6542 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6543 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6544 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6545 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6546 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6547 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6548 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6549 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6550 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6551 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6552 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6553 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6559 call mpistop(
"Special boundary is not defined for this region")
6562 end subroutine fixdivb_boundary
6571 double precision,
intent(in) :: qdt
6572 double precision,
intent(in) :: qt
6573 logical,
intent(inout) :: active
6576 integer,
parameter :: max_its = 50
6577 double precision :: residual_it(max_its), max_divb
6578 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6579 double precision :: res
6580 double precision,
parameter :: max_residual = 1
d-3
6581 double precision,
parameter :: residual_reduction = 1
d-10
6582 integer :: iigrid, igrid
6583 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6586 mg%operator_type = mg_laplacian
6594 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6595 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6598 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6599 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6601 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6602 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6605 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6606 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6610 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6611 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6612 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6620 do iigrid = 1, igridstail
6621 igrid = igrids(iigrid);
6624 lvl =
mg%boxes(id)%lvl
6625 nc =
mg%box_size_lvl(lvl)
6631 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6633 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6634 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6639 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6642 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6643 if (
mype == 0) print *,
"iteration vs residual"
6646 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6647 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6648 if (residual_it(n) < residual_reduction * max_divb)
exit
6650 if (
mype == 0 .and. n > max_its)
then
6651 print *,
"divb_multigrid warning: not fully converged"
6652 print *,
"current amplitude of divb: ", residual_it(max_its)
6653 print *,
"multigrid smallest grid: ", &
6654 mg%domain_size_lvl(:,
mg%lowest_lvl)
6655 print *,
"note: smallest grid ideally has <= 8 cells"
6656 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6657 print *,
"note: dx/dy/dz should be similar"
6661 call mg_fas_vcycle(
mg, max_res=res)
6662 if (res < max_residual)
exit
6664 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6669 do iigrid = 1, igridstail
6670 igrid = igrids(iigrid);
6679 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6683 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6685 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
6687 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6700 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6701 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6704 if(total_energy)
then
6706 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6709 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6718 subroutine mhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
6721 integer,
intent(in) :: ixi^
l, ixo^
l
6722 double precision,
intent(in) :: qt,qdt
6724 double precision,
intent(in) :: wprim(ixi^s,1:nw)
6725 type(state) :: sct, s
6726 type(ct_velocity) :: vcts
6727 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6728 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6732 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
6734 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
6736 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
6738 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6741 end subroutine mhd_update_faces
6744 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
6748 integer,
intent(in) :: ixi^
l, ixo^
l
6749 double precision,
intent(in) :: qt, qdt
6750 type(state) :: sct, s
6751 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6752 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6754 double precision :: circ(ixi^s,1:
ndim)
6756 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6757 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
6758 integer :: idim1,idim2,idir,iwdim1,iwdim2
6760 associate(bfaces=>s%ws,x=>s%x)
6767 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6774 i1kr^
d=
kr(idim1,^
d);
6777 i2kr^
d=
kr(idim2,^
d);
6780 if (
lvc(idim1,idim2,idir)==1)
then
6782 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6784 {
do ix^db=ixcmin^db,ixcmax^db\}
6785 fe(ix^
d,idir)=quarter*&
6786 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
6787 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
6789 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
6794 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
6802 if(
associated(usr_set_electric_field)) &
6803 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6805 circ(ixi^s,1:ndim)=zero
6810 ixcmin^d=ixomin^d-kr(idim1,^d);
6812 ixa^l=ixc^l-kr(idim2,^d);
6815 if(lvc(idim1,idim2,idir)==1)
then
6817 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6820 else if(lvc(idim1,idim2,idir)==-1)
then
6822 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
6829 where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
6830 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
6832 circ(ixc^s,idim1)=zero
6835 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
6840 end subroutine update_faces_average
6843 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
6848 integer,
intent(in) :: ixi^
l, ixo^
l
6849 double precision,
intent(in) :: qt, qdt
6851 double precision,
intent(in) :: wp(ixi^s,1:nw)
6852 type(state) :: sct, s
6853 type(ct_velocity) :: vcts
6854 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6855 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6857 double precision :: circ(ixi^s,1:
ndim)
6859 double precision :: ecc(ixi^s,
sdim:3)
6860 double precision :: ein(ixi^s,
sdim:3)
6862 double precision :: el(ixi^s),er(ixi^s)
6864 double precision :: elc,erc
6866 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6868 double precision :: jce(ixi^s,
sdim:3)
6870 double precision :: xs(ixgs^t,1:
ndim)
6871 double precision :: gradi(ixgs^t)
6872 integer :: ixc^
l,ixa^
l
6873 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
6875 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
6878 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6884 {
do ix^db=iximin^db,iximax^db\}
6887 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_)
6888 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_)
6889 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_)
6892 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
6899 {
do ix^db=iximin^db,iximax^db\}
6902 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
6903 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
6904 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6907 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
6921 i1kr^d=kr(idim1,^d);
6924 i2kr^d=kr(idim2,^d);
6927 if (lvc(idim1,idim2,idir)==1)
then
6929 ixcmin^d=ixomin^d+kr(idir,^d)-1;
6932 {
do ix^db=ixcmin^db,ixcmax^db\}
6933 fe(ix^d,idir)=quarter*&
6934 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
6935 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
6940 ixamax^d=ixcmax^d+i1kr^d;
6941 {
do ix^db=ixamin^db,ixamax^db\}
6942 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
6943 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
6946 do ix^db=ixcmin^db,ixcmax^db\}
6947 if(vnorm(ix^d,idim1)>0.d0)
then
6949 else if(vnorm(ix^d,idim1)<0.d0)
then
6950 elc=el({ix^d+i1kr^d})
6952 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
6954 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
6956 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
6957 erc=er({ix^d+i1kr^d})
6959 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
6961 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
6966 ixamax^d=ixcmax^d+i2kr^d;
6967 {
do ix^db=ixamin^db,ixamax^db\}
6968 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
6969 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
6972 do ix^db=ixcmin^db,ixcmax^db\}
6973 if(vnorm(ix^d,idim2)>0.d0)
then
6975 else if(vnorm(ix^d,idim2)<0.d0)
then
6976 elc=el({ix^d+i2kr^d})
6978 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
6980 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
6982 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
6983 erc=er({ix^d+i2kr^d})
6985 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
6987 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
6991 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
6996 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7010 if (lvc(idim1,idim2,idir)==0) cycle
7012 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7013 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7016 xs(ixa^s,:)=x(ixa^s,:)
7017 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7018 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi,.false.)
7019 if (lvc(idim1,idim2,idir)==1)
then
7020 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7022 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7029 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7031 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7035 {
do ix^db=ixomin^db,ixomax^db\}
7036 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7037 +ein(ix1,ix2-1,ix3-1,idir))
7038 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7039 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7041 else if(idir==2)
then
7042 {
do ix^db=ixomin^db,ixomax^db\}
7043 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7044 +ein(ix1-1,ix2,ix3-1,idir))
7045 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7046 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7049 {
do ix^db=ixomin^db,ixomax^db\}
7050 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7051 +ein(ix1-1,ix2-1,ix3,idir))
7052 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7053 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7059 {
do ix^db=ixomin^db,ixomax^db\}
7060 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7061 +ein(ix1-1,ix2-1,idir))
7062 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7063 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7068 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7074 if(
associated(usr_set_electric_field)) &
7075 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7077 circ(ixi^s,1:ndim)=zero
7082 ixcmin^d=ixomin^d-kr(idim1,^d);
7084 ixa^l=ixc^l-kr(idim2,^d);
7087 if(lvc(idim1,idim2,idir)==1)
then
7089 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7092 else if(lvc(idim1,idim2,idir)==-1)
then
7094 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7101 where(s%surfaceC(ixc^s,idim1) > smalldouble)
7102 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
7104 circ(ixc^s,idim1)=zero
7107 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
7112 end subroutine update_faces_contact
7115 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
7120 integer,
intent(in) :: ixi^
l, ixo^
l
7121 double precision,
intent(in) :: qt, qdt
7122 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7123 type(state) :: sct, s
7124 type(ct_velocity) :: vcts
7126 double precision :: vtill(ixi^s,2)
7127 double precision :: vtilr(ixi^s,2)
7128 double precision :: bfacetot(ixi^s,
ndim)
7129 double precision :: btill(ixi^s,
ndim)
7130 double precision :: btilr(ixi^s,
ndim)
7131 double precision :: cp(ixi^s,2)
7132 double precision :: cm(ixi^s,2)
7133 double precision :: circ(ixi^s,1:
ndim)
7135 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7136 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7137 integer :: idim1,idim2,idir
7139 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7140 cbarmax=>vcts%cbarmax)
7153 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
7169 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7173 idim2=mod(idir+1,3)+1
7175 jxc^
l=ixc^
l+
kr(idim1,^
d);
7176 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7180 vtill(ixi^s,2),vtilr(ixi^s,2))
7183 vtill(ixi^s,1),vtilr(ixi^s,1))
7189 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7190 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7192 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7193 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7196 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7199 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7203 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7204 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7206 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7207 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7211 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7212 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7213 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7214 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7215 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7216 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7217 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7218 /(cp(ixc^s,2)+cm(ixc^s,2))
7221 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7225 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7239 circ(ixi^s,1:
ndim)=zero
7244 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7248 if(
lvc(idim1,idim2,idir)/=0)
then
7249 hxc^
l=ixc^
l-
kr(idim2,^
d);
7251 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7252 +
lvc(idim1,idim2,idir)&
7259 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
7260 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
7262 circ(ixc^s,idim1)=zero
7265 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
7269 end subroutine update_faces_hll
7272 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
7277 integer,
intent(in) :: ixi^
l, ixo^
l
7278 type(state),
intent(in) :: sct, s
7280 double precision :: jce(ixi^s,
sdim:3)
7283 double precision :: jcc(ixi^s,7-2*
ndir:3)
7285 double precision :: xs(ixgs^t,1:
ndim)
7287 double precision :: eta(ixi^s)
7288 double precision :: gradi(ixgs^t)
7289 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7291 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7297 if (
lvc(idim1,idim2,idir)==0) cycle
7299 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7300 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7303 xs(ixb^s,:)=x(ixb^s,:)
7304 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7305 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,.true.)
7306 if (
lvc(idim1,idim2,idir)==1)
then
7307 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7309 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7316 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7324 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7325 jcc(ixc^s,idir)=0.d0
7327 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7328 ixamin^
d=ixcmin^
d+ix^
d;
7329 ixamax^
d=ixcmax^
d+ix^
d;
7330 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7332 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7333 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7338 end subroutine get_resistive_electric_field
7341 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7344 integer,
intent(in) :: ixi^
l, ixo^
l
7345 double precision,
intent(in) :: w(ixi^s,1:nw)
7346 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7347 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7349 double precision :: jxbxb(ixi^s,1:3)
7350 integer :: idir,ixa^
l,ixc^
l,ix^
d
7353 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7360 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7363 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7364 ixamin^
d=ixcmin^
d+ix^
d;
7365 ixamax^
d=ixcmax^
d+ix^
d;
7366 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7368 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7371 end subroutine get_ambipolar_electric_field
7377 integer,
intent(in) :: ixo^
l
7387 do ix^db=ixomin^db,ixomax^db\}
7389 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7390 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7391 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7392 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7393 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7394 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7397 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7398 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7399 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7400 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7443 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7444 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7445 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7447 double precision :: adummy(ixis^s,1:3)
7453 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7456 integer,
intent(in) :: ixi^
l, ixo^
l
7457 double precision,
intent(in) :: w(ixi^s,1:nw)
7458 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7459 double precision,
intent(out):: rfactor(ixi^s)
7461 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7465 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)
7467 end subroutine rfactor_from_temperature_ionization
7469 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7471 integer,
intent(in) :: ixi^
l, ixo^
l
7472 double precision,
intent(in) :: w(ixi^s,1:nw)
7473 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7474 double precision,
intent(out):: rfactor(ixi^s)
7478 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.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
subroutine, public get_divb(w, ixil, ixol, divb, fourthorder)
Calculate div B within ixO.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
subroutine gradient(q, ixil, ixol, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
integer, parameter cylindrical
subroutine gradients(q, ixil, ixol, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
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 gradientx(q, x, ixil, ixol, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
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)
logical, public, protected mhd_divb_4thorder
Whether divB is computed with a fourth order approximation.
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 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.