22 double precision,
public ::
mhd_eta = 0.0d0
30 double precision,
protected :: small_e
41 double precision :: divbdiff = 0.8d0
46 double precision,
public,
protected ::
h_ion_fr=1d0
49 double precision,
public,
protected ::
he_ion_fr=1d0
56 double precision,
public,
protected ::
rr=1d0
58 double precision :: gamma_1, inv_gamma_1
60 double precision :: inv_squared_c0, inv_squared_c
67 integer,
public,
protected ::
rho_
69 integer,
allocatable,
public,
protected ::
mom(:)
71 integer,
public,
protected :: ^
c&m^C_
73 integer,
public,
protected ::
e_
75 integer,
public,
protected :: ^
c&b^C_
77 integer,
public,
protected ::
p_
79 integer,
public,
protected ::
q_
81 integer,
public,
protected ::
psi_
83 integer,
public,
protected ::
te_
88 integer,
allocatable,
public,
protected ::
tracer(:)
96 integer,
parameter :: divb_none = 0
97 integer,
parameter :: divb_multigrid = -1
98 integer,
parameter :: divb_glm = 1
99 integer,
parameter :: divb_powel = 2
100 integer,
parameter :: divb_janhunen = 3
101 integer,
parameter :: divb_linde = 4
102 integer,
parameter :: divb_lindejanhunen = 5
103 integer,
parameter :: divb_lindepowel = 6
104 integer,
parameter :: divb_lindeglm = 7
105 integer,
parameter :: divb_ct = 8
135 logical,
public,
protected ::
mhd_glm = .false.
170 logical :: compactres = .false.
184 logical :: total_energy = .true.
188 logical :: gravity_energy
190 logical :: gravity_rhov = .false.
192 character(len=std_len),
public,
protected ::
typedivbfix =
'linde'
194 character(len=std_len),
public,
protected ::
type_ct =
'uct_contact'
196 character(len=std_len) :: typedivbdiff =
'all'
207 subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
209 integer,
intent(in) :: ixi^
l, ixo^
l
210 double precision,
intent(in) :: x(ixi^s,1:
ndim)
211 double precision,
intent(in) :: w(ixi^s,1:nw)
212 double precision,
intent(inout) :: res(ixi^s)
213 end subroutine mask_subroutine
250 subroutine mhd_read_params(files)
253 character(len=*),
intent(in) :: files(:)
269 do n = 1,
size(files)
270 open(
unitpar, file=trim(files(n)), status=
"old")
271 read(
unitpar, mhd_list,
end=111)
275 end subroutine mhd_read_params
278 subroutine mhd_write_info(fh)
280 integer,
intent(in) :: fh
283 integer,
parameter :: n_par = 1
284 double precision :: values(n_par)
285 integer,
dimension(MPI_STATUS_SIZE) :: st
286 character(len=name_len) :: names(n_par)
288 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
292 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
293 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
294 end subroutine mhd_write_info
321 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
332 if(
mype==0)
write(*,*)
'WARNING: set mhd_internal_e=F when mhd_energy=F'
336 if(
mype==0)
write(*,*)
'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
340 if(
mype==0)
write(*,*)
'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
344 if(
mype==0)
write(*,*)
'WARNING: set mhd_hyperbolic_thermal_conduction=F when mhd_energy=F'
348 if(
mype==0)
write(*,*)
'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
352 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac=F when mhd_energy=F'
356 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
360 if(
mype==0)
write(*,*)
'WARNING: set B0field=F when mhd_energy=F'
364 if(
mype==0)
write(*,*)
'WARNING: set has_equi_rho0=F when mhd_energy=F'
368 if(
mype==0)
write(*,*)
'WARNING: set has_equi_pe0=F when mhd_energy=F'
374 if(
mype==0)
write(*,*)
'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
380 if(
mype==0)
write(*,*)
'WARNING: turn off parabolic TC when using hyperbolic TC'
405 phys_total_energy=total_energy
408 gravity_energy=.false.
410 gravity_energy=.true.
419 gravity_energy=.false.
425 if(
mype==0)
write(*,*)
'WARNING: reset mhd_trac_type=1 for 1D simulation'
430 if(
mype==0)
write(*,*)
'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
439 type_divb = divb_none
442 type_divb = divb_multigrid
444 mg%operator_type = mg_laplacian
451 case (
'powel',
'powell')
452 type_divb = divb_powel
454 type_divb = divb_janhunen
456 type_divb = divb_linde
457 case (
'lindejanhunen')
458 type_divb = divb_lindejanhunen
460 type_divb = divb_lindepowel
464 type_divb = divb_lindeglm
469 call mpistop(
'Unknown divB fix')
472 allocate(start_indices(number_species),stop_indices(number_species))
479 mom(:) = var_set_momentum(
ndir)
485 e_ = var_set_energy()
494 mag(:) = var_set_bfield(
ndir)
498 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
514 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
519 te_ = var_set_auxvar(
'Te',
'Te')
528 stop_indices(1)=nwflux
559 allocate(iw_vector(nvector))
560 iw_vector(1) =
mom(1) - 1
561 iw_vector(2) = mag(1) - 1
564 if (.not.
allocated(flux_type))
then
565 allocate(flux_type(
ndir, nwflux))
566 flux_type = flux_default
567 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
568 call mpistop(
"phys_check error: flux_type has wrong shape")
571 if(nwflux>mag(
ndir))
then
573 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
578 flux_type(:,
psi_)=flux_special
580 flux_type(idir,mag(idir))=flux_special
584 flux_type(idir,mag(idir))=flux_tvdlf
590 phys_get_dt => mhd_get_dt
593 phys_get_cmax => mhd_get_cmax_semirelati
595 phys_get_cmax => mhd_get_cmax_semirelati_noe
599 phys_get_cmax => mhd_get_cmax_origin
601 phys_get_cmax => mhd_get_cmax_origin_noe
604 phys_get_a2max => mhd_get_a2max
605 phys_get_tcutoff => mhd_get_tcutoff
606 phys_get_h_speed => mhd_get_h_speed
608 phys_get_cbounds => mhd_get_cbounds_split_rho
610 phys_get_cbounds => mhd_get_cbounds_semirelati
612 phys_get_cbounds => mhd_get_cbounds
615 phys_to_primitive => mhd_to_primitive_hde
617 phys_to_conserved => mhd_to_conserved_hde
621 phys_to_primitive => mhd_to_primitive_semirelati
623 phys_to_conserved => mhd_to_conserved_semirelati
626 phys_to_primitive => mhd_to_primitive_semirelati_noe
628 phys_to_conserved => mhd_to_conserved_semirelati_noe
633 phys_to_primitive => mhd_to_primitive_split_rho
635 phys_to_conserved => mhd_to_conserved_split_rho
638 phys_to_primitive => mhd_to_primitive_inte
640 phys_to_conserved => mhd_to_conserved_inte
643 phys_to_primitive => mhd_to_primitive_origin
645 phys_to_conserved => mhd_to_conserved_origin
648 phys_to_primitive => mhd_to_primitive_origin_noe
650 phys_to_conserved => mhd_to_conserved_origin_noe
655 phys_get_flux => mhd_get_flux_hde
658 phys_get_flux => mhd_get_flux_semirelati
660 phys_get_flux => mhd_get_flux_semirelati_noe
664 phys_get_flux => mhd_get_flux_split
666 phys_get_flux => mhd_get_flux
668 phys_get_flux => mhd_get_flux_noe
673 phys_add_source_geom => mhd_add_source_geom_semirelati
675 phys_add_source_geom => mhd_add_source_geom_split
677 phys_add_source_geom => mhd_add_source_geom
679 phys_add_source => mhd_add_source
680 phys_check_params => mhd_check_params
681 phys_write_info => mhd_write_info
684 phys_handle_small_values => mhd_handle_small_values_inte
685 mhd_handle_small_values => mhd_handle_small_values_inte
686 phys_check_w => mhd_check_w_inte
688 phys_handle_small_values => mhd_handle_small_values_hde
689 mhd_handle_small_values => mhd_handle_small_values_hde
690 phys_check_w => mhd_check_w_hde
692 phys_handle_small_values => mhd_handle_small_values_semirelati
693 mhd_handle_small_values => mhd_handle_small_values_semirelati
694 phys_check_w => mhd_check_w_semirelati
696 phys_handle_small_values => mhd_handle_small_values_split
697 mhd_handle_small_values => mhd_handle_small_values_split
698 phys_check_w => mhd_check_w_split
700 phys_handle_small_values => mhd_handle_small_values_origin
701 mhd_handle_small_values => mhd_handle_small_values_origin
702 phys_check_w => mhd_check_w_origin
704 phys_handle_small_values => mhd_handle_small_values_noe
705 mhd_handle_small_values => mhd_handle_small_values_noe
706 phys_check_w => mhd_check_w_noe
710 phys_get_pthermal => mhd_get_pthermal_inte
713 phys_get_pthermal => mhd_get_pthermal_hde
716 phys_get_pthermal => mhd_get_pthermal_semirelati
719 phys_get_pthermal => mhd_get_pthermal_origin
722 phys_get_pthermal => mhd_get_pthermal_noe
727 phys_set_equi_vars => set_equi_vars_grid
730 if(type_divb==divb_glm)
then
731 phys_modify_wlr => mhd_modify_wlr
736 mhd_get_rfactor=>rfactor_from_temperature_ionization
737 phys_update_temperature => mhd_update_temperature
741 mhd_get_rfactor=>rfactor_from_constant_ionization
764 phys_get_ct_velocity => mhd_get_ct_velocity
765 phys_update_faces => mhd_update_faces
767 phys_modify_wlr => mhd_modify_wlr
769 phys_boundary_adjust => mhd_boundary_adjust
778 call mhd_physical_units()
784 call mpistop(
"thermal conduction needs mhd_energy=T")
787 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
790 call mpistop(
"radiative cooling needs mhd_energy=T")
801 if(phys_internal_e)
then
803 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
805 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
809 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot_with_equi
811 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
815 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
817 tc_fl%has_equi = .true.
818 tc_fl%get_temperature_equi => mhd_get_temperature_equi
819 tc_fl%get_rho_equi => mhd_get_rho_equi
821 tc_fl%has_equi = .false.
824 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
848 rc_fl%get_var_Rfactor => mhd_get_rfactor
852 rc_fl%has_equi = .true.
853 rc_fl%get_rho_equi => mhd_get_rho_equi
854 rc_fl%get_pthermal_equi => mhd_get_pe_equi
856 rc_fl%has_equi = .false.
862 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
864 phys_te_images => mhd_te_images
880 if (particles_eta < zero) particles_eta =
mhd_eta
881 if (particles_etah < zero) particles_eta =
mhd_etah
883 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
884 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
898 phys_wider_stencil = 2
900 phys_wider_stencil = 1
908 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
920 phys_wider_stencil = 2
922 phys_wider_stencil = 1
936 subroutine mhd_te_images
941 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
943 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
945 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
947 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
950 call mpistop(
"Error in synthesize emission: Unknown convert_type")
952 end subroutine mhd_te_images
958 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
962 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
963 double precision,
intent(in) :: x(ixi^s,1:
ndim)
964 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
965 double precision,
intent(in) :: my_dt
966 logical,
intent(in) :: fix_conserve_at_step
968 end subroutine mhd_sts_set_source_tc_mhd
970 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
977 integer,
intent(in) :: ixi^
l, ixo^
l
978 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
979 double precision,
intent(in) :: w(ixi^s,1:nw)
980 double precision :: dtnew
983 end function mhd_get_tc_dt_mhd
985 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
988 integer,
intent(in) :: ixi^
l,ixo^
l
989 double precision,
intent(inout) :: w(ixi^s,1:nw)
990 double precision,
intent(in) :: x(ixi^s,1:
ndim)
991 integer,
intent(in) :: step
992 character(len=140) :: error_msg
994 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
995 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
996 end subroutine mhd_tc_handle_small_e
999 subroutine tc_params_read_mhd(fl)
1001 type(tc_fluid),
intent(inout) :: fl
1003 double precision :: tc_k_para=0d0
1004 double precision :: tc_k_perp=0d0
1007 logical :: tc_perpendicular=.false.
1008 logical :: tc_saturate=.false.
1009 character(len=std_len) :: tc_slope_limiter=
"MC"
1011 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1015 read(
unitpar, tc_list,
end=111)
1019 fl%tc_perpendicular = tc_perpendicular
1020 fl%tc_saturate = tc_saturate
1021 fl%tc_k_para = tc_k_para
1022 fl%tc_k_perp = tc_k_perp
1023 select case(tc_slope_limiter)
1025 fl%tc_slope_limiter = 0
1028 fl%tc_slope_limiter = 1
1031 fl%tc_slope_limiter = 2
1034 fl%tc_slope_limiter = 3
1037 fl%tc_slope_limiter = 4
1039 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1041 end subroutine tc_params_read_mhd
1045 subroutine rc_params_read(fl)
1048 type(rc_fluid),
intent(inout) :: fl
1050 double precision :: cfrac=0.1d0
1053 double precision :: rad_cut_hgt=0.5d0
1054 double precision :: rad_cut_dey=0.15d0
1057 integer :: ncool = 4000
1059 logical :: tfix=.false.
1061 logical :: rc_split=.false.
1062 logical :: rad_cut=.false.
1064 character(len=std_len) :: coolcurve=
'JCcorona'
1066 character(len=std_len) :: coolmethod=
'exact'
1068 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1072 read(
unitpar, rc_list,
end=111)
1077 fl%coolcurve=coolcurve
1078 fl%coolmethod=coolmethod
1081 fl%rc_split=rc_split
1084 fl%rad_cut_hgt=rad_cut_hgt
1085 fl%rad_cut_dey=rad_cut_dey
1086 end subroutine rc_params_read
1090 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1093 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1094 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1096 double precision :: delx(ixi^s,1:
ndim)
1097 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1098 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1104 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1108 hxo^
l=ixo^
l-
kr(idims,^
d);
1114 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1117 xshift^
d=half*(one-
kr(^
d,idims));
1124 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1128 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1131 end subroutine set_equi_vars_grid_faces
1134 subroutine set_equi_vars_grid(igrid)
1138 integer,
intent(in) :: igrid
1144 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1146 end subroutine set_equi_vars_grid
1149 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1151 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1152 double precision,
intent(in) :: w(ixi^s, 1:nw)
1153 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1154 double precision :: wnew(ixo^s, 1:nwc)
1161 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1167 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1171 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1175 if(
b0field .and. total_energy)
then
1176 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1177 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1181 end function convert_vars_splitting
1183 subroutine mhd_check_params
1191 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1192 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1196 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1197 inv_gamma_1=1.d0/gamma_1
1202 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1207 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1212 end subroutine mhd_check_params
1214 subroutine mhd_physical_units()
1216 double precision :: mp,kb,miu0,c_lightspeed
1217 double precision :: a,
b
1228 c_lightspeed=const_c
1369 end subroutine mhd_physical_units
1371 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1374 logical,
intent(in) :: primitive
1375 logical,
intent(inout) :: flag(ixi^s,1:nw)
1376 integer,
intent(in) :: ixi^
l, ixo^
l
1377 double precision,
intent(in) :: w(ixi^s,nw)
1379 double precision :: tmp,b2,
b(ixo^s,1:
ndir)
1380 double precision :: v(ixo^s,1:
ndir),gamma2,inv_rho
1391 {
do ix^db=ixomin^db,ixomax^db \}
1392 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1395 {
do ix^db=ixomin^db,ixomax^db \}
1396 b2=(^
c&w(ix^d,
b^
c_)**2+)
1397 if(b2>smalldouble)
then
1402 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1403 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1404 inv_rho = 1d0/w(ix^d,
rho_)
1406 b2=b2*inv_rho*inv_squared_c
1408 gamma2=1.d0/(1.d0+b2)
1410 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1413 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1414 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1415 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1420 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1426 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1427 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1428 if(tmp<small_e) flag(ix^d,
e_)=.true.
1434 end subroutine mhd_check_w_semirelati
1436 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1439 logical,
intent(in) :: primitive
1440 integer,
intent(in) :: ixi^
l, ixo^
l
1441 double precision,
intent(in) :: w(ixi^s,nw)
1442 logical,
intent(inout) :: flag(ixi^s,1:nw)
1447 {
do ix^db=ixomin^db,ixomax^db\}
1453 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1457 end subroutine mhd_check_w_origin
1459 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1462 logical,
intent(in) :: primitive
1463 integer,
intent(in) :: ixi^
l, ixo^
l
1464 double precision,
intent(in) :: w(ixi^s,nw)
1465 logical,
intent(inout) :: flag(ixi^s,1:nw)
1467 double precision :: tmp
1471 {
do ix^db=ixomin^db,ixomax^db\}
1477 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1482 end subroutine mhd_check_w_split
1484 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1487 logical,
intent(in) :: primitive
1488 integer,
intent(in) :: ixi^
l, ixo^
l
1489 double precision,
intent(in) :: w(ixi^s,nw)
1490 logical,
intent(inout) :: flag(ixi^s,1:nw)
1495 {
do ix^db=ixomin^db,ixomax^db\}
1499 end subroutine mhd_check_w_noe
1501 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1504 logical,
intent(in) :: primitive
1505 integer,
intent(in) :: ixi^
l, ixo^
l
1506 double precision,
intent(in) :: w(ixi^s,nw)
1507 logical,
intent(inout) :: flag(ixi^s,1:nw)
1512 {
do ix^db=ixomin^db,ixomax^db\}
1517 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1521 end subroutine mhd_check_w_inte
1523 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1526 logical,
intent(in) :: primitive
1527 integer,
intent(in) :: ixi^
l, ixo^
l
1528 double precision,
intent(in) :: w(ixi^s,nw)
1529 logical,
intent(inout) :: flag(ixi^s,1:nw)
1534 {
do ix^db=ixomin^db,ixomax^db\}
1539 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1543 end subroutine mhd_check_w_hde
1546 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1548 integer,
intent(in) :: ixi^
l, ixo^
l
1549 double precision,
intent(inout) :: w(ixi^s, nw)
1550 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1554 {
do ix^db=ixomin^db,ixomax^db\}
1556 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1558 +(^
c&w(ix^
d,
b^
c_)**2+))
1563 end subroutine mhd_to_conserved_origin
1566 subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
1568 integer,
intent(in) :: ixi^
l, ixo^
l
1569 double precision,
intent(inout) :: w(ixi^s, nw)
1570 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1574 {
do ix^db=ixomin^db,ixomax^db\}
1579 end subroutine mhd_to_conserved_origin_noe
1582 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1584 integer,
intent(in) :: ixi^
l, ixo^
l
1585 double precision,
intent(inout) :: w(ixi^s, nw)
1586 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1590 {
do ix^db=ixomin^db,ixomax^db\}
1592 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1598 end subroutine mhd_to_conserved_hde
1601 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1603 integer,
intent(in) :: ixi^
l, ixo^
l
1604 double precision,
intent(inout) :: w(ixi^s, nw)
1605 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1609 {
do ix^db=ixomin^db,ixomax^db\}
1611 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1616 end subroutine mhd_to_conserved_inte
1619 subroutine mhd_to_conserved_split_rho(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 :: rho
1628 {
do ix^db=ixomin^db,ixomax^db\}
1631 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1632 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1633 +(^
c&w(ix^
d,
b^
c_)**2+))
1638 end subroutine mhd_to_conserved_split_rho
1641 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1643 integer,
intent(in) :: ixi^
l, ixo^
l
1644 double precision,
intent(inout) :: w(ixi^s, nw)
1645 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1647 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1650 {
do ix^db=ixomin^db,ixomax^db\}
1652 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1653 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1654 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1655 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1656 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1657 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1662 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1663 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1664 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1672 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1676 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1678 +(^
c&w(ix^
d,
b^
c_)**2+)&
1679 +(^
c&e(ix^
d,^
c)**2+)*inv_squared_c)
1687 end subroutine mhd_to_conserved_semirelati
1689 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1691 integer,
intent(in) :: ixi^
l, ixo^
l
1692 double precision,
intent(inout) :: w(ixi^s, nw)
1693 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1695 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1698 {
do ix^db=ixomin^db,ixomax^db\}
1700 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1701 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1702 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1703 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1704 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1705 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1710 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1711 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1712 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1722 end subroutine mhd_to_conserved_semirelati_noe
1725 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1727 integer,
intent(in) :: ixi^
l, ixo^
l
1728 double precision,
intent(inout) :: w(ixi^s, nw)
1729 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1731 double precision :: inv_rho
1736 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1739 {
do ix^db=ixomin^db,ixomax^db\}
1740 inv_rho = 1.d0/w(ix^
d,
rho_)
1744 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1746 +(^
c&w(ix^
d,
b^
c_)**2+)))
1749 end subroutine mhd_to_primitive_origin
1752 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1754 integer,
intent(in) :: ixi^
l, ixo^
l
1755 double precision,
intent(inout) :: w(ixi^s, nw)
1756 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1758 double precision :: inv_rho
1763 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1766 {
do ix^db=ixomin^db,ixomax^db\}
1767 inv_rho = 1.d0/w(ix^
d,
rho_)
1772 end subroutine mhd_to_primitive_origin_noe
1775 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1777 integer,
intent(in) :: ixi^
l, ixo^
l
1778 double precision,
intent(inout) :: w(ixi^s, nw)
1779 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1781 double precision :: inv_rho
1786 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1789 {
do ix^db=ixomin^db,ixomax^db\}
1790 inv_rho = 1d0/w(ix^
d,
rho_)
1797 end subroutine mhd_to_primitive_hde
1800 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1802 integer,
intent(in) :: ixi^
l, ixo^
l
1803 double precision,
intent(inout) :: w(ixi^s, nw)
1804 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1806 double precision :: inv_rho
1811 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1814 {
do ix^db=ixomin^db,ixomax^db\}
1816 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1818 inv_rho = 1.d0/w(ix^
d,
rho_)
1822 end subroutine mhd_to_primitive_inte
1825 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1827 integer,
intent(in) :: ixi^
l, ixo^
l
1828 double precision,
intent(inout) :: w(ixi^s, nw)
1829 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1831 double precision :: inv_rho
1836 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1839 {
do ix^db=ixomin^db,ixomax^db\}
1844 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1846 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1849 end subroutine mhd_to_primitive_split_rho
1852 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1854 integer,
intent(in) :: ixi^
l, ixo^
l
1855 double precision,
intent(inout) :: w(ixi^s, nw)
1856 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1858 double precision ::
b(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
1863 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1866 {
do ix^db=ixomin^db,ixomax^db\}
1867 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1868 if(b2>smalldouble)
then
1876 inv_rho=1.d0/w(ix^
d,
rho_)
1878 b2=b2*inv_rho*inv_squared_c
1880 gamma2=1.d0/(1.d0+b2)
1882 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1886 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1890 b(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1891 b(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1892 b(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1896 b(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1902 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1904 +(^
c&w(ix^
d,
b^
c_)**2+)&
1905 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
1909 end subroutine mhd_to_primitive_semirelati
1912 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1914 integer,
intent(in) :: ixi^
l, ixo^
l
1915 double precision,
intent(inout) :: w(ixi^s, nw)
1916 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1918 double precision ::
b(ixo^s,1:
ndir),tmp,b2,gamma2,inv_rho
1919 integer :: ix^
d, idir
1923 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1926 {
do ix^db=ixomin^db,ixomax^db\}
1927 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1928 if(b2>smalldouble)
then
1936 inv_rho=1.d0/w(ix^
d,
rho_)
1938 b2=b2*inv_rho*inv_squared_c
1940 gamma2=1.d0/(1.d0+b2)
1942 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1945 end subroutine mhd_to_primitive_semirelati_noe
1950 integer,
intent(in) :: ixi^
l, ixo^
l
1951 double precision,
intent(inout) :: w(ixi^s, nw)
1952 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1957 {
do ix^db=ixomin^db,ixomax^db\}
1960 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1962 +(^
c&w(ix^
d,
b^
c_)**2+))
1965 {
do ix^db=ixomin^db,ixomax^db\}
1967 w(ix^d,
e_)=w(ix^d,
e_)&
1968 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1969 +(^
c&w(ix^d,
b^
c_)**2+))
1976 subroutine mhd_ei_to_e_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)
1984 {
do ix^db=ixomin^db,ixomax^db\}
1990 end subroutine mhd_ei_to_e_hde
1993 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
1995 integer,
intent(in) :: ixi^
l, ixo^
l
1996 double precision,
intent(inout) :: w(ixi^s, nw)
1997 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1999 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2000 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2002 end subroutine mhd_ei_to_e_semirelati
2007 integer,
intent(in) :: ixi^
l, ixo^
l
2008 double precision,
intent(inout) :: w(ixi^s, nw)
2009 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2014 {
do ix^db=ixomin^db,ixomax^db\}
2017 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2019 +(^
c&w(ix^
d,
b^
c_)**2+))
2022 {
do ix^db=ixomin^db,ixomax^db\}
2024 w(ix^d,
e_)=w(ix^d,
e_)&
2025 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2026 +(^
c&w(ix^d,
b^
c_)**2+))
2030 if(fix_small_values)
then
2031 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2037 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2039 integer,
intent(in) :: ixi^
l, ixo^
l
2040 double precision,
intent(inout) :: w(ixi^s, nw)
2041 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2045 {
do ix^db=ixomin^db,ixomax^db\}
2051 if(fix_small_values)
then
2052 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2055 end subroutine mhd_e_to_ei_hde
2058 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2060 integer,
intent(in) :: ixi^
l, ixo^
l
2061 double precision,
intent(inout) :: w(ixi^s, nw)
2062 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2064 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2065 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2067 end subroutine mhd_e_to_ei_semirelati
2069 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2072 logical,
intent(in) :: primitive
2073 integer,
intent(in) :: ixi^
l,ixo^
l
2074 double precision,
intent(inout) :: w(ixi^s,1:nw)
2075 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2076 character(len=*),
intent(in) :: subname
2078 double precision ::
b(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2079 double precision :: tmp, b2, gamma2, inv_rho
2081 logical :: flag(ixi^s,1:nw)
2090 {
do ix^db=ixomin^db,ixomax^db\}
2091 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2092 if(b2>smalldouble)
then
2099 inv_rho=1.d0/w(ix^
d,
rho_)
2101 b2=b2*inv_rho*inv_squared_c
2103 gamma2=1.d0/(1.d0+b2)
2105 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
2108 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2109 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2110 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2114 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2120 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2121 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2122 +(^
c&w(ix^
d,
b^
c_)**2+)&
2123 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
2130 select case (small_values_method)
2132 {
do ix^db=ixomin^db,ixomax^db\}
2133 if(flag(ix^d,
rho_))
then
2134 w(ix^d,
rho_) = small_density
2135 ^
c&w(ix^d,
m^
c_)=0.d0\
2139 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2141 if(flag(ix^d,
e_))
then
2142 w(ix^d,
e_)=small_pressure*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2143 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2150 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2153 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2155 w(ixo^s,
e_)=pressure(ixo^s)
2156 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2157 {
do ix^db=ixomin^db,ixomax^db\}
2158 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2159 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2164 if(.not.primitive)
then
2166 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2167 w(ixo^s,
e_)=pressure(ixo^s)
2169 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2173 end subroutine mhd_handle_small_values_semirelati
2175 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2178 logical,
intent(in) :: primitive
2179 integer,
intent(in) :: ixi^
l,ixo^
l
2180 double precision,
intent(inout) :: w(ixi^s,1:nw)
2181 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2182 character(len=*),
intent(in) :: subname
2185 logical :: flag(ixi^s,1:nw)
2187 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2192 {
do ix^db=ixomin^db,ixomax^db\}
2196 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2208 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2210 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2213 {
do ix^db=iximin^db,iximax^db\}
2214 w(ix^d,
e_)=w(ix^d,
e_)&
2215 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2217 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2219 {
do ix^db=iximin^db,iximax^db\}
2220 w(ix^d,
e_)=w(ix^d,
e_)&
2221 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2225 if(.not.primitive)
then
2227 {
do ix^db=ixomin^db,ixomax^db\}
2229 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2230 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2233 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2237 end subroutine mhd_handle_small_values_origin
2239 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2242 logical,
intent(in) :: primitive
2243 integer,
intent(in) :: ixi^
l,ixo^
l
2244 double precision,
intent(inout) :: w(ixi^s,1:nw)
2245 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2246 character(len=*),
intent(in) :: subname
2248 double precision :: rho
2250 logical :: flag(ixi^s,1:nw)
2252 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2257 {
do ix^db=ixomin^db,ixomax^db\}
2262 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2269 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2275 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2277 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2280 {
do ix^db=iximin^db,iximax^db\}
2282 w(ix^d,
e_)=w(ix^d,
e_)&
2283 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2285 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2287 {
do ix^db=iximin^db,iximax^db\}
2289 w(ix^d,
e_)=w(ix^d,
e_)&
2290 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2294 if(.not.primitive)
then
2296 {
do ix^db=ixomin^db,ixomax^db\}
2298 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2299 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2300 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2303 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2307 end subroutine mhd_handle_small_values_split
2309 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2312 logical,
intent(in) :: primitive
2313 integer,
intent(in) :: ixi^
l,ixo^
l
2314 double precision,
intent(inout) :: w(ixi^s,1:nw)
2315 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2316 character(len=*),
intent(in) :: subname
2319 logical :: flag(ixi^s,1:nw)
2321 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2326 {
do ix^db=ixomin^db,ixomax^db\}
2327 if(flag(ix^
d,
rho_))
then
2329 ^
c&w(ix^
d,
m^
c_)=0.d0\
2334 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2339 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2341 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2343 if(.not.primitive)
then
2345 {
do ix^db=ixomin^db,ixomax^db\}
2347 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2350 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2354 end subroutine mhd_handle_small_values_inte
2356 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2359 logical,
intent(in) :: primitive
2360 integer,
intent(in) :: ixi^
l,ixo^
l
2361 double precision,
intent(inout) :: w(ixi^s,1:nw)
2362 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2363 character(len=*),
intent(in) :: subname
2366 logical :: flag(ixi^s,1:nw)
2368 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2373 {
do ix^db=ixomin^db,ixomax^db\}
2377 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2383 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2385 if(.not.primitive)
then
2387 {
do ix^db=ixomin^db,ixomax^db\}
2391 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2395 end subroutine mhd_handle_small_values_noe
2397 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2400 logical,
intent(in) :: primitive
2401 integer,
intent(in) :: ixi^
l,ixo^
l
2402 double precision,
intent(inout) :: w(ixi^s,1:nw)
2403 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2404 character(len=*),
intent(in) :: subname
2407 logical :: flag(ixi^s,1:nw)
2409 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2414 {
do ix^db=ixomin^db,ixomax^db\}
2415 if(flag(ix^
d,
rho_))
then
2417 ^
c&w(ix^
d,
m^
c_)=0.d0\
2422 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2427 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2429 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2431 if(.not.primitive)
then
2433 {
do ix^db=ixomin^db,ixomax^db\}
2435 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2438 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2442 end subroutine mhd_handle_small_values_hde
2448 integer,
intent(in) :: ixi^
l, ixo^
l
2449 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2450 double precision,
intent(out) :: v(ixi^s,
ndir)
2452 double precision :: rho(ixi^s)
2457 rho(ixo^s)=1.d0/rho(ixo^s)
2460 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2466 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2469 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2470 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2471 double precision,
intent(inout) :: cmax(ixi^s)
2473 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2479 {
do ix^db=ixomin^db,ixomax^db \}
2490 cfast2=b2*inv_rho+cmax(ix^
d)
2491 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2492 if(avmincs2<zero) avmincs2=zero
2493 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2497 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2499 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2502 {
do ix^db=ixomin^db,ixomax^db \}
2512 b2=(^
c&w(ix^d,
b^
c_)**2+)
2513 cfast2=b2*inv_rho+cmax(ix^d)
2514 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2515 if(avmincs2<zero) avmincs2=zero
2516 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2520 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2522 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2526 end subroutine mhd_get_cmax_origin
2529 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2532 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2533 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2534 double precision,
intent(inout) :: cmax(ixi^s)
2536 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2542 {
do ix^db=ixomin^db,ixomax^db \}
2553 cfast2=b2*inv_rho+cmax(ix^
d)
2554 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2555 if(avmincs2<zero) avmincs2=zero
2556 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2560 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2562 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2565 {
do ix^db=ixomin^db,ixomax^db \}
2575 b2=(^
c&w(ix^d,
b^
c_)**2+)
2576 cfast2=b2*inv_rho+cmax(ix^d)
2577 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2578 if(avmincs2<zero) avmincs2=zero
2579 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2583 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2585 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2589 end subroutine mhd_get_cmax_origin_noe
2592 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2595 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2596 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2597 double precision,
intent(inout):: cmax(ixi^s)
2599 double precision :: csound, avmincs2, idim_alfven_speed2
2600 double precision :: inv_rho, alfven_speed2, gamma2
2603 {
do ix^db=ixomin^db,ixomax^db \}
2604 inv_rho=1.d0/w(ix^
d,
rho_)
2605 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2606 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2607 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2610 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2613 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2614 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2615 if(avmincs2<zero) avmincs2=zero
2617 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2618 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2621 end subroutine mhd_get_cmax_semirelati
2624 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2627 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2628 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2629 double precision,
intent(inout):: cmax(ixi^s)
2631 double precision :: csound, avmincs2, idim_alfven_speed2
2632 double precision :: inv_rho, alfven_speed2, gamma2
2635 {
do ix^db=ixomin^db,ixomax^db \}
2636 inv_rho=1.d0/w(ix^
d,
rho_)
2637 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2638 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2639 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2641 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2644 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2645 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2646 if(avmincs2<zero) avmincs2=zero
2648 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2649 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2652 end subroutine mhd_get_cmax_semirelati_noe
2654 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2657 integer,
intent(in) :: ixi^
l, ixo^
l
2658 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2659 double precision,
intent(inout) :: a2max(
ndim)
2660 double precision :: a2(ixi^s,
ndim,nw)
2661 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2666 hxo^
l=ixo^
l-
kr(i,^
d);
2667 gxo^
l=hxo^
l-
kr(i,^
d);
2668 jxo^
l=ixo^
l+
kr(i,^
d);
2669 kxo^
l=jxo^
l+
kr(i,^
d);
2670 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2671 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2672 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2674 end subroutine mhd_get_a2max
2677 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2680 integer,
intent(in) :: ixi^
l,ixo^
l
2681 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2683 double precision,
intent(inout) :: w(ixi^s,1:nw)
2684 double precision,
intent(out) :: tco_local,tmax_local
2686 double precision,
parameter :: trac_delta=0.25d0
2687 double precision :: te(ixi^s),lts(ixi^s)
2688 double precision,
dimension(1:ndim) :: bdir, bunitvec
2689 double precision,
dimension(ixI^S,1:ndim) :: gradt
2690 double precision :: ltrc,ltrp,altr
2691 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2692 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2695 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2697 call mhd_get_rfactor(w,x,ixi^
l,ixi^
l,te)
2698 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2701 tmax_local=maxval(te(ixo^s))
2709 do ix1=ixomin1,ixomax1
2710 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2711 if(lts(ix1)>trac_delta)
then
2712 tco_local=max(tco_local,te(ix1))
2724 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2725 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2726 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2727 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2729 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2740 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2747 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2752 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2753 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2757 if(bdir(1)/=0.d0)
then
2758 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2760 block%special_values(3)=0.d0
2762 if(bdir(2)/=0.d0)
then
2763 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2765 block%special_values(4)=0.d0
2769 if(bdir(1)/=0.d0)
then
2770 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2771 (bdir(3)/bdir(1))**2)
2773 block%special_values(3)=0.d0
2775 if(bdir(2)/=0.d0)
then
2776 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2777 (bdir(3)/bdir(2))**2)
2779 block%special_values(4)=0.d0
2781 if(bdir(3)/=0.d0)
then
2782 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2783 (bdir(2)/bdir(3))**2)
2785 block%special_values(5)=0.d0
2790 block%special_values(1)=zero
2791 {
do ix^db=ixomin^db,ixomax^db\}
2793 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2795 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2798 if(bdir(1)/=0.d0)
then
2799 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2803 if(bdir(2)/=0.d0)
then
2804 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2809 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2810 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2813 if(bdir(1)/=0.d0)
then
2814 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2818 if(bdir(2)/=0.d0)
then
2819 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2823 if(bdir(3)/=0.d0)
then
2824 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2829 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2830 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2832 if(lts(ix^d)>trac_delta)
then
2833 block%special_values(1)=max(block%special_values(1),te(ix^d))
2836 block%special_values(2)=tmax_local
2855 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2856 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2857 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2861 {
do ix^db=ixpmin^db,ixpmax^db\}
2862 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2864 if(bdir(1)/=0.d0)
then
2865 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2869 if(bdir(2)/=0.d0)
then
2870 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2876 if(bdir(1)/=0.d0)
then
2877 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2881 if(bdir(2)/=0.d0)
then
2882 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2886 if(bdir(3)/=0.d0)
then
2887 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2893 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2895 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2896 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2899 {
do ix^db=ixpmin^db,ixpmax^db\}
2901 if(w(ix^d,iw_mag(1))/=0.d0)
then
2902 bunitvec(1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2)
2906 if(w(ix^d,iw_mag(2))/=0.d0)
then
2907 bunitvec(2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2)
2913 if(w(ix^d,iw_mag(1))/=0.d0)
then
2914 bunitvec(1)=sign(1.d0,w(ix^d,iw_mag(1)))/dsqrt(1.d0+(w(ix^d,iw_mag(2))/w(ix^d,iw_mag(1)))**2+&
2915 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
2919 if(w(ix^d,iw_mag(2))/=0.d0)
then
2920 bunitvec(2)=sign(1.d0,w(ix^d,iw_mag(2)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(2)))**2+&
2921 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
2925 if(w(ix^d,iw_mag(3))/=0.d0)
then
2926 bunitvec(3)=sign(1.d0,w(ix^d,iw_mag(3)))/dsqrt(1.d0+(w(ix^d,iw_mag(1))/w(ix^d,iw_mag(3)))**2+&
2927 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
2933 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2935 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2936 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2942 {
do ix^db=ixpmin^db,ixpmax^db\}
2944 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
2945 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
2946 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2949 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
2950 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
2951 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
2952 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2958 call mpistop(
"unknown mhd_trac_type")
2961 end subroutine mhd_get_tcutoff
2964 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2967 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2968 double precision,
intent(in) :: wprim(ixi^s, nw)
2969 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2970 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2972 double precision :: csound(ixi^s,
ndim)
2973 double precision,
allocatable :: tmp(:^
d&)
2974 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
2978 allocate(tmp(ixa^s))
2980 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
2981 csound(ixa^s,id)=tmp(ixa^s)
2984 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
2985 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
2986 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
2987 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))
2991 ixamax^
d=ixcmax^
d+
kr(id,^
d);
2992 ixamin^
d=ixcmin^
d+
kr(id,^
d);
2993 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)))
2994 ixamax^
d=ixcmax^
d-
kr(id,^
d);
2995 ixamin^
d=ixcmin^
d-
kr(id,^
d);
2996 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)))
3001 ixamax^
d=jxcmax^
d+
kr(id,^
d);
3002 ixamin^
d=jxcmin^
d+
kr(id,^
d);
3003 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)))
3004 ixamax^
d=jxcmax^
d-
kr(id,^
d);
3005 ixamin^
d=jxcmin^
d-
kr(id,^
d);
3006 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)))
3010 end subroutine mhd_get_h_speed
3013 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3016 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3017 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3018 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3019 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3020 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3021 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3022 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3024 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3025 double precision :: umean, dmean, tmp1, tmp2, tmp3
3032 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3033 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3034 if(
present(cmin))
then
3035 {
do ix^db=ixomin^db,ixomax^db\}
3036 tmp1=sqrt(wlp(ix^
d,
rho_))
3037 tmp2=sqrt(wrp(ix^
d,
rho_))
3038 tmp3=1.d0/(tmp1+tmp2)
3039 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3040 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3041 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3042 cmin(ix^
d,1)=umean-dmean
3043 cmax(ix^
d,1)=umean+dmean
3045 if(h_correction)
then
3046 {
do ix^db=ixomin^db,ixomax^db\}
3047 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3048 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3052 {
do ix^db=ixomin^db,ixomax^db\}
3053 tmp1=sqrt(wlp(ix^d,
rho_))
3054 tmp2=sqrt(wrp(ix^d,
rho_))
3055 tmp3=1.d0/(tmp1+tmp2)
3056 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3057 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3058 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3059 cmax(ix^d,1)=abs(umean)+dmean
3063 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3064 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3065 if(
present(cmin))
then
3066 {
do ix^db=ixomin^db,ixomax^db\}
3067 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3068 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3070 if(h_correction)
then
3071 {
do ix^db=ixomin^db,ixomax^db\}
3072 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3073 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3077 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3081 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3082 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3083 if(
present(cmin))
then
3084 {
do ix^db=ixomin^db,ixomax^db\}
3085 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3086 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3087 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3089 if(h_correction)
then
3090 {
do ix^db=ixomin^db,ixomax^db\}
3091 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3092 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3096 {
do ix^db=ixomin^db,ixomax^db\}
3097 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3098 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3103 end subroutine mhd_get_cbounds
3106 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3109 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3110 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3111 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3112 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3113 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3114 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3115 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3117 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3122 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3123 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3125 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3126 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3128 if(
present(cmin))
then
3129 {
do ix^db=ixomin^db,ixomax^db\}
3130 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3131 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3132 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3135 {
do ix^db=ixomin^db,ixomax^db\}
3136 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3137 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3141 end subroutine mhd_get_cbounds_semirelati
3144 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3147 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3148 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3149 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3150 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3151 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3152 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3153 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3155 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3156 double precision :: umean, dmean, tmp1, tmp2, tmp3
3163 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3164 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3165 if(
present(cmin))
then
3166 {
do ix^db=ixomin^db,ixomax^db\}
3169 tmp3=1.d0/(tmp1+tmp2)
3170 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3171 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3172 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3173 cmin(ix^
d,1)=umean-dmean
3174 cmax(ix^
d,1)=umean+dmean
3176 if(h_correction)
then
3177 {
do ix^db=ixomin^db,ixomax^db\}
3178 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3179 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3183 {
do ix^db=ixomin^db,ixomax^db\}
3186 tmp3=1.d0/(tmp1+tmp2)
3187 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3188 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3189 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3190 cmax(ix^d,1)=abs(umean)+dmean
3194 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3195 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3196 if(
present(cmin))
then
3197 {
do ix^db=ixomin^db,ixomax^db\}
3198 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3199 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3201 if(h_correction)
then
3202 {
do ix^db=ixomin^db,ixomax^db\}
3203 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3204 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3208 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3212 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3213 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3214 if(
present(cmin))
then
3215 {
do ix^db=ixomin^db,ixomax^db\}
3216 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3217 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3218 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3220 if(h_correction)
then
3221 {
do ix^db=ixomin^db,ixomax^db\}
3222 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3223 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3227 {
do ix^db=ixomin^db,ixomax^db\}
3228 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3229 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3234 end subroutine mhd_get_cbounds_split_rho
3237 subroutine mhd_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3240 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3241 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3242 double precision,
intent(in) :: cmax(ixi^s)
3243 double precision,
intent(in),
optional :: cmin(ixi^s)
3244 type(ct_velocity),
intent(inout):: vcts
3246 integer :: idime,idimn
3252 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3254 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3256 if(.not.
allocated(vcts%vbarC))
then
3257 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3258 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3261 if(
present(cmin))
then
3262 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3263 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3265 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3266 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3269 idimn=mod(idim,
ndir)+1
3270 idime=mod(idim+1,
ndir)+1
3272 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3273 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3274 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3275 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3276 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3278 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3279 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3280 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3281 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3282 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3284 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
3287 end subroutine mhd_get_ct_velocity
3290 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3293 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3294 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3295 double precision,
intent(out):: csound(ixo^s)
3297 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3304 {
do ix^db=ixomin^db,ixomax^db \}
3305 inv_rho=1.d0/w(ix^
d,
rho_)
3312 cfast2=b2*inv_rho+csound(ix^
d)
3313 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3315 if(avmincs2<zero) avmincs2=zero
3316 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3318 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3322 {
do ix^db=ixomin^db,ixomax^db \}
3323 inv_rho=1.d0/w(ix^d,
rho_)
3329 b2=(^
c&w(ix^d,
b^
c_)**2+)
3330 cfast2=b2*inv_rho+csound(ix^d)
3331 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3332 if(avmincs2<zero) avmincs2=zero
3333 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3335 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3340 end subroutine mhd_get_csound_prim
3343 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3346 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3347 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3348 double precision,
intent(out):: csound(ixo^s)
3350 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3357 {
do ix^db=ixomin^db,ixomax^db \}
3364 cfast2=b2*inv_rho+csound(ix^
d)
3365 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3367 if(avmincs2<zero) avmincs2=zero
3368 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3370 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3374 {
do ix^db=ixomin^db,ixomax^db \}
3380 b2=(^
c&w(ix^d,
b^
c_)**2+)
3381 cfast2=b2*inv_rho+csound(ix^d)
3382 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3383 if(avmincs2<zero) avmincs2=zero
3384 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3386 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3391 end subroutine mhd_get_csound_prim_split
3394 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3397 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3399 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3400 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3402 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3405 {
do ix^db=ixomin^db,ixomax^db\}
3406 inv_rho = 1.d0/w(ix^
d,
rho_)
3409 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3410 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3411 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3412 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3415 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3416 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3417 if(avmincs2<zero) avmincs2=zero
3419 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3422 end subroutine mhd_get_csound_semirelati
3425 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3428 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3430 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3431 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3433 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3436 {
do ix^db=ixomin^db,ixomax^db\}
3437 inv_rho = 1.d0/w(ix^
d,
rho_)
3440 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3441 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3442 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3443 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3446 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3447 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3448 if(avmincs2<zero) avmincs2=zero
3450 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3453 end subroutine mhd_get_csound_semirelati_noe
3456 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3459 integer,
intent(in) :: ixi^
l, ixo^
l
3460 double precision,
intent(in) :: w(ixi^s,nw)
3461 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3462 double precision,
intent(out):: pth(ixi^s)
3470 end subroutine mhd_get_pthermal_noe
3473 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3477 integer,
intent(in) :: ixi^
l, ixo^
l
3478 double precision,
intent(in) :: w(ixi^s,nw)
3479 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3480 double precision,
intent(out):: pth(ixi^s)
3484 {
do ix^db= ixomin^db,ixomax^db\}
3488 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3493 if(check_small_values.and..not.fix_small_values)
then
3494 {
do ix^db= ixomin^db,ixomax^db\}
3495 if(pth(ix^d)<small_pressure)
then
3496 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3497 " encountered when call mhd_get_pthermal_inte"
3498 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3499 write(*,*)
"Location: ", x(ix^d,:)
3500 write(*,*)
"Cell number: ", ix^d
3502 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3505 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3506 write(*,*)
"Saving status at the previous time step"
3512 end subroutine mhd_get_pthermal_inte
3515 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3519 integer,
intent(in) :: ixi^
l, ixo^
l
3520 double precision,
intent(in) :: w(ixi^s,nw)
3521 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3522 double precision,
intent(out):: pth(ixi^s)
3526 {
do ix^db=ixomin^db,ixomax^db\}
3531 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3532 +(^
c&w(ix^
d,
b^
c_)**2+)))
3537 if(check_small_values.and..not.fix_small_values)
then
3538 {
do ix^db=ixomin^db,ixomax^db\}
3539 if(pth(ix^d)<small_pressure)
then
3540 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3541 " encountered when call mhd_get_pthermal"
3542 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3543 write(*,*)
"Location: ", x(ix^d,:)
3544 write(*,*)
"Cell number: ", ix^d
3546 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3549 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3550 write(*,*)
"Saving status at the previous time step"
3556 end subroutine mhd_get_pthermal_origin
3559 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3563 integer,
intent(in) :: ixi^
l, ixo^
l
3564 double precision,
intent(in) :: w(ixi^s,nw)
3565 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3566 double precision,
intent(out):: pth(ixi^s)
3568 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3571 {
do ix^db=ixomin^db,ixomax^db\}
3572 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3573 if(b2>smalldouble)
then
3581 inv_rho=1.d0/w(ix^
d,
rho_)
3583 b2=b2*inv_rho*inv_squared_c
3585 gamma2=1.d0/(1.d0+b2)
3587 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3591 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3592 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3593 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3597 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3603 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3604 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3605 +(^
c&w(ix^
d,
b^
c_)**2+)&
3606 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3610 if(check_small_values.and..not.fix_small_values)
then
3611 {
do ix^db=ixomin^db,ixomax^db\}
3612 if(pth(ix^d)<small_pressure)
then
3613 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3614 " encountered when call mhd_get_pthermal_semirelati"
3615 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3616 write(*,*)
"Location: ", x(ix^d,:)
3617 write(*,*)
"Cell number: ", ix^d
3619 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3622 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3623 write(*,*)
"Saving status at the previous time step"
3629 end subroutine mhd_get_pthermal_semirelati
3632 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3636 integer,
intent(in) :: ixi^
l, ixo^
l
3637 double precision,
intent(in) :: w(ixi^s,nw)
3638 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3639 double precision,
intent(out):: pth(ixi^s)
3643 {
do ix^db= ixomin^db,ixomax^db\}
3644 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3647 if(check_small_values.and..not.fix_small_values)
then
3648 {
do ix^db= ixomin^db,ixomax^db\}
3649 if(pth(ix^d)<small_pressure)
then
3650 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3651 " encountered when call mhd_get_pthermal_hde"
3652 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3653 write(*,*)
"Location: ", x(ix^d,:)
3654 write(*,*)
"Cell number: ", ix^d
3656 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3659 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3660 write(*,*)
"Saving status at the previous time step"
3666 end subroutine mhd_get_pthermal_hde
3669 subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3671 integer,
intent(in) :: ixi^
l, ixo^
l
3672 double precision,
intent(in) :: w(ixi^s, 1:nw)
3673 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3674 double precision,
intent(out):: res(ixi^s)
3675 res(ixo^s) = w(ixo^s,
te_)
3676 end subroutine mhd_get_temperature_from_te
3679 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3681 integer,
intent(in) :: ixi^
l, ixo^
l
3682 double precision,
intent(in) :: w(ixi^s, 1:nw)
3683 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3684 double precision,
intent(out):: res(ixi^s)
3686 double precision :: r(ixi^s)
3688 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3689 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3690 end subroutine mhd_get_temperature_from_eint
3693 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3695 integer,
intent(in) :: ixi^
l, ixo^
l
3696 double precision,
intent(in) :: w(ixi^s, 1:nw)
3697 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3698 double precision,
intent(out):: res(ixi^s)
3700 double precision :: r(ixi^s)
3702 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3704 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3706 end subroutine mhd_get_temperature_from_etot
3708 subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
3710 integer,
intent(in) :: ixi^
l, ixo^
l
3711 double precision,
intent(in) :: w(ixi^s, 1:nw)
3712 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3713 double precision,
intent(out):: res(ixi^s)
3715 double precision :: r(ixi^s)
3717 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3721 end subroutine mhd_get_temperature_from_etot_with_equi
3723 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3725 integer,
intent(in) :: ixi^
l, ixo^
l
3726 double precision,
intent(in) :: w(ixi^s, 1:nw)
3727 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3728 double precision,
intent(out):: res(ixi^s)
3730 double precision :: r(ixi^s)
3732 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3736 end subroutine mhd_get_temperature_from_eint_with_equi
3738 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3740 integer,
intent(in) :: ixi^
l, ixo^
l
3741 double precision,
intent(in) :: w(ixi^s, 1:nw)
3742 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3743 double precision,
intent(out):: res(ixi^s)
3745 double precision :: r(ixi^s)
3747 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3750 end subroutine mhd_get_temperature_equi
3752 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3754 integer,
intent(in) :: ixi^
l, ixo^
l
3755 double precision,
intent(in) :: w(ixi^s, 1:nw)
3756 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3757 double precision,
intent(out):: res(ixi^s)
3759 end subroutine mhd_get_rho_equi
3761 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3763 integer,
intent(in) :: ixi^
l, ixo^
l
3764 double precision,
intent(in) :: w(ixi^s, 1:nw)
3765 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3766 double precision,
intent(out):: res(ixi^s)
3768 end subroutine mhd_get_pe_equi
3771 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3775 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3777 double precision,
intent(in) :: wc(ixi^s,nw)
3779 double precision,
intent(in) :: w(ixi^s,nw)
3780 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3781 double precision,
intent(out) :: f(ixi^s,nwflux)
3783 double precision :: vhall(ixi^s,1:
ndir)
3784 double precision :: ptotal
3788 {
do ix^db=ixomin^db,ixomax^db\}
3801 {
do ix^db=ixomin^db,ixomax^db\}
3805 ^
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_)\
3806 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3808 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3811 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3812 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3814 ^
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_)\
3818 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3819 {
do ix^db=ixomin^db,ixomax^db\}
3820 if(total_energy)
then
3822 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3823 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3826 ^
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))\
3830 {
do ix^db=ixomin^db,ixomax^db\}
3831 f(ix^d,mag(idim))=w(ix^d,
psi_)
3833 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3838 {
do ix^db=ixomin^db,ixomax^db\}
3844 {
do ix^db=ixomin^db,ixomax^db\}
3845 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)
3850 end subroutine mhd_get_flux
3853 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3857 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3859 double precision,
intent(in) :: wc(ixi^s,nw)
3861 double precision,
intent(in) :: w(ixi^s,nw)
3862 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3863 double precision,
intent(out) :: f(ixi^s,nwflux)
3865 double precision :: vhall(ixi^s,1:
ndir)
3868 {
do ix^db=ixomin^db,ixomax^db\}
3879 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3880 {
do ix^db=ixomin^db,ixomax^db\}
3882 ^
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))\
3886 {
do ix^db=ixomin^db,ixomax^db\}
3887 f(ix^d,mag(idim))=w(ix^d,
psi_)
3889 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3894 {
do ix^db=ixomin^db,ixomax^db\}
3899 end subroutine mhd_get_flux_noe
3902 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3906 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3908 double precision,
intent(in) :: wc(ixi^s,nw)
3910 double precision,
intent(in) :: w(ixi^s,nw)
3911 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3912 double precision,
intent(out) :: f(ixi^s,nwflux)
3914 double precision :: vhall(ixi^s,1:
ndir)
3917 {
do ix^db=ixomin^db,ixomax^db\}
3930 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3931 {
do ix^db=ixomin^db,ixomax^db\}
3933 ^
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))\
3937 {
do ix^db=ixomin^db,ixomax^db\}
3938 f(ix^d,mag(idim))=w(ix^d,
psi_)
3940 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3945 {
do ix^db=ixomin^db,ixomax^db\}
3951 {
do ix^db=ixomin^db,ixomax^db\}
3952 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)
3957 end subroutine mhd_get_flux_hde
3960 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
3964 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3966 double precision,
intent(in) :: wc(ixi^s,nw)
3968 double precision,
intent(in) :: w(ixi^s,nw)
3969 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3970 double precision,
intent(out) :: f(ixi^s,nwflux)
3972 double precision :: vhall(ixi^s,1:
ndir)
3973 double precision :: ptotal, btotal(ixo^s,1:
ndir)
3976 {
do ix^db=ixomin^db,ixomax^db\}
3984 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
3988 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
3992 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
3993 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
3995 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
3999 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4002 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
4009 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
4010 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
4015 {
do ix^db=ixomin^db,ixomax^db\}
4016 f(ix^d,mag(idim))=w(ix^d,
psi_)
4018 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4023 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4024 {
do ix^db=ixomin^db,ixomax^db\}
4026 ^
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))\
4027 if(total_energy)
then
4029 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
4030 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
4036 {
do ix^db=ixomin^db,ixomax^db\}
4041 {
do ix^db=ixomin^db,ixomax^db\}
4042 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
4047 end subroutine mhd_get_flux_split
4050 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
4054 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4056 double precision,
intent(in) :: wc(ixi^s,nw)
4058 double precision,
intent(in) :: w(ixi^s,nw)
4059 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4060 double precision,
intent(out) :: f(ixi^s,nwflux)
4062 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
4065 {
do ix^db=ixomin^db,ixomax^db\}
4070 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4071 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4072 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4077 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4082 e2=(^
c&e(ix^
d,^
c)**2+)
4089 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4090 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4091 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4094 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4095 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4108 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4110 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)
4117 {
do ix^db=ixomin^db,ixomax^db\}
4118 f(ix^d,mag(idim))=w(ix^d,
psi_)
4120 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4125 {
do ix^db=ixomin^db,ixomax^db\}
4130 {
do ix^db=ixomin^db,ixomax^db\}
4131 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)
4136 end subroutine mhd_get_flux_semirelati
4138 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4142 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4144 double precision,
intent(in) :: wc(ixi^s,nw)
4146 double precision,
intent(in) :: w(ixi^s,nw)
4147 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4148 double precision,
intent(out) :: f(ixi^s,nwflux)
4150 double precision :: e(ixo^s,1:
ndir),e2
4153 {
do ix^db=ixomin^db,ixomax^db\}
4158 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4159 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4160 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4161 e2=(^
c&e(ix^
d,^
c)**2+)
4166 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4175 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4184 {
do ix^db=ixomin^db,ixomax^db\}
4185 f(ix^d,mag(idim))=w(ix^d,
psi_)
4187 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4192 {
do ix^db=ixomin^db,ixomax^db\}
4197 end subroutine mhd_get_flux_semirelati_noe
4203 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4205 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4206 double precision,
intent(in) :: qdt
4207 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4208 double precision,
intent(inout) :: w(ixi^s,1:nw)
4209 double precision :: tmp(ixi^s)
4210 double precision :: jxbxb(ixi^s,1:3)
4212 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4215 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4217 end subroutine add_source_ambipolar_internal_energy
4219 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4222 integer,
intent(in) :: ixi^
l, ixo^
l
4223 double precision,
intent(in) :: w(ixi^s,nw)
4224 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4225 double precision,
intent(out) :: res(:^
d&,:)
4227 double precision :: btot(ixi^s,1:3)
4228 double precision :: current(ixi^s,7-2*
ndir:3)
4229 double precision :: tmp(ixi^s),b2(ixi^s)
4230 integer :: idir, idirmin
4239 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4242 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4245 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4246 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4248 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4251 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4253 end subroutine mhd_get_jxbxb
4260 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4264 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4265 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4266 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4267 double precision,
intent(in) :: my_dt
4268 logical,
intent(in) :: fix_conserve_at_step
4270 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4271 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4272 double precision :: fe(ixi^s,
sdim:3)
4273 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4274 integer :: i, ixa^
l, ie_
4280 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4290 btot(ixa^s,1:3)=0.d0
4296 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4299 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4300 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4302 wres(ixo^s,
e_)=-tmp2(ixo^s)
4308 ff(ixa^s,1) = tmp(ixa^s,2)
4309 ff(ixa^s,2) = -tmp(ixa^s,1)
4311 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4312 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4313 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4316 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4318 ixamin^
d=ixomin^
d-1;
4319 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4328 ff(ixa^s,2) = tmp(ixa^s,3)
4329 ff(ixa^s,3) = -tmp(ixa^s,2)
4330 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4331 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4333 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4335 ff(ixa^s,1) = -tmp(ixa^s,3)
4337 ff(ixa^s,3) = tmp(ixa^s,1)
4338 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4339 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4340 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4344 ff(ixa^s,1) = tmp(ixa^s,2)
4345 ff(ixa^s,2) = -tmp(ixa^s,1)
4347 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4348 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4349 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4354 if(fix_conserve_at_step)
then
4355 fluxall=my_dt*fluxall
4362 end subroutine sts_set_source_ambipolar
4365 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4368 integer,
intent(in) :: ixi^
l, ixo^
l
4369 double precision,
intent(in) :: w(ixi^s,1:nw)
4370 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4372 double precision,
intent(in) :: ecc(ixi^s,1:3)
4373 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4374 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4376 integer :: hxc^
l,ixc^
l,ixa^
l
4377 integer :: idim1,idim2,idir,ix^
d
4383 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4385 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4386 ixamin^
d=ixcmin^
d+ix^
d;
4387 ixamax^
d=ixcmax^
d+ix^
d;
4388 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4390 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4396 ixcmin^d=ixomin^d-1;
4404 hxc^l=ixc^l-kr(idim2,^d);
4406 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4407 +lvc(idim1,idim2,idir)&
4412 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4415 end subroutine update_faces_ambipolar
4419 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4422 integer,
intent(in) :: ixi^
l, ixo^
l
4423 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4424 double precision,
intent(out) :: src(ixi^s)
4426 double precision :: ffc(ixi^s,1:
ndim)
4427 double precision :: dxinv(
ndim)
4428 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4434 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4436 ixbmin^
d=ixcmin^
d+ix^
d;
4437 ixbmax^
d=ixcmax^
d+ix^
d;
4440 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4442 ff(ixi^s,1:ndim)=0.d0
4444 ixb^l=ixo^l-kr(idims,^d);
4445 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4447 if({ ix^d==0 .and. ^d==idims | .or.})
then
4448 ixbmin^d=ixcmin^d-ix^d;
4449 ixbmax^d=ixcmax^d-ix^d;
4450 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4453 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4456 if(slab_uniform)
then
4458 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4459 ixb^l=ixo^l-kr(idims,^d);
4460 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4464 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4465 ixb^l=ixo^l-kr(idims,^d);
4466 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4468 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4470 end subroutine get_flux_on_cell_face
4474 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4477 integer,
intent(in) :: ixi^
l, ixo^
l
4478 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4479 double precision,
intent(in) :: w(ixi^s,1:nw)
4480 double precision :: dtnew
4482 double precision :: coef
4483 double precision :: dxarr(
ndim)
4484 double precision :: tmp(ixi^s)
4489 coef = maxval(abs(tmp(ixo^s)))
4496 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4498 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4501 end function get_ambipolar_dt
4509 integer,
intent(in) :: ixi^
l, ixo^
l
4510 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4511 double precision,
intent(inout) :: res(ixi^s)
4512 double precision :: tmp(ixi^s)
4513 double precision :: rho(ixi^s)
4522 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4526 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4533 integer,
intent(in) :: ixi^
l, ixo^
l
4534 double precision,
intent(in) :: qdt,dtfactor
4535 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4536 double precision,
intent(inout) :: w(ixi^s,1:nw)
4537 logical,
intent(in) :: qsourcesplit
4538 logical,
intent(inout) :: active
4545 if (.not. qsourcesplit)
then
4549 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4553 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4558 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4564 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4568 if (abs(
mhd_eta)>smalldouble)
then
4570 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4575 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4581 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4585 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4592 select case (type_divb)
4597 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4600 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4603 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4604 case (divb_janhunen)
4606 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4607 case (divb_lindejanhunen)
4609 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4610 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4611 case (divb_lindepowel)
4613 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4614 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4615 case (divb_lindeglm)
4617 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4618 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4619 case (divb_multigrid)
4624 call mpistop(
'Unknown divB fix')
4631 w,x,qsourcesplit,active,
rc_fl)
4641 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4650 if(.not.qsourcesplit)
then
4652 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4656 end subroutine mhd_add_source
4658 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4662 integer,
intent(in) :: ixi^
l, ixo^
l
4663 double precision,
intent(in) :: qdt,dtfactor
4664 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4665 double precision,
intent(inout) :: w(ixi^s,1:nw)
4666 double precision :: divv(ixi^s)
4682 end subroutine add_pe0_divv
4684 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4686 integer,
intent(in) :: ixi^
l,ixo^
l
4687 double precision,
intent(in) :: qdt
4688 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4689 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4690 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4692 double precision,
dimension(ixI^S) :: r,te,rho_loc
4693 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,bdir(
ndim),bunitvec(
ndim)
4697 call mhd_get_rfactor(wctprim,x,ixi^
l,ixi^
l,r)
4698 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*rho_loc(ixi^s))
4702 do ix2=ixomin2,ixomax2
4703 do ix1=ixomin1,ixomax1
4710 sigma_t7=sigma_t5*te(ix^
d)
4714 sigma_t7=sigma_t5*te(ix^
d)
4717 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))+
block%B0({ix^
d},^
d,0)\
4719 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))\
4721 if(bdir(1)/=0.d0)
then
4722 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
4726 if(bdir(2)/=0.d0)
then
4727 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
4731 sigmat5_bgradt=sigma_t5*(&
4732 bunitvec(1)*((8.d0*(te(ix1+1,ix2)-te(ix1-1,ix2))-te(ix1+2,ix2)+te(ix1-2,ix2))/12.d0)/
block%ds(ix^
d,1)&
4733 +bunitvec(2)*((8.d0*(te(ix1,ix2+1)-te(ix1,ix2-1))-te(ix1,ix2+2)+te(ix1,ix2-2))/12.d0)/
block%ds(ix^
d,2))
4735 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^
d)*(
mhd_gamma*wctprim(ix^
d,
p_)/rho_loc(ix^
d))**1.5d0)
4737 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^
d,
q_))/tau
4739 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(sigmat5_bgradt+wct(ix^
d,
q_))/&
4746 do ix3=ixomin3,ixomax3
4747 do ix2=ixomin2,ixomax2
4748 do ix1=ixomin1,ixomax1
4755 sigma_t7=sigma_t5*te(ix^
d)
4759 sigma_t7=sigma_t5*te(ix^
d)
4762 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))+
block%B0({ix^
d},^
d,0)\
4764 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))\
4766 if(bdir(1)/=0.d0)
then
4767 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
4771 if(bdir(2)/=0.d0)
then
4772 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
4776 if(bdir(3)/=0.d0)
then
4777 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
4781 sigmat5_bgradt=sigma_t5*(&
4782 bunitvec(1)*((8.d0*(te(ix1+1,ix2,ix3)-te(ix1-1,ix2,ix3))-te(ix1+2,ix2,ix3)+te(ix1-2,ix2,ix3))/12.d0)/
block%ds(ix^
d,1)&
4783 +bunitvec(2)*((8.d0*(te(ix1,ix2+1,ix3)-te(ix1,ix2-1,ix3))-te(ix1,ix2+2,ix3)+te(ix1,ix2-2,ix3))/12.d0)/
block%ds(ix^
d,2)&
4784 +bunitvec(3)*((8.d0*(te(ix1,ix2,ix3+1)-te(ix1,ix2,ix3-1))-te(ix1,ix2,ix3+2)+te(ix1,ix2,ix3-2))/12.d0)/
block%ds(ix^
d,3))
4786 f_sat=one/(one+abs(sigmat5_bgradt))/(1.5d0*rho_loc(ix^
d)*(
mhd_gamma*wctprim(ix^
d,
p_)/rho_loc(ix^
d))**1.5d0)
4788 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^
d,
q_))/tau
4790 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(sigmat5_bgradt+wct(ix^
d,
q_))/&
4797 end subroutine add_hypertc_source
4800 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4802 integer,
intent(in) :: ixi^
l, ixo^
l
4803 double precision,
intent(in) :: w(ixi^s,1:nw)
4804 double precision,
intent(inout) :: jxb(ixi^s,3)
4805 double precision :: a(ixi^s,3),
b(ixi^s,3)
4807 double precision :: current(ixi^s,7-2*
ndir:3)
4808 integer :: idir, idirmin
4813 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4817 b(ixo^s, idir) = w(ixo^s,mag(idir))
4826 a(ixo^s,idir)=current(ixo^s,idir)
4830 end subroutine get_lorentz_force
4834 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4836 integer,
intent(in) :: ixi^
l, ixo^
l
4837 double precision,
intent(in) :: w(ixi^s, nw)
4838 double precision,
intent(out) :: gamma_a2(ixo^s)
4839 double precision :: rho(ixi^s)
4845 rho(ixo^s) = w(ixo^s,
rho_)
4848 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4849 end subroutine mhd_gamma2_alfven
4853 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4855 integer,
intent(in) :: ixi^
l, ixo^
l
4856 double precision,
intent(in) :: w(ixi^s, nw)
4857 double precision :: gamma_a(ixo^s)
4859 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4860 gamma_a = sqrt(gamma_a)
4861 end function mhd_gamma_alfven
4865 integer,
intent(in) :: ixi^
l, ixo^
l
4866 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4867 double precision,
intent(out) :: rho(ixi^s)
4872 rho(ixo^s) = w(ixo^s,
rho_)
4878 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4881 integer,
intent(in) :: ixi^
l,ixo^
l, ie
4882 double precision,
intent(inout) :: w(ixi^s,1:nw)
4883 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4884 character(len=*),
intent(in) :: subname
4886 double precision :: rho(ixi^s)
4888 logical :: flag(ixi^s,1:nw)
4892 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4893 flag(ixo^s,ie)=.true.
4895 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4897 if(any(flag(ixo^s,ie)))
then
4901 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4904 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4910 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4913 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4919 end subroutine mhd_handle_small_ei
4921 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
4925 integer,
intent(in) :: ixi^
l, ixo^
l
4926 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4927 double precision,
intent(inout) :: w(ixi^s,1:nw)
4929 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
4938 end subroutine mhd_update_temperature
4941 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4944 integer,
intent(in) :: ixi^
l, ixo^
l
4945 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4946 double precision,
intent(inout) :: w(ixi^s,1:nw)
4948 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4960 a(ixo^s,idir)=
block%J0(ixo^s,idir)
4965 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4968 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4974 if(total_energy)
then
4977 b(ixo^s,:)=wct(ixo^s,mag(:))
4986 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
4989 axb(ixo^s,:)=axb(ixo^s,:)*qdt
4993 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
4997 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5002 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5007 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5009 end subroutine add_source_b0split
5012 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5016 integer,
intent(in) :: ixi^
l, ixo^
l
5017 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5018 double precision,
intent(inout) :: w(ixi^s,1:nw)
5019 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5021 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5022 integer :: idir, idirmin, ix^
d
5026 {
do ix^db=iximin^db,iximax^db\}
5028 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5029 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5030 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5032 call divvector(e,ixi^l,ixo^l,dive)
5034 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5037 {
do ix^db=ixomin^db,ixomax^db\}
5038 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5039 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5040 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5041 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5042 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5043 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5047 end subroutine add_source_semirelativistic
5050 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5054 integer,
intent(in) :: ixi^
l, ixo^
l
5055 double precision,
intent(in) :: qdt
5056 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5057 double precision,
intent(inout) :: w(ixi^s,1:nw)
5058 double precision,
intent(in) :: wctprim(ixi^s,1:nw)
5060 double precision :: divv(ixi^s), tmp
5072 {
do ix^db=ixomin^db,ixomax^db\}
5074 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5075 if(w(ix^
d,
e_)<small_e)
then
5080 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
5083 if(fix_small_values)
then
5084 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5086 end subroutine add_source_internal_e
5089 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5094 integer,
intent(in) :: ixi^
l, ixo^
l
5095 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5096 double precision,
intent(inout) :: w(ixi^s,1:nw)
5097 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5099 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5100 double precision :: current(ixi^s,7-2*
ndir:3)
5101 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5102 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5103 integer :: idir, idirmin, idims, ix^
d
5108 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5116 j(ixo^s,idir)=current(ixo^s,idir)
5134 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
5140 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5147 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)
5151 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5155 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5156 tmp(ixo^s)=sqrt(b2(ixo^s))
5157 where(tmp(ixo^s)>smalldouble)
5158 tmp(ixo^s)=1.d0/tmp(ixo^s)
5164 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5169 {
do ix^db=ixomin^db,ixomax^db\}
5171 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5173 b2(ix^
d)=vaoc/(1.d0+vaoc)
5176 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5179 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5183 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5186 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5187 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5190 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5193 end subroutine add_source_hydrodynamic_e
5199 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5204 integer,
intent(in) :: ixi^
l, ixo^
l
5205 double precision,
intent(in) :: qdt
5206 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5207 double precision,
intent(inout) :: w(ixi^s,1:nw)
5208 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5209 integer :: lxo^
l, kxo^
l
5211 double precision :: tmp(ixi^s),tmp2(ixi^s)
5214 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5215 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5224 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5225 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5232 gradeta(ixo^s,1:
ndim)=zero
5238 gradeta(ixo^s,idim)=tmp(ixo^s)
5245 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5252 tmp2(ixi^s)=bf(ixi^s,idir)
5254 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5255 jxo^
l=ixo^
l+
kr(idim,^
d);
5256 hxo^
l=ixo^
l-
kr(idim,^
d);
5257 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5258 tmp(ixo^s)=tmp(ixo^s)+&
5259 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5264 tmp2(ixi^s)=bf(ixi^s,idir)
5266 jxo^
l=ixo^
l+
kr(idim,^
d);
5267 hxo^
l=ixo^
l-
kr(idim,^
d);
5268 tmp(ixo^s)=tmp(ixo^s)+&
5269 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5274 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5278 do jdir=1,
ndim;
do kdir=idirmin,3
5279 if (
lvc(idir,jdir,kdir)/=0)
then
5280 if (
lvc(idir,jdir,kdir)==1)
then
5281 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5283 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5290 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5291 if(total_energy)
then
5292 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5298 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5301 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5303 end subroutine add_source_res1
5307 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5312 integer,
intent(in) :: ixi^
l, ixo^
l
5313 double precision,
intent(in) :: qdt
5314 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5315 double precision,
intent(inout) :: w(ixi^s,1:nw)
5318 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5319 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5320 integer :: ixa^
l,idir,idirmin,idirmin1
5324 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5325 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5335 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5340 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5349 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5352 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5357 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5359 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5361 if(total_energy)
then
5364 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5365 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5368 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5372 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5373 end subroutine add_source_res2
5377 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5381 integer,
intent(in) :: ixi^
l, ixo^
l
5382 double precision,
intent(in) :: qdt
5383 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5384 double precision,
intent(inout) :: w(ixi^s,1:nw)
5386 double precision :: current(ixi^s,7-2*
ndir:3)
5387 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5388 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5391 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5392 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5395 tmpvec(ixa^s,1:
ndir)=zero
5397 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5401 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5404 tmpvec(ixa^s,1:
ndir)=zero
5405 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5409 tmpvec2(ixa^s,1:
ndir)=zero
5410 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5413 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5416 if(total_energy)
then
5419 tmpvec2(ixa^s,1:
ndir)=zero
5420 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5421 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5422 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5423 end do;
end do;
end do
5425 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5426 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5429 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5431 end subroutine add_source_hyperres
5433 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5440 integer,
intent(in) :: ixi^
l, ixo^
l
5441 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5442 double precision,
intent(inout) :: w(ixi^s,1:nw)
5444 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5465 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5468 if(total_energy)
then
5477 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5486 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5490 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5492 end subroutine add_source_glm
5495 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5498 integer,
intent(in) :: ixi^
l, ixo^
l
5499 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5500 double precision,
intent(inout) :: w(ixi^s,1:nw)
5502 double precision :: divb(ixi^s), ba(1:
ndir)
5503 integer :: idir, ix^
d
5509 {
do ix^db=ixomin^db,ixomax^db\}
5514 if (total_energy)
then
5520 {
do ix^db=ixomin^db,ixomax^db\}
5522 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5524 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5525 if (total_energy)
then
5527 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5532 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5534 end subroutine add_source_powel
5536 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5541 integer,
intent(in) :: ixi^
l, ixo^
l
5542 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5543 double precision,
intent(inout) :: w(ixi^s,1:nw)
5545 double precision :: divb(ixi^s)
5546 integer :: idir, ix^
d
5551 {
do ix^db=ixomin^db,ixomax^db\}
5556 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5558 end subroutine add_source_janhunen
5560 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5565 integer,
intent(in) :: ixi^
l, ixo^
l
5566 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5567 double precision,
intent(inout) :: w(ixi^s,1:nw)
5569 double precision :: divb(ixi^s),graddivb(ixi^s)
5570 integer :: idim, idir, ixp^
l, i^
d, iside
5571 logical,
dimension(-1:1^D&) :: leveljump
5579 if(i^
d==0|.and.) cycle
5580 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5581 leveljump(i^
d)=.true.
5583 leveljump(i^
d)=.false.
5592 i^dd=kr(^dd,^d)*(2*iside-3);
5593 if (leveljump(i^dd))
then
5595 ixpmin^d=ixomin^d-i^d
5597 ixpmax^d=ixomax^d-i^d
5608 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5610 {
do i^db=ixpmin^db,ixpmax^db\}
5612 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5614 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5616 if (typedivbdiff==
'all' .and. total_energy)
then
5618 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5623 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5625 end subroutine add_source_linde
5632 integer,
intent(in) :: ixi^
l, ixo^
l
5633 double precision,
intent(in) :: w(ixi^s,1:nw)
5634 double precision :: divb(ixi^s), dsurface(ixi^s)
5636 double precision :: invb(ixo^s)
5637 integer :: ixa^
l,idims
5639 call get_divb(w,ixi^
l,ixo^
l,divb)
5641 where(invb(ixo^s)/=0.d0)
5642 invb(ixo^s)=1.d0/invb(ixo^s)
5645 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5647 ixamin^
d=ixomin^
d-1;
5648 ixamax^
d=ixomax^
d-1;
5649 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5651 ixa^
l=ixo^
l-
kr(idims,^
d);
5652 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5654 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5655 block%dvolume(ixo^s)/dsurface(ixo^s)
5666 integer,
intent(in) :: ixo^
l, ixi^
l
5667 double precision,
intent(in) :: w(ixi^s,1:nw)
5668 integer,
intent(out) :: idirmin
5671 double precision :: current(ixi^s,7-2*
ndir:3)
5672 integer :: idir, idirmin0
5678 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5679 block%J0(ixo^s,idirmin0:3)
5683 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5691 integer,
intent(in) :: ixi^
l, ixo^
l
5692 double precision,
intent(inout) :: dtnew
5693 double precision,
intent(in) ::
dx^
d
5694 double precision,
intent(in) :: w(ixi^s,1:nw)
5695 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5697 double precision :: dxarr(
ndim)
5698 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5699 integer :: idirmin,idim
5713 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5716 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5742 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5749 end subroutine mhd_get_dt
5752 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5757 integer,
intent(in) :: ixi^
l, ixo^
l
5758 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5759 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5761 double precision :: tmp,tmp1,invr,cot
5763 integer :: mr_,mphi_
5764 integer :: br_,bphi_
5767 br_=mag(1); bphi_=mag(1)-1+
phi_
5772 {
do ix^db=ixomin^db,ixomax^db\}
5775 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5780 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5785 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5786 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5787 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5788 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5789 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5791 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5792 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5793 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5796 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5801 {
do ix^db=ixomin^db,ixomax^db\}
5803 if(local_timestep)
then
5804 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5809 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5815 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5818 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5819 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5823 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5829 cot=1.d0/tan(x(ix^d,2))
5833 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5834 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5836 if(.not.stagger_grid)
then
5837 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5839 tmp=tmp+wprim(ix^d,
psi_)*cot
5841 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5846 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5847 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5848 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5850 if(.not.stagger_grid)
then
5851 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5853 tmp=tmp+wprim(ix^d,
psi_)*cot
5855 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5858 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5859 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5860 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5861 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5862 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5864 if(.not.stagger_grid)
then
5865 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5866 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5867 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5868 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5869 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5876 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5879 end subroutine mhd_add_source_geom
5882 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5887 integer,
intent(in) :: ixi^
l, ixo^
l
5888 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5889 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5891 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
5893 integer :: mr_,mphi_
5894 integer :: br_,bphi_
5897 br_=mag(1); bphi_=mag(1)-1+
phi_
5902 {
do ix^db=ixomin^db,ixomax^db\}
5905 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5916 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
5917 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
5918 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5923 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5929 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
5930 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
5931 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
5932 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5933 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5934 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
5936 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5937 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5938 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5941 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
5942 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
5947 {
do ix^db=ixomin^db,ixomax^db\}
5949 if(local_timestep)
then
5950 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
5956 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5957 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5958 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5962 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
5969 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
5975 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
5978 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
5979 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
5980 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
5984 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
5990 cot=1.d0/tan(x(ix^d,2))
5994 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
5995 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
5997 if(.not.stagger_grid)
then
5998 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6000 tmp=tmp+wprim(ix^d,
psi_)*cot
6002 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6008 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6009 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
6010 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6011 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
6013 if(.not.stagger_grid)
then
6014 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6016 tmp=tmp+wprim(ix^d,
psi_)*cot
6018 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6021 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6022 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6023 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6024 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
6025 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6026 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6027 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
6029 if(.not.stagger_grid)
then
6030 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6031 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6032 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6033 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6034 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6041 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6044 end subroutine mhd_add_source_geom_semirelati
6047 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6052 integer,
intent(in) :: ixi^
l, ixo^
l
6053 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6054 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6056 double precision :: tmp,tmp1,tmp2,invr,cot
6058 integer :: mr_,mphi_
6059 integer :: br_,bphi_
6062 br_=mag(1); bphi_=mag(1)-1+
phi_
6067 {
do ix^db=ixomin^db,ixomax^db\}
6070 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6075 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6080 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6081 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6082 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6083 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6084 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6086 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6087 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6088 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6091 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6096 {
do ix^db=ixomin^db,ixomax^db\}
6098 if(local_timestep)
then
6099 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6103 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6104 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6107 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6111 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6112 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6113 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6115 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6116 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6121 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6127 cot=1.d0/tan(x(ix^d,2))
6132 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6133 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6134 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6136 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6137 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6140 if(.not.stagger_grid)
then
6142 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6143 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6145 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6148 tmp=tmp+wprim(ix^d,
psi_)*cot
6150 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6156 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6157 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6158 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6159 +(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)
6161 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6162 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6163 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6166 if(.not.stagger_grid)
then
6168 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6169 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6171 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6174 tmp=tmp+wprim(ix^d,
psi_)*cot
6176 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6180 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6181 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6182 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6183 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6184 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6185 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6186 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6187 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6188 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6190 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6191 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6192 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6193 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6194 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6197 if(.not.stagger_grid)
then
6199 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6200 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6201 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6202 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6203 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6204 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6205 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6206 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6207 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6209 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6210 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6211 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6212 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6213 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6221 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6224 end subroutine mhd_add_source_geom_split
6229 integer,
intent(in) :: ixi^
l, ixo^
l
6230 double precision,
intent(in) :: w(ixi^s, nw)
6231 double precision :: mge(ixo^s)
6234 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6236 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6240 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6243 integer,
intent(in) :: ixi^
l, ixo^
l
6244 double precision,
intent(in) :: w(ixi^s,nw)
6245 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6246 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6248 double precision :: current(ixi^s,7-2*
ndir:3)
6249 double precision :: rho(ixi^s)
6250 integer :: idir, idirmin, ix^
d
6255 do idir = idirmin,
ndir
6256 {
do ix^db=ixomin^db,ixomax^db\}
6257 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6261 end subroutine mhd_getv_hall
6263 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6266 integer,
intent(in) :: ixi^
l, ixo^
l
6267 double precision,
intent(in) :: w(ixi^s,nw)
6268 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6269 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6272 double precision :: current(ixi^s,7-2*
ndir:3)
6273 integer :: idir, idirmin
6280 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6281 do idir = idirmin, 3
6285 end subroutine mhd_get_jambi
6287 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6290 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6291 double precision,
intent(in) :: qt
6292 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6293 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6296 double precision :: db(ixo^s), dpsi(ixo^s)
6300 {
do ix^db=ixomin^db,ixomax^db\}
6301 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6302 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6303 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6304 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6313 {
do ix^db=ixomin^db,ixomax^db\}
6314 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6315 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6316 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6317 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6318 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6320 if(total_energy)
then
6321 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6322 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6324 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6326 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6329 if(total_energy)
then
6330 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6331 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6336 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6338 end subroutine mhd_modify_wlr
6340 subroutine mhd_boundary_adjust(igrid,psb)
6342 integer,
intent(in) :: igrid
6345 integer :: ib, idims, iside, ixo^
l, i^
d
6354 i^
d=
kr(^
d,idims)*(2*iside-3);
6355 if (neighbor_type(i^
d,igrid)/=1) cycle
6356 ib=(idims-1)*2+iside
6374 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6379 end subroutine mhd_boundary_adjust
6381 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6384 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6385 double precision,
intent(inout) :: w(ixg^s,1:nw)
6386 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6388 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6389 integer :: ix^
d,ixf^
l
6402 do ix1=ixfmax1,ixfmin1,-1
6403 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6404 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6405 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6408 do ix1=ixfmax1,ixfmin1,-1
6409 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6410 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6411 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6412 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6413 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6414 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6415 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6429 do ix1=ixfmax1,ixfmin1,-1
6430 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6431 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6432 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6433 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6434 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6435 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6438 do ix1=ixfmax1,ixfmin1,-1
6439 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6440 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6441 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6442 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6443 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6444 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6445 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6446 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6447 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6448 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6449 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6450 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6451 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6452 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6453 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6454 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6455 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6456 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6470 do ix1=ixfmin1,ixfmax1
6471 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6472 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6473 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6476 do ix1=ixfmin1,ixfmax1
6477 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6478 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6479 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6480 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6481 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6482 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6483 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6497 do ix1=ixfmin1,ixfmax1
6498 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6499 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6500 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6501 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6502 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6503 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6506 do ix1=ixfmin1,ixfmax1
6507 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6508 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6509 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6510 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6511 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6512 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6513 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6514 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6515 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6516 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6517 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6518 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6519 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6520 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6521 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6522 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6523 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6524 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6538 do ix2=ixfmax2,ixfmin2,-1
6539 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6540 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6541 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6544 do ix2=ixfmax2,ixfmin2,-1
6545 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6546 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6547 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6548 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6549 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6550 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6551 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6565 do ix2=ixfmax2,ixfmin2,-1
6566 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6567 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6568 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6569 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6570 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6571 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6574 do ix2=ixfmax2,ixfmin2,-1
6575 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6576 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6577 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6578 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6579 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6580 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6581 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6582 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6583 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6584 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6585 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6586 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6587 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6588 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6589 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6590 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6591 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6592 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6606 do ix2=ixfmin2,ixfmax2
6607 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6608 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6609 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6612 do ix2=ixfmin2,ixfmax2
6613 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6614 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6615 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6616 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6617 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6618 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6619 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6633 do ix2=ixfmin2,ixfmax2
6634 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6635 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6636 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6637 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6638 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6639 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6642 do ix2=ixfmin2,ixfmax2
6643 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6644 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6645 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6646 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6647 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6648 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6649 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6650 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6651 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6652 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6653 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6654 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6655 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6656 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6657 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6658 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6659 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6660 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6677 do ix3=ixfmax3,ixfmin3,-1
6678 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6679 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6680 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6681 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6682 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6683 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6686 do ix3=ixfmax3,ixfmin3,-1
6687 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6688 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6689 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6690 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6691 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6692 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6693 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6694 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6695 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6696 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6697 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6698 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6699 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6700 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6701 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6702 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6703 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6704 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6719 do ix3=ixfmin3,ixfmax3
6720 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6721 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6722 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6723 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6724 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6725 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6728 do ix3=ixfmin3,ixfmax3
6729 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6730 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6731 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6732 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6733 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6734 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6735 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6736 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6737 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6738 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6739 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6740 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6741 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6742 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6743 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6744 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6745 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6746 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6752 call mpistop(
"Special boundary is not defined for this region")
6755 end subroutine fixdivb_boundary
6764 double precision,
intent(in) :: qdt
6765 double precision,
intent(in) :: qt
6766 logical,
intent(inout) :: active
6769 integer,
parameter :: max_its = 50
6770 double precision :: residual_it(max_its), max_divb
6771 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6772 double precision :: res
6773 double precision,
parameter :: max_residual = 1
d-3
6774 double precision,
parameter :: residual_reduction = 1
d-10
6775 integer :: iigrid, igrid
6776 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6779 mg%operator_type = mg_laplacian
6787 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6788 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6791 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6792 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6794 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6795 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6798 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6799 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6803 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6804 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6805 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6813 do iigrid = 1, igridstail
6814 igrid = igrids(iigrid);
6817 lvl =
mg%boxes(id)%lvl
6818 nc =
mg%box_size_lvl(lvl)
6824 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6826 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6827 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6832 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6835 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6836 if (
mype == 0) print *,
"iteration vs residual"
6839 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6840 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6841 if (residual_it(n) < residual_reduction * max_divb)
exit
6843 if (
mype == 0 .and. n > max_its)
then
6844 print *,
"divb_multigrid warning: not fully converged"
6845 print *,
"current amplitude of divb: ", residual_it(max_its)
6846 print *,
"multigrid smallest grid: ", &
6847 mg%domain_size_lvl(:,
mg%lowest_lvl)
6848 print *,
"note: smallest grid ideally has <= 8 cells"
6849 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6850 print *,
"note: dx/dy/dz should be similar"
6854 call mg_fas_vcycle(
mg, max_res=res)
6855 if (res < max_residual)
exit
6857 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6862 do iigrid = 1, igridstail
6863 igrid = igrids(iigrid);
6872 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6876 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6878 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
6880 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6893 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6894 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6897 if(total_energy)
then
6899 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6902 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6911 subroutine mhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
6914 integer,
intent(in) :: ixi^
l, ixo^
l
6915 double precision,
intent(in) :: qt,qdt
6917 double precision,
intent(in) :: wprim(ixi^s,1:nw)
6918 type(state) :: sct, s
6919 type(ct_velocity) :: vcts
6920 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6921 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6925 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
6927 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
6929 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
6931 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
6934 end subroutine mhd_update_faces
6937 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
6941 integer,
intent(in) :: ixi^
l, ixo^
l
6942 double precision,
intent(in) :: qt, qdt
6943 type(state) :: sct, s
6944 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6945 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6947 double precision :: circ(ixi^s,1:
ndim)
6949 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6950 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
6951 integer :: idim1,idim2,idir,iwdim1,iwdim2
6953 associate(bfaces=>s%ws,x=>s%x)
6960 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
6967 i1kr^
d=
kr(idim1,^
d);
6970 i2kr^
d=
kr(idim2,^
d);
6973 if (
lvc(idim1,idim2,idir)==1)
then
6975 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6977 {
do ix^db=ixcmin^db,ixcmax^db\}
6978 fe(ix^
d,idir)=quarter*&
6979 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
6980 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
6982 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
6987 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
6995 if(
associated(usr_set_electric_field)) &
6996 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
6998 circ(ixi^s,1:ndim)=zero
7003 ixcmin^d=ixomin^d-kr(idim1,^d);
7005 ixa^l=ixc^l-kr(idim2,^d);
7008 if(lvc(idim1,idim2,idir)==1)
then
7010 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7013 else if(lvc(idim1,idim2,idir)==-1)
then
7015 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7021 {
do ix^db=ixcmin^db,ixcmax^db\}
7023 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7025 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7032 end subroutine update_faces_average
7035 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7040 integer,
intent(in) :: ixi^
l, ixo^
l
7041 double precision,
intent(in) :: qt, qdt
7043 double precision,
intent(in) :: wp(ixi^s,1:nw)
7044 type(state) :: sct, s
7045 type(ct_velocity) :: vcts
7046 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7047 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7049 double precision :: circ(ixi^s,1:
ndim)
7051 double precision :: ecc(ixi^s,
sdim:3)
7052 double precision :: ein(ixi^s,
sdim:3)
7054 double precision :: el(ixi^s),er(ixi^s)
7056 double precision :: elc,erc
7058 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7060 double precision :: jce(ixi^s,
sdim:3)
7062 double precision :: xs(ixgs^t,1:
ndim)
7063 double precision :: gradi(ixgs^t)
7064 integer :: ixc^
l,ixa^
l
7065 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7067 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7070 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
7076 {
do ix^db=iximin^db,iximax^db\}
7079 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_)
7080 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_)
7081 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_)
7084 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7091 {
do ix^db=iximin^db,iximax^db\}
7094 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7095 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7096 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7099 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7113 i1kr^d=kr(idim1,^d);
7116 i2kr^d=kr(idim2,^d);
7119 if (lvc(idim1,idim2,idir)==1)
then
7121 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7124 {
do ix^db=ixcmin^db,ixcmax^db\}
7125 fe(ix^d,idir)=quarter*&
7126 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7127 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7132 ixamax^d=ixcmax^d+i1kr^d;
7133 {
do ix^db=ixamin^db,ixamax^db\}
7134 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7135 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7138 do ix^db=ixcmin^db,ixcmax^db\}
7139 if(vnorm(ix^d,idim1)>0.d0)
then
7141 else if(vnorm(ix^d,idim1)<0.d0)
then
7142 elc=el({ix^d+i1kr^d})
7144 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7146 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7148 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7149 erc=er({ix^d+i1kr^d})
7151 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7153 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7158 ixamax^d=ixcmax^d+i2kr^d;
7159 {
do ix^db=ixamin^db,ixamax^db\}
7160 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7161 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7164 do ix^db=ixcmin^db,ixcmax^db\}
7165 if(vnorm(ix^d,idim2)>0.d0)
then
7167 else if(vnorm(ix^d,idim2)<0.d0)
then
7168 elc=el({ix^d+i2kr^d})
7170 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7172 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7174 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7175 erc=er({ix^d+i2kr^d})
7177 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7179 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7183 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7188 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7202 if (lvc(idim1,idim2,idir)==0) cycle
7204 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7205 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7208 xs(ixa^s,:)=x(ixa^s,:)
7209 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7210 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7211 if (lvc(idim1,idim2,idir)==1)
then
7212 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7214 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7221 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7223 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7227 {
do ix^db=ixomin^db,ixomax^db\}
7228 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7229 +ein(ix1,ix2-1,ix3-1,idir))
7230 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7231 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7233 else if(idir==2)
then
7234 {
do ix^db=ixomin^db,ixomax^db\}
7235 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7236 +ein(ix1-1,ix2,ix3-1,idir))
7237 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7238 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7241 {
do ix^db=ixomin^db,ixomax^db\}
7242 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7243 +ein(ix1-1,ix2-1,ix3,idir))
7244 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7245 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7251 {
do ix^db=ixomin^db,ixomax^db\}
7252 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7253 +ein(ix1-1,ix2-1,idir))
7254 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7255 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7260 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7266 if(
associated(usr_set_electric_field)) &
7267 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7269 circ(ixi^s,1:ndim)=zero
7274 ixcmin^d=ixomin^d-kr(idim1,^d);
7276 ixa^l=ixc^l-kr(idim2,^d);
7279 if(lvc(idim1,idim2,idir)==1)
then
7281 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7284 else if(lvc(idim1,idim2,idir)==-1)
then
7286 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7292 {
do ix^db=ixcmin^db,ixcmax^db\}
7294 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7296 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7303 end subroutine update_faces_contact
7306 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
7311 integer,
intent(in) :: ixi^
l, ixo^
l
7312 double precision,
intent(in) :: qt, qdt
7313 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7314 type(state) :: sct, s
7315 type(ct_velocity) :: vcts
7317 double precision :: vtill(ixi^s,2)
7318 double precision :: vtilr(ixi^s,2)
7319 double precision :: bfacetot(ixi^s,
ndim)
7320 double precision :: btill(ixi^s,
ndim)
7321 double precision :: btilr(ixi^s,
ndim)
7322 double precision :: cp(ixi^s,2)
7323 double precision :: cm(ixi^s,2)
7324 double precision :: circ(ixi^s,1:
ndim)
7326 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7327 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7328 integer :: idim1,idim2,idir,ix^
d
7330 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7331 cbarmax=>vcts%cbarmax)
7344 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
7360 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7364 idim2=mod(idir+1,3)+1
7366 jxc^
l=ixc^
l+
kr(idim1,^
d);
7367 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7371 vtill(ixi^s,2),vtilr(ixi^s,2))
7374 vtill(ixi^s,1),vtilr(ixi^s,1))
7380 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7381 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7383 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7384 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7387 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7390 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7394 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7395 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7397 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7398 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7402 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7403 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7404 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7405 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7406 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7407 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7408 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7409 /(cp(ixc^s,2)+cm(ixc^s,2))
7412 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7416 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7430 circ(ixi^s,1:
ndim)=zero
7435 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7439 if(
lvc(idim1,idim2,idir)/=0)
then
7440 hxc^
l=ixc^
l-
kr(idim2,^
d);
7442 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7443 +
lvc(idim1,idim2,idir)&
7449 {
do ix^db=ixcmin^db,ixcmax^db\}
7451 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7453 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7459 end subroutine update_faces_hll
7462 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
7467 integer,
intent(in) :: ixi^
l, ixo^
l
7468 type(state),
intent(in) :: sct, s
7470 double precision :: jce(ixi^s,
sdim:3)
7473 double precision :: jcc(ixi^s,7-2*
ndir:3)
7475 double precision :: xs(ixgs^t,1:
ndim)
7477 double precision :: eta(ixi^s)
7478 double precision :: gradi(ixgs^t)
7479 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7481 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7487 if (
lvc(idim1,idim2,idir)==0) cycle
7489 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7490 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7493 xs(ixb^s,:)=x(ixb^s,:)
7494 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7495 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7496 if (
lvc(idim1,idim2,idir)==1)
then
7497 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7499 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7506 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7514 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7515 jcc(ixc^s,idir)=0.d0
7517 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7518 ixamin^
d=ixcmin^
d+ix^
d;
7519 ixamax^
d=ixcmax^
d+ix^
d;
7520 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7522 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7523 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7528 end subroutine get_resistive_electric_field
7531 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7534 integer,
intent(in) :: ixi^
l, ixo^
l
7535 double precision,
intent(in) :: w(ixi^s,1:nw)
7536 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7537 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7539 double precision :: jxbxb(ixi^s,1:3)
7540 integer :: idir,ixa^
l,ixc^
l,ix^
d
7543 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7550 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7553 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7554 ixamin^
d=ixcmin^
d+ix^
d;
7555 ixamax^
d=ixcmax^
d+ix^
d;
7556 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7558 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7561 end subroutine get_ambipolar_electric_field
7567 integer,
intent(in) :: ixo^
l
7577 do ix^db=ixomin^db,ixomax^db\}
7579 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7580 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7581 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7582 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7583 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7584 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7587 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7588 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7589 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7590 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7633 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7634 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7635 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7637 double precision :: adummy(ixis^s,1:3)
7643 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7646 integer,
intent(in) :: ixi^
l, ixo^
l
7647 double precision,
intent(in) :: w(ixi^s,1:nw)
7648 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7649 double precision,
intent(out):: rfactor(ixi^s)
7651 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7655 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)
7657 end subroutine rfactor_from_temperature_ionization
7659 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7661 integer,
intent(in) :: ixi^
l, ixo^
l
7662 double precision,
intent(in) :: w(ixi^s,1:nw)
7663 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7664 double precision,
intent(out):: rfactor(ixi^s)
7668 end subroutine rfactor_from_constant_ionization
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
subroutine cak_get_dt(w, ixil, ixol, dtnew, dxd, x)
Check time step for total radiation contribution.
subroutine cak_init(phys_gamma)
Initialize the module.
subroutine cak_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine reconstruct(ixil, ixcl, idir, q, ql, qr)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Module for flux conservation near refinement boundaries.
subroutine, public store_flux(igrid, fc, idimlim, nwfluxin)
subroutine, public store_edge(igrid, ixil, fe, idimlim)
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvector(qvec, ixil, ixol, divq, nth_in)
integer, parameter cylindrical
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
pure subroutine cross_product(ixil, ixol, a, b, axb)
Cross product of two vectors.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
logical b0fieldalloccoarse
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision phys_trac_mask
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
double precision dt
global time step
integer ndir
Number of spatial dimensions (components) for vector variables.
double precision courantpar
The Courant (CFL) number used for the simulation.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_velocity
Physical scaling factor for velocity.
double precision c_norm
Normalised speed of light.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
spatial steps for all dimensions at all levels
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
integer, parameter bc_symm
logical phys_trac
Use TRAC for MHD or 1D HD.
logical need_global_cmax
need global maximal wave speed
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
logical fix_small_values
fix small values with average or replace methods
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision small_density
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer phys_trac_finegrid
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for including gravity in (magneto)hydrodynamics simulations.
subroutine gravity_get_dt(w, ixil, ixol, dtnew, dxd, x)
subroutine gravity_init()
Initialize the module.
subroutine gravity_add_source(qdt, ixil, ixol, wct, wctprim, w, x, energy, rhov, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
subroutine ionization_degree_init()
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
integer, public, protected c_
logical, public, protected mhd_gravity
Whether gravity is added.
logical, public has_equi_rho0
whether split off equilibrium density
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
logical, public, protected mhd_glm_extended
Whether extended GLM-MHD is used with additional sources.
character(len=std_len), public, protected type_ct
Method type of constrained transport.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
logical, public, protected mhd_hyperbolic_thermal_conduction
Whether thermal conduction is used.
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
subroutine, public mhd_e_to_ei(ixil, ixol, w, x)
Transform total energy to internal energy.
double precision, public mhd_adiab
The adiabatic constant.
logical, public, protected mhd_partial_ionization
Whether plasma is partially ionized.
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
double precision, public, protected rr
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
double precision, public mhd_gamma
The adiabatic index.
integer, public, protected mhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
logical, public, protected mhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected mhd_semirelativistic
Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved.
integer, public, protected mhd_divb_nth
Whether divB is computed with a fourth order approximation.
integer, public, protected q_
Index of the heat flux q.
integer, public, protected mhd_n_tracer
Number of tracer species.
integer, public, protected te_
Indices of temperature.
integer, public, protected m
integer, public equi_rho0_
equi vars indices in the stateequi_vars array
integer, public, protected mhd_trac_type
Which TRAC method is used.
logical, public, protected mhd_cak_force
Whether CAK radiation line force is activated.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
logical, public, protected mhd_hall
Whether Hall-MHD is used.
type(te_fluid), allocatable, public te_fl_mhd
type of fluid for thermal emission synthesis
logical, public, protected mhd_ambipolar
Whether Ambipolar term is used.
double precision, public hypertc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixil, ixol)
Compute 2 times total magnetic energy.
subroutine, public multiplyambicoef(ixil, ixol, res, w, x)
multiply res by the ambipolar coefficient The ambipolar coefficient is calculated as -mhd_eta_ambi/rh...
logical, public partial_energy
Whether an internal or hydrodynamic energy equation is used.
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
logical, public, protected mhd_viscosity
Whether viscosity is added.
double precision, public, protected mhd_reduced_c
Reduced speed of light for semirelativistic MHD: 2% of light speed.
logical, public, protected mhd_energy
Whether an energy equation is used.
logical, public, protected mhd_ambipolar_exp
Whether Ambipolar term is implemented explicitly.
logical, public, protected mhd_htc_sat
Wheterh saturation is considered for hyperbolic TC.
logical, public, protected mhd_glm
Whether GLM-MHD is used to control div B.
logical, public clean_initial_divb
clean initial divB
procedure(sub_convert), pointer, public mhd_to_conserved
double precision, public mhd_eta
The MHD resistivity.
logical, public divbwave
Add divB wave in Roe solver.
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
double precision, public, protected mhd_trac_mask
Height of the mask used in the TRAC method.
procedure(mask_subroutine), pointer, public usr_mask_ambipolar
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
procedure(sub_get_pthermal), pointer, public mhd_get_temperature
integer, public equi_pe0_
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
procedure(sub_convert), pointer, public mhd_to_primitive
logical, public has_equi_pe0
whether split off equilibrium thermal pressure
logical, public, protected mhd_dump_full_vars
whether dump full variables (when splitting is used) in a separate dat file
logical, public, protected mhd_particles
Whether particles module is added.
integer, public, protected b
subroutine, public mhd_face_to_center(ixol, s)
calculate cell-center values from face-center values
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
subroutine, public get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
double precision, public mhd_etah
Hall resistivity.
subroutine, public mhd_get_v(w, x, ixil, ixol, v)
Calculate v vector.
double precision, public mhd_eta_ambi
The MHD ambipolar coefficient.
logical, public, protected mhd_hydrodynamic_e
Whether hydrodynamic energy is solved instead of total energy.
subroutine, public mhd_phys_init()
logical, public, protected mhd_trac
Whether TRAC method is used.
logical, public, protected eq_state_units
type(rc_fluid), allocatable, public rc_fl
type of fluid for radiative cooling
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
integer, public, protected rho_
Index of the density (in the w array)
logical, public, protected b0field_forcefree
B0 field is force-free.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
integer, public, protected tweight_
logical, public, protected mhd_ambipolar_sts
Whether Ambipolar term is implemented using supertimestepping.
procedure(sub_get_pthermal), pointer, public mhd_get_pthermal
subroutine, public mhd_ei_to_e(ixil, ixol, w, x)
Transform internal energy to total energy.
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected mhd_4th_order
MHD fourth order.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine, public mhd_get_rho(w, x, ixil, ixol, rho)
integer, public, protected psi_
Indices of the GLM psi.
logical, public mhd_equi_thermal
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_add_source(qdt, dtfactor, ixil, ixol, wct, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rotating_frame_init()
Initialize the module.
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startvar, nflux, startwbc, nwbc, evolve_b)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the 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.