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')
474 allocate(start_indices(number_species),stop_indices(number_species))
481 mom(:) = var_set_momentum(
ndir)
487 e_ = var_set_energy()
496 mag(:) = var_set_bfield(
ndir)
500 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
516 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
521 te_ = var_set_auxvar(
'Te',
'Te')
530 stop_indices(1)=nwflux
561 allocate(iw_vector(nvector))
562 iw_vector(1) =
mom(1) - 1
563 iw_vector(2) = mag(1) - 1
566 if (.not.
allocated(flux_type))
then
567 allocate(flux_type(
ndir, nwflux))
568 flux_type = flux_default
569 else if (any(shape(flux_type) /= [
ndir, nwflux]))
then
570 call mpistop(
"phys_check error: flux_type has wrong shape")
573 if(nwflux>mag(
ndir))
then
575 flux_type(:,mag(
ndir)+1:nwflux)=flux_hll
580 flux_type(:,
psi_)=flux_special
582 flux_type(idir,mag(idir))=flux_special
586 flux_type(idir,mag(idir))=flux_tvdlf
592 phys_get_dt => mhd_get_dt
595 phys_get_cmax => mhd_get_cmax_semirelati
597 phys_get_cmax => mhd_get_cmax_semirelati_noe
601 phys_get_cmax => mhd_get_cmax_origin
603 phys_get_cmax => mhd_get_cmax_origin_noe
606 phys_get_a2max => mhd_get_a2max
607 phys_get_tcutoff => mhd_get_tcutoff
608 phys_get_h_speed => mhd_get_h_speed
610 phys_get_cbounds => mhd_get_cbounds_split_rho
612 phys_get_cbounds => mhd_get_cbounds_semirelati
614 phys_get_cbounds => mhd_get_cbounds
617 phys_to_primitive => mhd_to_primitive_hde
619 phys_to_conserved => mhd_to_conserved_hde
623 phys_to_primitive => mhd_to_primitive_semirelati
625 phys_to_conserved => mhd_to_conserved_semirelati
628 phys_to_primitive => mhd_to_primitive_semirelati_noe
630 phys_to_conserved => mhd_to_conserved_semirelati_noe
635 phys_to_primitive => mhd_to_primitive_split_rho
637 phys_to_conserved => mhd_to_conserved_split_rho
640 phys_to_primitive => mhd_to_primitive_inte
642 phys_to_conserved => mhd_to_conserved_inte
645 phys_to_primitive => mhd_to_primitive_origin
647 phys_to_conserved => mhd_to_conserved_origin
650 phys_to_primitive => mhd_to_primitive_origin_noe
652 phys_to_conserved => mhd_to_conserved_origin_noe
657 phys_get_flux => mhd_get_flux_hde
660 phys_get_flux => mhd_get_flux_semirelati
662 phys_get_flux => mhd_get_flux_semirelati_noe
666 phys_get_flux => mhd_get_flux_split
668 phys_get_flux => mhd_get_flux
670 phys_get_flux => mhd_get_flux_noe
675 phys_add_source_geom => mhd_add_source_geom_semirelati
677 phys_add_source_geom => mhd_add_source_geom_split
679 phys_add_source_geom => mhd_add_source_geom
681 phys_add_source => mhd_add_source
682 phys_check_params => mhd_check_params
683 phys_write_info => mhd_write_info
686 phys_handle_small_values => mhd_handle_small_values_inte
687 mhd_handle_small_values => mhd_handle_small_values_inte
688 phys_check_w => mhd_check_w_inte
690 phys_handle_small_values => mhd_handle_small_values_hde
691 mhd_handle_small_values => mhd_handle_small_values_hde
692 phys_check_w => mhd_check_w_hde
694 phys_handle_small_values => mhd_handle_small_values_semirelati
695 mhd_handle_small_values => mhd_handle_small_values_semirelati
696 phys_check_w => mhd_check_w_semirelati
698 phys_handle_small_values => mhd_handle_small_values_split
699 mhd_handle_small_values => mhd_handle_small_values_split
700 phys_check_w => mhd_check_w_split
702 phys_handle_small_values => mhd_handle_small_values_origin
703 mhd_handle_small_values => mhd_handle_small_values_origin
704 phys_check_w => mhd_check_w_origin
706 phys_handle_small_values => mhd_handle_small_values_noe
707 mhd_handle_small_values => mhd_handle_small_values_noe
708 phys_check_w => mhd_check_w_noe
712 phys_get_pthermal => mhd_get_pthermal_inte
715 phys_get_pthermal => mhd_get_pthermal_hde
718 phys_get_pthermal => mhd_get_pthermal_semirelati
721 phys_get_pthermal => mhd_get_pthermal_origin
724 phys_get_pthermal => mhd_get_pthermal_noe
729 phys_set_equi_vars => set_equi_vars_grid
732 if(type_divb==divb_glm)
then
733 phys_modify_wlr => mhd_modify_wlr
738 mhd_get_rfactor=>rfactor_from_temperature_ionization
739 phys_update_temperature => mhd_update_temperature
743 mhd_get_rfactor=>rfactor_from_constant_ionization
768 transverse_ghost_cells = 1
769 phys_get_ct_velocity => mhd_get_ct_velocity_average
770 phys_update_faces => mhd_update_faces_average
772 transverse_ghost_cells = 1
773 phys_get_ct_velocity => mhd_get_ct_velocity_contact
774 phys_update_faces => mhd_update_faces_contact
776 transverse_ghost_cells = 2
777 phys_get_ct_velocity => mhd_get_ct_velocity_hll
778 phys_update_faces => mhd_update_faces_hll
780 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
783 phys_modify_wlr => mhd_modify_wlr
785 phys_boundary_adjust => mhd_boundary_adjust
794 call mhd_physical_units()
800 call mpistop(
"thermal conduction needs mhd_energy=T")
803 call mpistop(
"hyperbolic thermal conduction needs mhd_energy=T")
806 call mpistop(
"radiative cooling needs mhd_energy=T")
817 if(phys_internal_e)
then
819 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
821 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
825 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot_with_equi
827 tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
831 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
833 tc_fl%has_equi = .true.
834 tc_fl%get_temperature_equi => mhd_get_temperature_equi
835 tc_fl%get_rho_equi => mhd_get_rho_equi
837 tc_fl%has_equi = .false.
840 tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
864 rc_fl%get_var_Rfactor => mhd_get_rfactor
868 rc_fl%has_equi = .true.
869 rc_fl%get_rho_equi => mhd_get_rho_equi
870 rc_fl%get_pthermal_equi => mhd_get_pe_equi
872 rc_fl%has_equi = .false.
878 te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
880 phys_te_images => mhd_te_images
896 if (particles_eta < zero) particles_eta =
mhd_eta
897 if (particles_etah < zero) particles_eta =
mhd_etah
899 write(*,*)
'*****Using particles: with mhd_eta, mhd_etah :',
mhd_eta,
mhd_etah
900 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
914 phys_wider_stencil = 2
916 phys_wider_stencil = 1
924 call add_sts_method(get_ambipolar_dt,sts_set_source_ambipolar,mag(1),&
936 phys_wider_stencil = 2
938 phys_wider_stencil = 1
952 subroutine mhd_te_images
957 case(
'EIvtiCCmpi',
'EIvtuCCmpi')
959 case(
'ESvtiCCmpi',
'ESvtuCCmpi')
961 case(
'SIvtiCCmpi',
'SIvtuCCmpi')
963 case(
'WIvtiCCmpi',
'WIvtuCCmpi')
966 call mpistop(
"Error in synthesize emission: Unknown convert_type")
968 end subroutine mhd_te_images
974 subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
978 integer,
intent(in) :: ixi^
l, ixo^
l, igrid, nflux
979 double precision,
intent(in) :: x(ixi^s,1:
ndim)
980 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
981 double precision,
intent(in) :: my_dt
982 logical,
intent(in) :: fix_conserve_at_step
984 end subroutine mhd_sts_set_source_tc_mhd
986 function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
993 integer,
intent(in) :: ixi^
l, ixo^
l
994 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
995 double precision,
intent(in) :: w(ixi^s,1:nw)
996 double precision :: dtnew
999 end function mhd_get_tc_dt_mhd
1001 subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
1004 integer,
intent(in) :: ixi^
l,ixo^
l
1005 double precision,
intent(inout) :: w(ixi^s,1:nw)
1006 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1007 integer,
intent(in) :: step
1008 character(len=140) :: error_msg
1010 write(error_msg,
"(a,i3)")
"Thermal conduction step ", step
1011 call mhd_handle_small_ei(w,x,ixi^
l,ixo^
l,
e_,error_msg)
1012 end subroutine mhd_tc_handle_small_e
1015 subroutine tc_params_read_mhd(fl)
1017 type(tc_fluid),
intent(inout) :: fl
1019 double precision :: tc_k_para=0d0
1020 double precision :: tc_k_perp=0d0
1023 logical :: tc_perpendicular=.false.
1024 logical :: tc_saturate=.false.
1025 character(len=std_len) :: tc_slope_limiter=
"MC"
1027 namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1031 read(
unitpar, tc_list,
end=111)
1035 fl%tc_perpendicular = tc_perpendicular
1036 fl%tc_saturate = tc_saturate
1037 fl%tc_k_para = tc_k_para
1038 fl%tc_k_perp = tc_k_perp
1039 select case(tc_slope_limiter)
1041 fl%tc_slope_limiter = 0
1044 fl%tc_slope_limiter = 1
1047 fl%tc_slope_limiter = 2
1050 fl%tc_slope_limiter = 3
1053 fl%tc_slope_limiter = 4
1055 call mpistop(
"Unknown tc_slope_limiter, choose MC, minmod")
1057 end subroutine tc_params_read_mhd
1061 subroutine rc_params_read(fl)
1064 type(rc_fluid),
intent(inout) :: fl
1066 double precision :: cfrac=0.1d0
1069 double precision :: rad_cut_hgt=0.5d0
1070 double precision :: rad_cut_dey=0.15d0
1073 integer :: ncool = 4000
1075 logical :: tfix=.false.
1077 logical :: rc_split=.false.
1078 logical :: rad_cut=.false.
1080 character(len=std_len) :: coolcurve=
'JCcorona'
1082 character(len=std_len) :: coolmethod=
'exact'
1084 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split,rad_cut,rad_cut_hgt,rad_cut_dey
1088 read(
unitpar, rc_list,
end=111)
1093 fl%coolcurve=coolcurve
1094 fl%coolmethod=coolmethod
1097 fl%rc_split=rc_split
1100 fl%rad_cut_hgt=rad_cut_hgt
1101 fl%rad_cut_dey=rad_cut_dey
1102 end subroutine rc_params_read
1106 subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1109 integer,
intent(in) :: igrid, ixi^
l, ixo^
l
1110 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1112 double precision :: delx(ixi^s,1:
ndim)
1113 double precision :: xc(ixi^s,1:
ndim),xshift^
d
1114 integer :: idims, ixc^
l, hxo^
l, ix, idims2
1120 delx(ixi^s,1:
ndim)=ps(igrid)%dx(ixi^s,1:
ndim)
1124 hxo^
l=ixo^
l-
kr(idims,^
d);
1130 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
1133 xshift^
d=half*(one-
kr(^
d,idims));
1140 xc(ix^
d%ixC^s,^
d)=x(ix^
d%ixC^s,^
d)+(half-xshift^
d)*delx(ix^
d%ixC^s,^
d)
1144 call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1147 end subroutine set_equi_vars_grid_faces
1150 subroutine set_equi_vars_grid(igrid)
1154 integer,
intent(in) :: igrid
1160 call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^
ll,
ixm^
ll)
1162 end subroutine set_equi_vars_grid
1165 function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc)
result(wnew)
1167 integer,
intent(in) :: ixi^
l,ixo^
l, nwc
1168 double precision,
intent(in) :: w(ixi^s, 1:nw)
1169 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1170 double precision :: wnew(ixo^s, 1:nwc)
1177 wnew(ixo^s,
mom(:))=w(ixo^s,
mom(:))
1183 wnew(ixo^s,mag(1:
ndir))=w(ixo^s,mag(1:
ndir))
1187 wnew(ixo^s,
e_)=w(ixo^s,
e_)
1191 if(
b0field .and. total_energy)
then
1192 wnew(ixo^s,
e_)=wnew(ixo^s,
e_)+0.5d0*sum(
block%B0(ixo^s,:,0)**2,dim=
ndim+1) &
1193 + sum(w(ixo^s,mag(:))*
block%B0(ixo^s,:,0),dim=
ndim+1)
1197 end function convert_vars_splitting
1199 subroutine mhd_check_params
1207 if (
mhd_gamma <= 0.0d0)
call mpistop (
"Error: mhd_gamma <= 0")
1208 if (
mhd_adiab < 0.0d0)
call mpistop (
"Error: mhd_adiab < 0")
1212 call mpistop (
"Error: mhd_gamma <= 0 or mhd_gamma == 1")
1213 inv_gamma_1=1.d0/gamma_1
1218 call mpistop(
"usr_set_equi_vars has to be implemented in the user file")
1223 if(
mype .eq. 0) print*,
" add conversion method: split -> full "
1228 end subroutine mhd_check_params
1230 subroutine mhd_physical_units()
1232 double precision :: mp,kb,miu0,c_lightspeed
1233 double precision :: a,
b
1244 c_lightspeed=const_c
1385 end subroutine mhd_physical_units
1387 subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1390 logical,
intent(in) :: primitive
1391 logical,
intent(inout) :: flag(ixi^s,1:nw)
1392 integer,
intent(in) :: ixi^
l, ixo^
l
1393 double precision,
intent(in) :: w(ixi^s,nw)
1395 double precision :: tmp,b2,
b(ixo^s,1:
ndir)
1396 double precision :: v(ixo^s,1:
ndir),gamma2,inv_rho
1407 {
do ix^db=ixomin^db,ixomax^db \}
1408 if(w(ix^
d,
e_) < small_e) flag(ix^
d,
e_) = .true.
1411 {
do ix^db=ixomin^db,ixomax^db \}
1412 b2=(^
c&w(ix^d,
b^
c_)**2+)
1413 if(b2>smalldouble)
then
1418 ^
c&
b(ix^d,^
c)=w(ix^d,
b^
c_)*tmp\
1419 tmp=(^
c&
b(ix^d,^
c)*w(ix^d,
m^
c_)+)
1420 inv_rho = 1d0/w(ix^d,
rho_)
1422 b2=b2*inv_rho*inv_squared_c
1424 gamma2=1.d0/(1.d0+b2)
1426 ^
c&v(ix^d,^
c)=gamma2*(w(ix^d,
m^
c_)+b2*
b(ix^d,^
c)*tmp*inv_rho)\
1429 b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
1430 b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
1431 b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1436 b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
1442 tmp=w(ix^d,
e_)-half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
1443 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
1444 if(tmp<small_e) flag(ix^d,
e_)=.true.
1450 end subroutine mhd_check_w_semirelati
1452 subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1455 logical,
intent(in) :: primitive
1456 integer,
intent(in) :: ixi^
l, ixo^
l
1457 double precision,
intent(in) :: w(ixi^s,nw)
1458 logical,
intent(inout) :: flag(ixi^s,1:nw)
1463 {
do ix^db=ixomin^db,ixomax^db\}
1469 (^
c&w(ix^
d,
b^
c_)**2+))<small_e) flag(ix^
d,
e_) = .true.
1473 end subroutine mhd_check_w_origin
1475 subroutine mhd_check_w_split(primitive,ixI^L,ixO^L,w,flag)
1478 logical,
intent(in) :: primitive
1479 integer,
intent(in) :: ixi^
l, ixo^
l
1480 double precision,
intent(in) :: w(ixi^s,nw)
1481 logical,
intent(inout) :: flag(ixi^s,1:nw)
1483 double precision :: tmp
1487 {
do ix^db=ixomin^db,ixomax^db\}
1493 tmp=w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/tmp+(^
c&w(ix^
d,
b^
c_)**2+))
1498 end subroutine mhd_check_w_split
1500 subroutine mhd_check_w_noe(primitive,ixI^L,ixO^L,w,flag)
1503 logical,
intent(in) :: primitive
1504 integer,
intent(in) :: ixi^
l, ixo^
l
1505 double precision,
intent(in) :: w(ixi^s,nw)
1506 logical,
intent(inout) :: flag(ixi^s,1:nw)
1511 {
do ix^db=ixomin^db,ixomax^db\}
1515 end subroutine mhd_check_w_noe
1517 subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1520 logical,
intent(in) :: primitive
1521 integer,
intent(in) :: ixi^
l, ixo^
l
1522 double precision,
intent(in) :: w(ixi^s,nw)
1523 logical,
intent(inout) :: flag(ixi^s,1:nw)
1528 {
do ix^db=ixomin^db,ixomax^db\}
1533 if(w(ix^
d,
e_)<small_e) flag(ix^
d,
e_) = .true.
1537 end subroutine mhd_check_w_inte
1539 subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1542 logical,
intent(in) :: primitive
1543 integer,
intent(in) :: ixi^
l, ixo^
l
1544 double precision,
intent(in) :: w(ixi^s,nw)
1545 logical,
intent(inout) :: flag(ixi^s,1:nw)
1550 {
do ix^db=ixomin^db,ixomax^db\}
1555 if(w(ix^
d,
e_)-half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)<small_e) flag(ix^
d,
e_) = .true.
1559 end subroutine mhd_check_w_hde
1562 subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1564 integer,
intent(in) :: ixi^
l, ixo^
l
1565 double precision,
intent(inout) :: w(ixi^s, nw)
1566 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1570 {
do ix^db=ixomin^db,ixomax^db\}
1572 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1574 +(^
c&w(ix^
d,
b^
c_)**2+))
1579 end subroutine mhd_to_conserved_origin
1582 subroutine mhd_to_conserved_origin_noe(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\}
1595 end subroutine mhd_to_conserved_origin_noe
1598 subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1600 integer,
intent(in) :: ixi^
l, ixo^
l
1601 double precision,
intent(inout) :: w(ixi^s, nw)
1602 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1606 {
do ix^db=ixomin^db,ixomax^db\}
1608 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1614 end subroutine mhd_to_conserved_hde
1617 subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1619 integer,
intent(in) :: ixi^
l, ixo^
l
1620 double precision,
intent(inout) :: w(ixi^s, nw)
1621 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1625 {
do ix^db=ixomin^db,ixomax^db\}
1627 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1632 end subroutine mhd_to_conserved_inte
1635 subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1637 integer,
intent(in) :: ixi^
l, ixo^
l
1638 double precision,
intent(inout) :: w(ixi^s, nw)
1639 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1641 double precision :: rho
1644 {
do ix^db=ixomin^db,ixomax^db\}
1647 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1648 +half*((^
c&w(ix^
d,
m^
c_)**2+)*rho&
1649 +(^
c&w(ix^
d,
b^
c_)**2+))
1654 end subroutine mhd_to_conserved_split_rho
1657 subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1659 integer,
intent(in) :: ixi^
l, ixo^
l
1660 double precision,
intent(inout) :: w(ixi^s, nw)
1661 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1663 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1666 {
do ix^db=ixomin^db,ixomax^db\}
1668 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1669 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1670 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1671 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1672 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1673 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1678 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1679 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1680 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1688 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1
1692 w(ix^
d,
e_)=w(ix^
d,
p_)*inv_gamma_1&
1694 +(^
c&w(ix^
d,
b^
c_)**2+)&
1695 +(^
c&e(ix^
d,^
c)**2+)*inv_squared_c)
1703 end subroutine mhd_to_conserved_semirelati
1705 subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
1707 integer,
intent(in) :: ixi^
l, ixo^
l
1708 double precision,
intent(inout) :: w(ixi^s, nw)
1709 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1711 double precision :: e(ixo^s,1:
ndir), s(ixo^s,1:
ndir)
1714 {
do ix^db=ixomin^db,ixomax^db\}
1716 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1717 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1718 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1719 s(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
1720 s(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
1721 s(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
1726 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1727 s(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
1728 s(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
1738 end subroutine mhd_to_conserved_semirelati_noe
1741 subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1743 integer,
intent(in) :: ixi^
l, ixo^
l
1744 double precision,
intent(inout) :: w(ixi^s, nw)
1745 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1747 double precision :: inv_rho
1752 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin')
1755 {
do ix^db=ixomin^db,ixomax^db\}
1756 inv_rho = 1.d0/w(ix^
d,
rho_)
1760 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1762 +(^
c&w(ix^
d,
b^
c_)**2+)))
1765 end subroutine mhd_to_primitive_origin
1768 subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
1770 integer,
intent(in) :: ixi^
l, ixo^
l
1771 double precision,
intent(inout) :: w(ixi^s, nw)
1772 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1774 double precision :: inv_rho
1779 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_origin_noe')
1782 {
do ix^db=ixomin^db,ixomax^db\}
1783 inv_rho = 1.d0/w(ix^
d,
rho_)
1788 end subroutine mhd_to_primitive_origin_noe
1791 subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1793 integer,
intent(in) :: ixi^
l, ixo^
l
1794 double precision,
intent(inout) :: w(ixi^s, nw)
1795 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1797 double precision :: inv_rho
1802 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_hde')
1805 {
do ix^db=ixomin^db,ixomax^db\}
1806 inv_rho = 1d0/w(ix^
d,
rho_)
1813 end subroutine mhd_to_primitive_hde
1816 subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1818 integer,
intent(in) :: ixi^
l, ixo^
l
1819 double precision,
intent(inout) :: w(ixi^s, nw)
1820 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1822 double precision :: inv_rho
1827 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_inte')
1830 {
do ix^db=ixomin^db,ixomax^db\}
1832 w(ix^
d,
p_)=w(ix^
d,
e_)*gamma_1
1834 inv_rho = 1.d0/w(ix^
d,
rho_)
1838 end subroutine mhd_to_primitive_inte
1841 subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1843 integer,
intent(in) :: ixi^
l, ixo^
l
1844 double precision,
intent(inout) :: w(ixi^s, nw)
1845 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1847 double precision :: inv_rho
1852 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_split_rho')
1855 {
do ix^db=ixomin^db,ixomax^db\}
1860 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1862 (^
c&w(ix^
d,
m^
c_)**2+)+(^
c&w(ix^
d,
b^
c_)**2+)))
1865 end subroutine mhd_to_primitive_split_rho
1868 subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1870 integer,
intent(in) :: ixi^
l, ixo^
l
1871 double precision,
intent(inout) :: w(ixi^s, nw)
1872 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1874 double precision ::
b(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
1879 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati')
1882 {
do ix^db=ixomin^db,ixomax^db\}
1883 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1884 if(b2>smalldouble)
then
1892 inv_rho=1.d0/w(ix^
d,
rho_)
1894 b2=b2*inv_rho*inv_squared_c
1896 gamma2=1.d0/(1.d0+b2)
1898 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1902 w(ix^
d,
p_)=gamma_1*w(ix^
d,
e_)
1906 b(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
1907 b(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
1908 b(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1912 b(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
1918 w(ix^
d,
p_)=gamma_1*(w(ix^
d,
e_)&
1920 +(^
c&w(ix^
d,
b^
c_)**2+)&
1921 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
1925 end subroutine mhd_to_primitive_semirelati
1928 subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
1930 integer,
intent(in) :: ixi^
l, ixo^
l
1931 double precision,
intent(inout) :: w(ixi^s, nw)
1932 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1934 double precision ::
b(ixo^s,1:
ndir),tmp,b2,gamma2,inv_rho
1935 integer :: ix^
d, idir
1939 call mhd_handle_small_values(.false., w, x, ixi^
l, ixo^
l,
'mhd_to_primitive_semirelati_noe')
1942 {
do ix^db=ixomin^db,ixomax^db\}
1943 b2=(^
c&w(ix^
d,
b^
c_)**2+)
1944 if(b2>smalldouble)
then
1952 inv_rho=1.d0/w(ix^
d,
rho_)
1954 b2=b2*inv_rho*inv_squared_c
1956 gamma2=1.d0/(1.d0+b2)
1958 ^
c&w(ix^
d,
m^
c_)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
1961 end subroutine mhd_to_primitive_semirelati_noe
1966 integer,
intent(in) :: ixi^
l, ixo^
l
1967 double precision,
intent(inout) :: w(ixi^s, nw)
1968 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
1973 {
do ix^db=ixomin^db,ixomax^db\}
1976 +half*((^
c&w(ix^
d,
m^
c_)**2+)/&
1978 +(^
c&w(ix^
d,
b^
c_)**2+))
1981 {
do ix^db=ixomin^db,ixomax^db\}
1983 w(ix^d,
e_)=w(ix^d,
e_)&
1984 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
1985 +(^
c&w(ix^d,
b^
c_)**2+))
1992 subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
1994 integer,
intent(in) :: ixi^
l, ixo^
l
1995 double precision,
intent(inout) :: w(ixi^s, nw)
1996 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2000 {
do ix^db=ixomin^db,ixomax^db\}
2006 end subroutine mhd_ei_to_e_hde
2009 subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
2011 integer,
intent(in) :: ixi^
l, ixo^
l
2012 double precision,
intent(inout) :: w(ixi^s, nw)
2013 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2015 w(ixo^s,
p_)=w(ixo^s,
e_)*gamma_1
2016 call mhd_to_conserved_semirelati(ixi^
l,ixo^
l,w,x)
2018 end subroutine mhd_ei_to_e_semirelati
2023 integer,
intent(in) :: ixi^
l, ixo^
l
2024 double precision,
intent(inout) :: w(ixi^s, nw)
2025 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2030 {
do ix^db=ixomin^db,ixomax^db\}
2033 -half*((^
c&w(ix^
d,
m^
c_)**2+)/&
2035 +(^
c&w(ix^
d,
b^
c_)**2+))
2038 {
do ix^db=ixomin^db,ixomax^db\}
2040 w(ix^d,
e_)=w(ix^d,
e_)&
2041 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)&
2042 +(^
c&w(ix^d,
b^
c_)**2+))
2046 if(fix_small_values)
then
2047 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei')
2053 subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
2055 integer,
intent(in) :: ixi^
l, ixo^
l
2056 double precision,
intent(inout) :: w(ixi^s, nw)
2057 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2061 {
do ix^db=ixomin^db,ixomax^db\}
2067 if(fix_small_values)
then
2068 call mhd_handle_small_ei(w,x,ixi^l,ixi^l,
e_,
'mhd_e_to_ei_hde')
2071 end subroutine mhd_e_to_ei_hde
2074 subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
2076 integer,
intent(in) :: ixi^
l, ixo^
l
2077 double precision,
intent(inout) :: w(ixi^s, nw)
2078 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
2080 call mhd_to_primitive_semirelati(ixi^
l,ixo^
l,w,x)
2081 w(ixo^s,
e_)=w(ixo^s,
p_)*inv_gamma_1
2083 end subroutine mhd_e_to_ei_semirelati
2085 subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
2088 logical,
intent(in) :: primitive
2089 integer,
intent(in) :: ixi^
l,ixo^
l
2090 double precision,
intent(inout) :: w(ixi^s,1:nw)
2091 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2092 character(len=*),
intent(in) :: subname
2094 double precision ::
b(ixi^s,1:
ndir), pressure(ixi^s), v(ixi^s,1:
ndir)
2095 double precision :: tmp, b2, gamma2, inv_rho
2097 logical :: flag(ixi^s,1:nw)
2106 {
do ix^db=ixomin^db,ixomax^db\}
2107 b2=(^
c&w(ix^
d,
b^
c_)**2+)
2108 if(b2>smalldouble)
then
2115 inv_rho=1.d0/w(ix^
d,
rho_)
2117 b2=b2*inv_rho*inv_squared_c
2119 gamma2=1.d0/(1.d0+b2)
2121 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
2124 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
2125 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
2126 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2130 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
2136 pressure(ix^
d)=gamma_1*(w(ix^
d,
e_)&
2137 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
2138 +(^
c&w(ix^
d,
b^
c_)**2+)&
2139 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
2146 select case (small_values_method)
2148 {
do ix^db=ixomin^db,ixomax^db\}
2149 if(flag(ix^d,
rho_))
then
2150 w(ix^d,
rho_) = small_density
2151 ^
c&w(ix^d,
m^
c_)=0.d0\
2155 if(flag(ix^d,
e_)) w(ix^d,
p_) = small_pressure
2157 if(flag(ix^d,
e_))
then
2158 w(ix^d,
e_)=small_pressure*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)
2166 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2169 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2171 w(ixo^s,
e_)=pressure(ixo^s)
2172 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2173 {
do ix^db=ixomin^db,ixomax^db\}
2174 w(ix^d,
e_)=w(ix^d,
p_)*inv_gamma_1+half*((^
c&v(ix^d,^
c)**2+)*w(ix^d,
rho_)&
2175 +(^
c&w(ix^d,
b^
c_)**2+)+(^
c&
b(ix^d,^
c)**2+)*inv_squared_c)
2180 if(.not.primitive)
then
2182 w(ixo^s,
mom(1:ndir))=v(ixo^s,1:ndir)
2183 w(ixo^s,
e_)=pressure(ixo^s)
2185 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2189 end subroutine mhd_handle_small_values_semirelati
2191 subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2194 logical,
intent(in) :: primitive
2195 integer,
intent(in) :: ixi^
l,ixo^
l
2196 double precision,
intent(inout) :: w(ixi^s,1:nw)
2197 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2198 character(len=*),
intent(in) :: subname
2201 logical :: flag(ixi^s,1:nw)
2203 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2208 {
do ix^db=ixomin^db,ixomax^db\}
2212 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2224 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2226 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2229 {
do ix^db=iximin^db,iximax^db\}
2230 w(ix^d,
e_)=w(ix^d,
e_)&
2231 -half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2233 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2235 {
do ix^db=iximin^db,iximax^db\}
2236 w(ix^d,
e_)=w(ix^d,
e_)&
2237 +half*((^
c&w(ix^d,
m^
c_)**2+)/w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+))
2241 if(.not.primitive)
then
2243 {
do ix^db=ixomin^db,ixomax^db\}
2245 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2246 -half*((^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_)+(^
c&w(ix^d,
b^
c_)**2+)))
2249 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2253 end subroutine mhd_handle_small_values_origin
2255 subroutine mhd_handle_small_values_split(primitive, w, x, ixI^L, ixO^L, subname)
2258 logical,
intent(in) :: primitive
2259 integer,
intent(in) :: ixi^
l,ixo^
l
2260 double precision,
intent(inout) :: w(ixi^s,1:nw)
2261 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2262 character(len=*),
intent(in) :: subname
2264 double precision :: rho
2266 logical :: flag(ixi^s,1:nw)
2268 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2273 {
do ix^db=ixomin^db,ixomax^db\}
2278 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2285 w(ix^
d,
e_)=small_e+half*((^
c&w(ix^
d,
m^
c_)**2+)/rho+(^
c&w(ix^
d,
b^
c_)**2+))&
2291 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2293 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2296 {
do ix^db=iximin^db,iximax^db\}
2298 w(ix^d,
e_)=w(ix^d,
e_)&
2299 -half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2301 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2303 {
do ix^db=iximin^db,iximax^db\}
2305 w(ix^d,
e_)=w(ix^d,
e_)&
2306 +half*((^
c&w(ix^d,
m^
c_)**2+)/rho+(^
c&w(ix^d,
b^
c_)**2+))
2310 if(.not.primitive)
then
2312 {
do ix^db=ixomin^db,ixomax^db\}
2314 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)/rho\
2315 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)&
2316 -half*((^
c&w(ix^d,
m^
c_)**2+)*rho+(^
c&w(ix^d,
b^
c_)**2+)))
2319 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2323 end subroutine mhd_handle_small_values_split
2325 subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2328 logical,
intent(in) :: primitive
2329 integer,
intent(in) :: ixi^
l,ixo^
l
2330 double precision,
intent(inout) :: w(ixi^s,1:nw)
2331 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2332 character(len=*),
intent(in) :: subname
2335 logical :: flag(ixi^s,1:nw)
2337 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2342 {
do ix^db=ixomin^db,ixomax^db\}
2343 if(flag(ix^
d,
rho_))
then
2345 ^
c&w(ix^
d,
m^
c_)=0.d0\
2350 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e
2355 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2357 call small_values_average(ixi^l, ixo^l, w, x, flag,
p_)
2359 if(.not.primitive)
then
2361 {
do ix^db=ixomin^db,ixomax^db\}
2363 w(ix^d,
p_)=gamma_1*w(ix^d,
e_)
2366 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2370 end subroutine mhd_handle_small_values_inte
2372 subroutine mhd_handle_small_values_noe(primitive, w, x, ixI^L, ixO^L, subname)
2375 logical,
intent(in) :: primitive
2376 integer,
intent(in) :: ixi^
l,ixo^
l
2377 double precision,
intent(inout) :: w(ixi^s,1:nw)
2378 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2379 character(len=*),
intent(in) :: subname
2382 logical :: flag(ixi^s,1:nw)
2384 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2389 {
do ix^db=ixomin^db,ixomax^db\}
2393 if(flag({ix^
d},
rho_)) w({ix^
d},
m^
c_)=0.0d0
2399 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2401 if(.not.primitive)
then
2403 {
do ix^db=ixomin^db,ixomax^db\}
2407 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2411 end subroutine mhd_handle_small_values_noe
2413 subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2416 logical,
intent(in) :: primitive
2417 integer,
intent(in) :: ixi^
l,ixo^
l
2418 double precision,
intent(inout) :: w(ixi^s,1:nw)
2419 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2420 character(len=*),
intent(in) :: subname
2423 logical :: flag(ixi^s,1:nw)
2425 call phys_check_w(primitive, ixi^
l, ixo^
l, w, flag)
2430 {
do ix^db=ixomin^db,ixomax^db\}
2431 if(flag(ix^
d,
rho_))
then
2433 ^
c&w(ix^
d,
m^
c_)=0.d0\
2438 if(flag(ix^
d,
e_)) w(ix^
d,
e_)=small_e+half*(^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)
2443 call small_values_average(ixi^l, ixo^l, w, x, flag,
rho_)
2445 call small_values_average(ixi^l, ixo^l, w, x, flag,
e_)
2447 if(.not.primitive)
then
2449 {
do ix^db=ixomin^db,ixomax^db\}
2451 w(ix^d,
p_)=gamma_1*(w(ix^d,
e_)-half*(^
c&w(ix^d,
m^
c_)**2+)*w(ix^d,
rho_))
2454 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2458 end subroutine mhd_handle_small_values_hde
2464 integer,
intent(in) :: ixi^
l, ixo^
l
2465 double precision,
intent(in) :: w(ixi^s,nw), x(ixi^s,1:
ndim)
2466 double precision,
intent(out) :: v(ixi^s,
ndir)
2468 double precision :: rho(ixi^s)
2473 rho(ixo^s)=1.d0/rho(ixo^s)
2476 v(ixo^s, idir) = w(ixo^s,
mom(idir))*rho(ixo^s)
2482 subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2485 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2486 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2487 double precision,
intent(inout) :: cmax(ixi^s)
2489 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2495 {
do ix^db=ixomin^db,ixomax^db \}
2510 cfast2=b2*inv_rho+cmax(ix^
d)
2511 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2512 if(avmincs2<zero) avmincs2=zero
2513 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2517 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2519 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2522 {
do ix^db=ixomin^db,ixomax^db \}
2536 b2=(^
c&w(ix^d,
b^
c_)**2+)
2537 cfast2=b2*inv_rho+cmax(ix^d)
2538 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2539 if(avmincs2<zero) avmincs2=zero
2540 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2544 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2546 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2550 end subroutine mhd_get_cmax_origin
2553 subroutine mhd_get_cmax_origin_noe(w,x,ixI^L,ixO^L,idim,cmax)
2556 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2557 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2558 double precision,
intent(inout) :: cmax(ixi^s)
2560 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
2566 {
do ix^db=ixomin^db,ixomax^db \}
2577 cfast2=b2*inv_rho+cmax(ix^
d)
2578 avmincs2=cfast2**2-4.0d0*cmax(ix^
d)*(w(ix^
d,mag(idim))+
block%B0(ix^
d,idim,
b0i))**2*inv_rho
2579 if(avmincs2<zero) avmincs2=zero
2580 cmax(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2584 cmax(ix^
d)=max(cmax(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2586 cmax(ix^
d)=abs(w(ix^
d,
mom(idim)))+cmax(ix^
d)
2589 {
do ix^db=ixomin^db,ixomax^db \}
2599 b2=(^
c&w(ix^d,
b^
c_)**2+)
2600 cfast2=b2*inv_rho+cmax(ix^d)
2601 avmincs2=cfast2**2-4.0d0*cmax(ix^d)*w(ix^d,mag(idim))**2*inv_rho
2602 if(avmincs2<zero) avmincs2=zero
2603 cmax(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
2607 cmax(ix^d)=max(cmax(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
2609 cmax(ix^d)=abs(w(ix^d,
mom(idim)))+cmax(ix^d)
2613 end subroutine mhd_get_cmax_origin_noe
2616 subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2619 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2620 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2621 double precision,
intent(inout):: cmax(ixi^s)
2623 double precision :: csound, avmincs2, idim_alfven_speed2
2624 double precision :: inv_rho, alfven_speed2, gamma2
2627 {
do ix^db=ixomin^db,ixomax^db \}
2628 inv_rho=1.d0/w(ix^
d,
rho_)
2629 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2630 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2631 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2634 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2637 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2638 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2639 if(avmincs2<zero) avmincs2=zero
2641 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2642 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2645 end subroutine mhd_get_cmax_semirelati
2648 subroutine mhd_get_cmax_semirelati_noe(w,x,ixI^L,ixO^L,idim,cmax)
2651 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2652 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2653 double precision,
intent(inout):: cmax(ixi^s)
2655 double precision :: csound, avmincs2, idim_alfven_speed2
2656 double precision :: inv_rho, alfven_speed2, gamma2
2659 {
do ix^db=ixomin^db,ixomax^db \}
2660 inv_rho=1.d0/w(ix^
d,
rho_)
2661 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
2662 gamma2=1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2663 cmax(ix^
d)=1.d0-gamma2*w(ix^
d,
mom(idim))**2*inv_squared_c
2665 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
2668 alfven_speed2=alfven_speed2*cmax(ix^
d)+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2669 avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ix^
d)
2670 if(avmincs2<zero) avmincs2=zero
2672 csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2673 cmax(ix^
d)=gamma2*abs(w(ix^
d,
mom(idim)))+csound
2676 end subroutine mhd_get_cmax_semirelati_noe
2678 subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2681 integer,
intent(in) :: ixi^
l, ixo^
l
2682 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
2683 double precision,
intent(inout) :: a2max(
ndim)
2684 double precision :: a2(ixi^s,
ndim,nw)
2685 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i,j
2690 hxo^
l=ixo^
l-
kr(i,^
d);
2691 gxo^
l=hxo^
l-
kr(i,^
d);
2692 jxo^
l=ixo^
l+
kr(i,^
d);
2693 kxo^
l=jxo^
l+
kr(i,^
d);
2694 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2695 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2696 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/
dxlevel(i)**2
2698 end subroutine mhd_get_a2max
2701 subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2704 integer,
intent(in) :: ixi^
l,ixo^
l
2705 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2707 double precision,
intent(inout) :: w(ixi^s,1:nw)
2708 double precision,
intent(out) :: tco_local,tmax_local
2710 double precision,
parameter :: trac_delta=0.25d0
2711 double precision :: te(ixi^s),lts(ixi^s)
2712 double precision,
dimension(1:ndim) :: bdir, bunitvec
2713 double precision,
dimension(ixI^S,1:ndim) :: gradt
2714 double precision :: ltrc,ltrp,altr
2715 integer :: idims,ix^
d,jxo^
l,hxo^
l,ixa^
d,ixb^
d
2716 integer :: jxp^
l,hxp^
l,ixp^
l,ixq^
l
2719 call mhd_get_temperature_from_te(w,x,ixi^
l,ixi^
l,te)
2721 call mhd_get_rfactor(w,x,ixi^
l,ixi^
l,te)
2722 te(ixi^s)=w(ixi^s,
p_)/(te(ixi^s)*w(ixi^s,
rho_))
2725 tmax_local=maxval(te(ixo^s))
2733 do ix1=ixomin1,ixomax1
2734 lts(ix1)=0.5d0*abs(te(ix1+1)-te(ix1-1))/te(ix1)
2735 if(lts(ix1)>trac_delta)
then
2736 tco_local=max(tco_local,te(ix1))
2748 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2749 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2750 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2751 block%wextra(ixo^s,
tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2753 call mpistop(
"mhd_trac_type not allowed for 1D simulation")
2764 call gradient(te,ixi^
l,ixo^
l,idims,gradt(ixi^s,idims))
2771 ixb^
d=(ixomin^
d+ixomax^
d-1)/2+ixa^
d;
2776 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2777 bdir(1:ndim)=bdir(1:ndim)+w(ixb^d,iw_mag(1:ndim))
2781 if(bdir(1)/=0.d0)
then
2782 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2784 block%special_values(3)=0.d0
2786 if(bdir(2)/=0.d0)
then
2787 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2789 block%special_values(4)=0.d0
2793 if(bdir(1)/=0.d0)
then
2794 block%special_values(3)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+&
2795 (bdir(3)/bdir(1))**2)
2797 block%special_values(3)=0.d0
2799 if(bdir(2)/=0.d0)
then
2800 block%special_values(4)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+&
2801 (bdir(3)/bdir(2))**2)
2803 block%special_values(4)=0.d0
2805 if(bdir(3)/=0.d0)
then
2806 block%special_values(5)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+&
2807 (bdir(2)/bdir(3))**2)
2809 block%special_values(5)=0.d0
2814 block%special_values(1)=zero
2815 {
do ix^db=ixomin^db,ixomax^db\}
2817 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2819 ^d&bdir(^d)=w({ix^d},iw_mag(^d))\
2822 if(bdir(1)/=0.d0)
then
2823 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2827 if(bdir(2)/=0.d0)
then
2828 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2833 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2))*&
2834 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2837 if(bdir(1)/=0.d0)
then
2838 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2842 if(bdir(2)/=0.d0)
then
2843 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2847 if(bdir(3)/=0.d0)
then
2848 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2853 lts(ix^d)=min(block%ds(ix^d,1),block%ds(ix^d,2),block%ds(ix^d,3))*&
2854 abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2856 if(lts(ix^d)>trac_delta)
then
2857 block%special_values(1)=max(block%special_values(1),te(ix^d))
2860 block%special_values(2)=tmax_local
2879 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2880 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
2881 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
2885 {
do ix^db=ixpmin^db,ixpmax^db\}
2886 ^d&bdir(^d)=w({ix^d},iw_mag(^d))+block%B0({ix^d},^d,0)\
2888 if(bdir(1)/=0.d0)
then
2889 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
2893 if(bdir(2)/=0.d0)
then
2894 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
2900 if(bdir(1)/=0.d0)
then
2901 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
2905 if(bdir(2)/=0.d0)
then
2906 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
2910 if(bdir(3)/=0.d0)
then
2911 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
2917 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2919 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2920 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2923 {
do ix^db=ixpmin^db,ixpmax^db\}
2925 if(w(ix^d,iw_mag(1))/=0.d0)
then
2926 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)
2930 if(w(ix^d,iw_mag(2))/=0.d0)
then
2931 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)
2937 if(w(ix^d,iw_mag(1))/=0.d0)
then
2938 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+&
2939 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(1)))**2)
2943 if(w(ix^d,iw_mag(2))/=0.d0)
then
2944 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+&
2945 (w(ix^d,iw_mag(3))/w(ix^d,iw_mag(2)))**2)
2949 if(w(ix^d,iw_mag(3))/=0.d0)
then
2950 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+&
2951 (w(ix^d,iw_mag(2))/w(ix^d,iw_mag(3)))**2)
2957 lts(ix^d)=abs(^d&gradt({ix^d},^d)*bunitvec(^d)+)/te(ix^d)
2959 lts(ix^d)=min(^d&block%ds({ix^d},^d))*lts(ix^d)
2960 lts(ix^d)=max(one,(exp(lts(ix^d))/ltrc)**ltrp)
2966 {
do ix^db=ixpmin^db,ixpmax^db\}
2968 altr=0.25d0*((lts(ix1-1,ix2)+two*lts(ix^d)+lts(ix1+1,ix2))*bunitvec(1)**2+&
2969 (lts(ix1,ix2-1)+two*lts(ix^d)+lts(ix1,ix2+1))*bunitvec(2)**2)
2970 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2973 altr=0.25d0*((lts(ix1-1,ix2,ix3)+two*lts(ix^d)+lts(ix1+1,ix2,ix3))*bunitvec(1)**2+&
2974 (lts(ix1,ix2-1,ix3)+two*lts(ix^d)+lts(ix1,ix2+1,ix3))*bunitvec(2)**2+&
2975 (lts(ix1,ix2,ix3-1)+two*lts(ix^d)+lts(ix1,ix2,ix3+1))*bunitvec(3)**2)
2976 block%wextra(ix^d,
tcoff_)=te(ix^d)*altr**0.4d0
2982 call mpistop(
"unknown mhd_trac_type")
2985 end subroutine mhd_get_tcutoff
2988 subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2991 integer,
intent(in) :: ixi^
l, ixo^
l, idim
2992 double precision,
intent(in) :: wprim(ixi^s, nw)
2993 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2994 double precision,
intent(out) :: hspeed(ixi^s,1:number_species)
2996 double precision :: csound(ixi^s,
ndim)
2997 double precision,
allocatable :: tmp(:^
d&)
2998 integer :: jxc^
l, ixc^
l, ixa^
l, id, ix^
d
3002 allocate(tmp(ixa^s))
3004 call mhd_get_csound_prim(wprim,x,ixi^
l,ixa^
l,id,tmp)
3005 csound(ixa^s,id)=tmp(ixa^s)
3008 ixcmin^
d=ixomin^
d+
kr(idim,^
d)-1;
3009 jxcmax^
d=ixcmax^
d+
kr(idim,^
d);
3010 jxcmin^
d=ixcmin^
d+
kr(idim,^
d);
3011 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))
3015 ixamax^
d=ixcmax^
d+
kr(id,^
d);
3016 ixamin^
d=ixcmin^
d+
kr(id,^
d);
3017 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)))
3018 ixamax^
d=ixcmax^
d-
kr(id,^
d);
3019 ixamin^
d=ixcmin^
d-
kr(id,^
d);
3020 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)))
3025 ixamax^
d=jxcmax^
d+
kr(id,^
d);
3026 ixamin^
d=jxcmin^
d+
kr(id,^
d);
3027 hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,
mom(id))+csound(ixa^s,id)-wprim(jxc^s,
mom(id))+csound(jxc^s,id)))
3028 ixamax^
d=jxcmax^
d-
kr(id,^
d);
3029 ixamin^
d=jxcmin^
d-
kr(id,^
d);
3030 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)))
3034 end subroutine mhd_get_h_speed
3037 subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3040 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3041 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3042 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3043 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3044 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3045 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3046 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3048 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3049 double precision :: umean, dmean, tmp1, tmp2, tmp3
3056 call mhd_get_csound_prim(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3057 call mhd_get_csound_prim(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3058 if(
present(cmin))
then
3059 {
do ix^db=ixomin^db,ixomax^db\}
3060 tmp1=sqrt(wlp(ix^
d,
rho_))
3061 tmp2=sqrt(wrp(ix^
d,
rho_))
3062 tmp3=1.d0/(tmp1+tmp2)
3063 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3064 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3065 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3066 cmin(ix^
d,1)=umean-dmean
3067 cmax(ix^
d,1)=umean+dmean
3069 if(h_correction)
then
3070 {
do ix^db=ixomin^db,ixomax^db\}
3071 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3072 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3076 {
do ix^db=ixomin^db,ixomax^db\}
3077 tmp1=sqrt(wlp(ix^d,
rho_))
3078 tmp2=sqrt(wrp(ix^d,
rho_))
3079 tmp3=1.d0/(tmp1+tmp2)
3080 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3081 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3082 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3083 cmax(ix^d,1)=abs(umean)+dmean
3087 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3088 call mhd_get_csound_prim(wmean,x,ixi^l,ixo^l,idim,csoundr)
3089 if(
present(cmin))
then
3090 {
do ix^db=ixomin^db,ixomax^db\}
3091 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3092 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3094 if(h_correction)
then
3095 {
do ix^db=ixomin^db,ixomax^db\}
3096 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3097 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3101 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3105 call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
3106 call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
3107 if(
present(cmin))
then
3108 {
do ix^db=ixomin^db,ixomax^db\}
3109 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3110 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3111 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3113 if(h_correction)
then
3114 {
do ix^db=ixomin^db,ixomax^db\}
3115 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3116 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3120 {
do ix^db=ixomin^db,ixomax^db\}
3121 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3122 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3127 end subroutine mhd_get_cbounds
3130 subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3133 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3134 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3135 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3136 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3137 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3138 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3139 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3141 double precision,
dimension(ixO^S) :: csoundl, csoundr, gamma2l, gamma2r
3146 call mhd_get_csound_semirelati(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3147 call mhd_get_csound_semirelati(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3149 call mhd_get_csound_semirelati_noe(wlp,x,ixi^
l,ixo^
l,idim,csoundl,gamma2l)
3150 call mhd_get_csound_semirelati_noe(wrp,x,ixi^
l,ixo^
l,idim,csoundr,gamma2r)
3152 if(
present(cmin))
then
3153 {
do ix^db=ixomin^db,ixomax^db\}
3154 csoundl(ix^
d)=max(csoundl(ix^
d),csoundr(ix^
d))
3155 cmin(ix^
d,1)=min(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))-csoundl(ix^
d)
3156 cmax(ix^
d,1)=max(gamma2l(ix^
d)*wlp(ix^
d,
mom(idim)),gamma2r(ix^
d)*wrp(ix^
d,
mom(idim)))+csoundl(ix^
d)
3159 {
do ix^db=ixomin^db,ixomax^db\}
3160 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3161 cmax(ix^d,1)=max(gamma2l(ix^d)*wlp(ix^d,
mom(idim)),gamma2r(ix^d)*wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3165 end subroutine mhd_get_cbounds_semirelati
3168 subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
3171 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3172 double precision,
intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
3173 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3174 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3175 double precision,
intent(inout) :: cmax(ixi^s,1:number_species)
3176 double precision,
intent(inout),
optional :: cmin(ixi^s,1:number_species)
3177 double precision,
intent(in) :: hspeed(ixi^s,1:number_species)
3179 double precision :: wmean(ixi^s,nw), csoundl(ixo^s), csoundr(ixo^s)
3180 double precision :: umean, dmean, tmp1, tmp2, tmp3
3187 call mhd_get_csound_prim_split(wlp,x,ixi^
l,ixo^
l,idim,csoundl)
3188 call mhd_get_csound_prim_split(wrp,x,ixi^
l,ixo^
l,idim,csoundr)
3189 if(
present(cmin))
then
3190 {
do ix^db=ixomin^db,ixomax^db\}
3193 tmp3=1.d0/(tmp1+tmp2)
3194 umean=(wlp(ix^
d,
mom(idim))*tmp1+wrp(ix^
d,
mom(idim))*tmp2)*tmp3
3195 dmean=sqrt((tmp1*csoundl(ix^
d)**2+tmp2*csoundr(ix^
d)**2)*tmp3+&
3196 half*tmp1*tmp2*tmp3**2*(wrp(ix^
d,
mom(idim))-wlp(ix^
d,
mom(idim)))**2)
3197 cmin(ix^
d,1)=umean-dmean
3198 cmax(ix^
d,1)=umean+dmean
3200 if(h_correction)
then
3201 {
do ix^db=ixomin^db,ixomax^db\}
3202 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3203 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3207 {
do ix^db=ixomin^db,ixomax^db\}
3210 tmp3=1.d0/(tmp1+tmp2)
3211 umean=(wlp(ix^d,
mom(idim))*tmp1+wrp(ix^d,
mom(idim))*tmp2)*tmp3
3212 dmean=sqrt((tmp1*csoundl(ix^d)**2+tmp2*csoundr(ix^d)**2)*tmp3+&
3213 half*tmp1*tmp2*tmp3**2*(wrp(ix^d,
mom(idim))-wlp(ix^d,
mom(idim)))**2)
3214 cmax(ix^d,1)=abs(umean)+dmean
3218 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
3219 call mhd_get_csound_prim_split(wmean,x,ixi^l,ixo^l,idim,csoundr)
3220 if(
present(cmin))
then
3221 {
do ix^db=ixomin^db,ixomax^db\}
3222 cmax(ix^d,1)=max(wmean(ix^d,
mom(idim))+csoundr(ix^d),zero)
3223 cmin(ix^d,1)=min(wmean(ix^d,
mom(idim))-csoundr(ix^d),zero)
3225 if(h_correction)
then
3226 {
do ix^db=ixomin^db,ixomax^db\}
3227 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3228 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3232 cmax(ixo^s,1)=abs(wmean(ixo^s,
mom(idim)))+csoundr(ixo^s)
3236 call mhd_get_csound_prim_split(wlp,x,ixi^l,ixo^l,idim,csoundl)
3237 call mhd_get_csound_prim_split(wrp,x,ixi^l,ixo^l,idim,csoundr)
3238 if(
present(cmin))
then
3239 {
do ix^db=ixomin^db,ixomax^db\}
3240 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3241 cmin(ix^d,1)=min(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))-csoundl(ix^d)
3242 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3244 if(h_correction)
then
3245 {
do ix^db=ixomin^db,ixomax^db\}
3246 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
3247 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
3251 {
do ix^db=ixomin^db,ixomax^db\}
3252 csoundl(ix^d)=max(csoundl(ix^d),csoundr(ix^d))
3253 cmax(ix^d,1)=max(wlp(ix^d,
mom(idim)),wrp(ix^d,
mom(idim)))+csoundl(ix^d)
3258 end subroutine mhd_get_cbounds_split_rho
3261 subroutine mhd_get_ct_velocity_average(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3264 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3265 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3266 double precision,
intent(in) :: cmax(ixi^s)
3267 double precision,
intent(in),
optional :: cmin(ixi^s)
3268 type(ct_velocity),
intent(inout):: vcts
3270 end subroutine mhd_get_ct_velocity_average
3272 subroutine mhd_get_ct_velocity_contact(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3275 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3276 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3277 double precision,
intent(in) :: cmax(ixi^s)
3278 double precision,
intent(in),
optional :: cmin(ixi^s)
3279 type(ct_velocity),
intent(inout):: vcts
3281 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
3283 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
3285 end subroutine mhd_get_ct_velocity_contact
3287 subroutine mhd_get_ct_velocity_hll(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
3290 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3291 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
3292 double precision,
intent(in) :: cmax(ixi^s)
3293 double precision,
intent(in),
optional :: cmin(ixi^s)
3294 type(ct_velocity),
intent(inout):: vcts
3296 integer :: idime,idimn
3298 if(.not.
allocated(vcts%vbarC))
then
3299 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
3300 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
3303 if(
present(cmin))
then
3304 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
3305 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3307 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
3308 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
3311 idimn=mod(idim,
ndir)+1
3312 idime=mod(idim+1,
ndir)+1
3314 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
3315 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
3316 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
3317 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3318 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3320 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
3321 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
3322 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
3323 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
3324 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
3326 end subroutine mhd_get_ct_velocity_hll
3329 subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3332 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3333 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3334 double precision,
intent(out):: csound(ixo^s)
3336 double precision :: inv_rho, cfast2, avmincs2, b2, kmax
3343 {
do ix^db=ixomin^db,ixomax^db \}
3344 inv_rho=1.d0/w(ix^
d,
rho_)
3351 cfast2=b2*inv_rho+csound(ix^
d)
3352 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3354 if(avmincs2<zero) avmincs2=zero
3355 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3357 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3361 {
do ix^db=ixomin^db,ixomax^db \}
3362 inv_rho=1.d0/w(ix^d,
rho_)
3368 b2=(^
c&w(ix^d,
b^
c_)**2+)
3369 cfast2=b2*inv_rho+csound(ix^d)
3370 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3371 if(avmincs2<zero) avmincs2=zero
3372 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3374 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3379 end subroutine mhd_get_csound_prim
3382 subroutine mhd_get_csound_prim_split(w,x,ixI^L,ixO^L,idim,csound)
3385 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3386 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3387 double precision,
intent(out):: csound(ixo^s)
3389 double precision :: rho, inv_rho, cfast2, avmincs2, b2, kmax
3396 {
do ix^db=ixomin^db,ixomax^db \}
3403 cfast2=b2*inv_rho+csound(ix^
d)
3404 avmincs2=cfast2**2-4.0d0*csound(ix^
d)*(w(ix^
d,mag(idim))+&
3406 if(avmincs2<zero) avmincs2=zero
3407 csound(ix^
d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3409 csound(ix^
d)=max(csound(ix^
d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3413 {
do ix^db=ixomin^db,ixomax^db \}
3419 b2=(^
c&w(ix^d,
b^
c_)**2+)
3420 cfast2=b2*inv_rho+csound(ix^d)
3421 avmincs2=cfast2**2-4.0d0*csound(ix^d)*w(ix^d,mag(idim))**2*inv_rho
3422 if(avmincs2<zero) avmincs2=zero
3423 csound(ix^d)=sqrt(half*(cfast2+sqrt(avmincs2)))
3425 csound(ix^d)=max(csound(ix^d),
mhd_etah*sqrt(b2)*inv_rho*kmax)
3430 end subroutine mhd_get_csound_prim_split
3433 subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3436 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3438 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3439 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3441 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3444 {
do ix^db=ixomin^db,ixomax^db\}
3445 inv_rho = 1.d0/w(ix^
d,
rho_)
3448 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3449 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3450 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3451 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3454 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3455 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3456 if(avmincs2<zero) avmincs2=zero
3458 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3461 end subroutine mhd_get_csound_semirelati
3464 subroutine mhd_get_csound_semirelati_noe(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3467 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3469 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
3470 double precision,
intent(out):: csound(ixo^s), gamma2(ixo^s)
3472 double precision :: avmincs2, inv_rho, alfven_speed2, idim_alfven_speed2
3475 {
do ix^db=ixomin^db,ixomax^db\}
3476 inv_rho = 1.d0/w(ix^
d,
rho_)
3479 alfven_speed2=(^
c&w(ix^
d,
b^
c_)**2+)*inv_rho
3480 gamma2(ix^
d) = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3481 avmincs2=1.d0-gamma2(ix^
d)*w(ix^
d,
mom(idim))**2*inv_squared_c
3482 idim_alfven_speed2=w(ix^
d,mag(idim))**2*inv_rho
3485 alfven_speed2=alfven_speed2*avmincs2+csound(ix^
d)*(1.d0+idim_alfven_speed2*inv_squared_c)
3486 avmincs2=(gamma2(ix^
d)*alfven_speed2)**2-4.0d0*gamma2(ix^
d)*csound(ix^
d)*idim_alfven_speed2*avmincs2
3487 if(avmincs2<zero) avmincs2=zero
3489 csound(ix^
d) = sqrt(half*(gamma2(ix^
d)*alfven_speed2+sqrt(avmincs2)))
3492 end subroutine mhd_get_csound_semirelati_noe
3495 subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
3498 integer,
intent(in) :: ixi^
l, ixo^
l
3499 double precision,
intent(in) :: w(ixi^s,nw)
3500 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3501 double precision,
intent(out):: pth(ixi^s)
3509 end subroutine mhd_get_pthermal_noe
3512 subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
3516 integer,
intent(in) :: ixi^
l, ixo^
l
3517 double precision,
intent(in) :: w(ixi^s,nw)
3518 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3519 double precision,
intent(out):: pth(ixi^s)
3523 {
do ix^db= ixomin^db,ixomax^db\}
3527 pth(ix^
d)=gamma_1*w(ix^
d,
e_)
3532 if(check_small_values.and..not.fix_small_values)
then
3533 {
do ix^db= ixomin^db,ixomax^db\}
3534 if(pth(ix^d)<small_pressure)
then
3535 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3536 " encountered when call mhd_get_pthermal_inte"
3537 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3538 write(*,*)
"Location: ", x(ix^d,:)
3539 write(*,*)
"Cell number: ", ix^d
3541 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3544 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3545 write(*,*)
"Saving status at the previous time step"
3551 end subroutine mhd_get_pthermal_inte
3554 subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3558 integer,
intent(in) :: ixi^
l, ixo^
l
3559 double precision,
intent(in) :: w(ixi^s,nw)
3560 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3561 double precision,
intent(out):: pth(ixi^s)
3565 {
do ix^db=ixomin^db,ixomax^db\}
3570 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)&
3571 +(^
c&w(ix^
d,
b^
c_)**2+)))
3576 if(check_small_values.and..not.fix_small_values)
then
3577 {
do ix^db=ixomin^db,ixomax^db\}
3578 if(pth(ix^d)<small_pressure)
then
3579 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3580 " encountered when call mhd_get_pthermal"
3581 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3582 write(*,*)
"Location: ", x(ix^d,:)
3583 write(*,*)
"Cell number: ", ix^d
3585 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3588 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3589 write(*,*)
"Saving status at the previous time step"
3595 end subroutine mhd_get_pthermal_origin
3598 subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3602 integer,
intent(in) :: ixi^
l, ixo^
l
3603 double precision,
intent(in) :: w(ixi^s,nw)
3604 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3605 double precision,
intent(out):: pth(ixi^s)
3607 double precision ::
b(ixo^s,1:
ndir), v(ixo^s,1:
ndir), tmp, b2, gamma2, inv_rho
3610 {
do ix^db=ixomin^db,ixomax^db\}
3611 b2=(^
c&w(ix^
d,
b^
c_)**2+)
3612 if(b2>smalldouble)
then
3620 inv_rho=1.d0/w(ix^
d,
rho_)
3622 b2=b2*inv_rho*inv_squared_c
3624 gamma2=1.d0/(1.d0+b2)
3626 ^
c&v(ix^
d,^
c)=gamma2*(w(ix^
d,
m^
c_)+b2*
b(ix^
d,^
c)*tmp)*inv_rho\
3630 b(ix^
d,1)=w(ix^
d,b2_)*v(ix^
d,3)-w(ix^
d,b3_)*v(ix^
d,2)
3631 b(ix^
d,2)=w(ix^
d,b3_)*v(ix^
d,1)-w(ix^
d,b1_)*v(ix^
d,3)
3632 b(ix^
d,3)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3636 b(ix^
d,2)=w(ix^
d,b1_)*v(ix^
d,2)-w(ix^
d,b2_)*v(ix^
d,1)
3642 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)&
3643 -half*((^
c&v(ix^
d,^
c)**2+)*w(ix^
d,
rho_)&
3644 +(^
c&w(ix^
d,
b^
c_)**2+)&
3645 +(^
c&
b(ix^
d,^
c)**2+)*inv_squared_c))
3649 if(check_small_values.and..not.fix_small_values)
then
3650 {
do ix^db=ixomin^db,ixomax^db\}
3651 if(pth(ix^d)<small_pressure)
then
3652 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3653 " encountered when call mhd_get_pthermal_semirelati"
3654 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3655 write(*,*)
"Location: ", x(ix^d,:)
3656 write(*,*)
"Cell number: ", ix^d
3658 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3661 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3662 write(*,*)
"Saving status at the previous time step"
3668 end subroutine mhd_get_pthermal_semirelati
3671 subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3675 integer,
intent(in) :: ixi^
l, ixo^
l
3676 double precision,
intent(in) :: w(ixi^s,nw)
3677 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3678 double precision,
intent(out):: pth(ixi^s)
3682 {
do ix^db= ixomin^db,ixomax^db\}
3683 pth(ix^
d)=gamma_1*(w(ix^
d,
e_)-half*((^
c&w(ix^
d,
m^
c_)**2+)/w(ix^
d,
rho_)))
3686 if(check_small_values.and..not.fix_small_values)
then
3687 {
do ix^db= ixomin^db,ixomax^db\}
3688 if(pth(ix^d)<small_pressure)
then
3689 write(*,*)
"Error: small value of gas pressure",pth(ix^d),&
3690 " encountered when call mhd_get_pthermal_hde"
3691 write(*,*)
"Iteration: ", it,
" Time: ", global_time
3692 write(*,*)
"Location: ", x(ix^d,:)
3693 write(*,*)
"Cell number: ", ix^d
3695 write(*,*) trim(cons_wnames(iw)),
": ",w(ix^d,iw)
3698 if(trace_small_values)
write(*,*) sqrt(pth(ix^d)-bigdouble)
3699 write(*,*)
"Saving status at the previous time step"
3705 end subroutine mhd_get_pthermal_hde
3708 subroutine mhd_get_temperature_from_te(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)
3714 res(ixo^s) = w(ixo^s,
te_)
3715 end subroutine mhd_get_temperature_from_te
3718 subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3720 integer,
intent(in) :: ixi^
l, ixo^
l
3721 double precision,
intent(in) :: w(ixi^s, 1:nw)
3722 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3723 double precision,
intent(out):: res(ixi^s)
3725 double precision :: r(ixi^s)
3727 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3728 res(ixo^s) = gamma_1 * w(ixo^s,
e_)/(w(ixo^s,
rho_)*r(ixo^s))
3729 end subroutine mhd_get_temperature_from_eint
3732 subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3734 integer,
intent(in) :: ixi^
l, ixo^
l
3735 double precision,
intent(in) :: w(ixi^s, 1:nw)
3736 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3737 double precision,
intent(out):: res(ixi^s)
3739 double precision :: r(ixi^s)
3741 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3743 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,
rho_))
3745 end subroutine mhd_get_temperature_from_etot
3747 subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
3749 integer,
intent(in) :: ixi^
l, ixo^
l
3750 double precision,
intent(in) :: w(ixi^s, 1:nw)
3751 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3752 double precision,
intent(out):: res(ixi^s)
3754 double precision :: r(ixi^s)
3756 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3760 end subroutine mhd_get_temperature_from_etot_with_equi
3762 subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3764 integer,
intent(in) :: ixi^
l, ixo^
l
3765 double precision,
intent(in) :: w(ixi^s, 1:nw)
3766 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3767 double precision,
intent(out):: res(ixi^s)
3769 double precision :: r(ixi^s)
3771 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3775 end subroutine mhd_get_temperature_from_eint_with_equi
3777 subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3779 integer,
intent(in) :: ixi^
l, ixo^
l
3780 double precision,
intent(in) :: w(ixi^s, 1:nw)
3781 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3782 double precision,
intent(out):: res(ixi^s)
3784 double precision :: r(ixi^s)
3786 call mhd_get_rfactor(w,x,ixi^
l,ixo^
l,r)
3789 end subroutine mhd_get_temperature_equi
3791 subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3793 integer,
intent(in) :: ixi^
l, ixo^
l
3794 double precision,
intent(in) :: w(ixi^s, 1:nw)
3795 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3796 double precision,
intent(out):: res(ixi^s)
3798 end subroutine mhd_get_rho_equi
3800 subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3802 integer,
intent(in) :: ixi^
l, ixo^
l
3803 double precision,
intent(in) :: w(ixi^s, 1:nw)
3804 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
3805 double precision,
intent(out):: res(ixi^s)
3807 end subroutine mhd_get_pe_equi
3810 subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3814 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3816 double precision,
intent(in) :: wc(ixi^s,nw)
3818 double precision,
intent(in) :: w(ixi^s,nw)
3819 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3820 double precision,
intent(out) :: f(ixi^s,nwflux)
3822 double precision :: vhall(ixi^s,1:
ndir)
3823 double precision :: ptotal
3827 {
do ix^db=ixomin^db,ixomax^db\}
3840 {
do ix^db=ixomin^db,ixomax^db\}
3844 ^
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_)\
3845 ptotal=w(ix^d,
p_)+half*(^
c&w(ix^d,
b^
c_)**2+)
3847 f(ix^d,
mom(idim))=f(ix^d,
mom(idim))+ptotal
3850 f(ix^d,
e_)=w(ix^d,
mom(idim))*(wc(ix^d,
e_)+ptotal)&
3851 -w(ix^d,mag(idim))*(^
c&w(ix^d,
b^
c_)*w(ix^d,
m^
c_)+)
3853 ^
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_)\
3857 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3858 {
do ix^db=ixomin^db,ixomax^db\}
3859 if(total_energy)
then
3861 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)**2+)&
3862 -w(ix^d,mag(idim))*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
3865 ^
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))\
3869 {
do ix^db=ixomin^db,ixomax^db\}
3870 f(ix^d,mag(idim))=w(ix^d,
psi_)
3872 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3877 {
do ix^db=ixomin^db,ixomax^db\}
3883 {
do ix^db=ixomin^db,ixomax^db\}
3884 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)
3889 end subroutine mhd_get_flux
3892 subroutine mhd_get_flux_noe(wC,w,x,ixI^L,ixO^L,idim,f)
3896 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3898 double precision,
intent(in) :: wc(ixi^s,nw)
3900 double precision,
intent(in) :: w(ixi^s,nw)
3901 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3902 double precision,
intent(out) :: f(ixi^s,nwflux)
3904 double precision :: vhall(ixi^s,1:
ndir)
3907 {
do ix^db=ixomin^db,ixomax^db\}
3918 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3919 {
do ix^db=ixomin^db,ixomax^db\}
3921 ^
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))\
3925 {
do ix^db=ixomin^db,ixomax^db\}
3926 f(ix^d,mag(idim))=w(ix^d,
psi_)
3928 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3933 {
do ix^db=ixomin^db,ixomax^db\}
3938 end subroutine mhd_get_flux_noe
3941 subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3945 integer,
intent(in) :: ixi^
l, ixo^
l, idim
3947 double precision,
intent(in) :: wc(ixi^s,nw)
3949 double precision,
intent(in) :: w(ixi^s,nw)
3950 double precision,
intent(in) :: x(ixi^s,1:
ndim)
3951 double precision,
intent(out) :: f(ixi^s,nwflux)
3953 double precision :: vhall(ixi^s,1:
ndir)
3956 {
do ix^db=ixomin^db,ixomax^db\}
3969 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3970 {
do ix^db=ixomin^db,ixomax^db\}
3972 ^
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))\
3976 {
do ix^db=ixomin^db,ixomax^db\}
3977 f(ix^d,mag(idim))=w(ix^d,
psi_)
3979 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
3984 {
do ix^db=ixomin^db,ixomax^db\}
3990 {
do ix^db=ixomin^db,ixomax^db\}
3991 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)
3996 end subroutine mhd_get_flux_hde
3999 subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
4003 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4005 double precision,
intent(in) :: wc(ixi^s,nw)
4007 double precision,
intent(in) :: w(ixi^s,nw)
4008 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4009 double precision,
intent(out) :: f(ixi^s,nwflux)
4011 double precision :: vhall(ixi^s,1:
ndir)
4012 double precision :: ptotal, btotal(ixo^s,1:
ndir)
4015 {
do ix^db=ixomin^db,ixomax^db\}
4023 ptotal=w(ix^
d,
p_)+half*(^
c&w(ix^
d,
b^
c_)**2+)
4027 ptotal=ptotal+(^
c&w(ix^
d,
b^
c_)*
block%B0(ix^
d,^
c,idim)+)
4031 btotal(ix^
d,idim)*w(ix^
d,
b^
c_)-w(ix^
d,mag(idim))*
block%B0(ix^
d,^
c,idim)\
4032 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4034 ^
c&btotal(ix^
d,^
c)=w(ix^
d,
b^
c_)\
4038 f(ix^
d,
mom(idim))=f(ix^
d,
mom(idim))+ptotal
4041 ^
c&f(ix^
d,
b^
c_)=w(ix^
d,
mom(idim))*btotal(ix^
d,^
c)-btotal(ix^
d,idim)*w(ix^
d,
m^
c_)\
4048 f(ix^
d,
e_)=w(ix^
d,
mom(idim))*(wc(ix^
d,
e_)+ptotal)&
4049 -btotal(ix^
d,idim)*(^
c&w(ix^
d,
b^
c_)*w(ix^
d,
m^
c_)+)
4054 {
do ix^db=ixomin^db,ixomax^db\}
4055 f(ix^d,mag(idim))=w(ix^d,
psi_)
4057 f(ix^d,
psi_) = cmax_global**2*w(ix^d,mag(idim))
4062 call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
4063 {
do ix^db=ixomin^db,ixomax^db\}
4065 ^
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))\
4066 if(total_energy)
then
4068 f(ix^d,
e_)=f(ix^d,
e_)+vhall(ix^d,idim)*(^
c&w(ix^d,
b^
c_)*btotal(ix^d,^
c)+)&
4069 -btotal(ix^d,idim)*(^
c&vhall(ix^d,^
c)*w(ix^d,
b^
c_)+)
4075 {
do ix^db=ixomin^db,ixomax^db\}
4080 {
do ix^db=ixomin^db,ixomax^db\}
4081 f(ix^d,
e_)=f(ix^d,
e_)+w(ix^d,
q_)*btotal(ix^d,idim)/(dsqrt(^
c&btotal(ix^d,^
c)**2+)+smalldouble)
4086 end subroutine mhd_get_flux_split
4089 subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
4093 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4095 double precision,
intent(in) :: wc(ixi^s,nw)
4097 double precision,
intent(in) :: w(ixi^s,nw)
4098 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4099 double precision,
intent(out) :: f(ixi^s,nwflux)
4101 double precision :: sa(ixo^s,1:
ndir),e(ixo^s,1:
ndir),e2
4104 {
do ix^db=ixomin^db,ixomax^db\}
4109 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4110 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4111 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4116 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4121 e2=(^
c&e(ix^
d,^
c)**2+)
4128 sa(ix^
d,1)=e(ix^
d,2)*w(ix^
d,b3_)-e(ix^
d,3)*w(ix^
d,b2_)
4129 sa(ix^
d,2)=e(ix^
d,3)*w(ix^
d,b1_)-e(ix^
d,1)*w(ix^
d,b3_)
4130 sa(ix^
d,3)=e(ix^
d,1)*w(ix^
d,b2_)-e(ix^
d,2)*w(ix^
d,b1_)
4133 sa(ix^
d,1)=-e(ix^
d,2)*w(ix^
d,b2_)
4134 sa(ix^
d,2)=e(ix^
d,2)*w(ix^
d,b1_)
4147 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4149 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)
4156 {
do ix^db=ixomin^db,ixomax^db\}
4157 f(ix^d,mag(idim))=w(ix^d,
psi_)
4159 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4164 {
do ix^db=ixomin^db,ixomax^db\}
4169 {
do ix^db=ixomin^db,ixomax^db\}
4170 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)
4175 end subroutine mhd_get_flux_semirelati
4177 subroutine mhd_get_flux_semirelati_noe(wC,w,x,ixI^L,ixO^L,idim,f)
4181 integer,
intent(in) :: ixi^
l, ixo^
l, idim
4183 double precision,
intent(in) :: wc(ixi^s,nw)
4185 double precision,
intent(in) :: w(ixi^s,nw)
4186 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4187 double precision,
intent(out) :: f(ixi^s,nwflux)
4189 double precision :: e(ixo^s,1:
ndir),e2
4192 {
do ix^db=ixomin^db,ixomax^db\}
4197 e(ix^
d,1)=w(ix^
d,b2_)*w(ix^
d,m3_)-w(ix^
d,b3_)*w(ix^
d,m2_)
4198 e(ix^
d,2)=w(ix^
d,b3_)*w(ix^
d,m1_)-w(ix^
d,b1_)*w(ix^
d,m3_)
4199 e(ix^
d,3)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4200 e2=(^
c&e(ix^
d,^
c)**2+)
4205 e(ix^
d,2)=w(ix^
d,b1_)*w(ix^
d,m2_)-w(ix^
d,b2_)*w(ix^
d,m1_)
4214 -w(ix^
d,mag(idim))*w(ix^
d,
b^
c_)-e(ix^
d,idim)*e(ix^
d,^
c)*inv_squared_c\
4223 {
do ix^db=ixomin^db,ixomax^db\}
4224 f(ix^d,mag(idim))=w(ix^d,
psi_)
4226 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
4231 {
do ix^db=ixomin^db,ixomax^db\}
4236 end subroutine mhd_get_flux_semirelati_noe
4242 subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4244 integer,
intent(in) :: ixi^
l, ixo^
l,ie
4245 double precision,
intent(in) :: qdt
4246 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4247 double precision,
intent(inout) :: w(ixi^s,1:nw)
4248 double precision :: tmp(ixi^s)
4249 double precision :: jxbxb(ixi^s,1:3)
4251 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,jxbxb)
4254 w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
4256 end subroutine add_source_ambipolar_internal_energy
4258 subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
4261 integer,
intent(in) :: ixi^
l, ixo^
l
4262 double precision,
intent(in) :: w(ixi^s,nw)
4263 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4264 double precision,
intent(out) :: res(:^
d&,:)
4266 double precision :: btot(ixi^s,1:3)
4267 double precision :: current(ixi^s,7-2*
ndir:3)
4268 double precision :: tmp(ixi^s),b2(ixi^s)
4269 integer :: idir, idirmin
4278 btot(ixo^s, idir) = w(ixo^s,mag(idir)) +
block%B0(ixo^s,idir,
b0i)
4281 btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4284 tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=
ndim+1)
4285 b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=
ndim+1)
4287 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4290 res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4292 end subroutine mhd_get_jxbxb
4299 subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4303 integer,
intent(in) :: ixi^
l, ixo^
l,igrid,nflux
4304 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4305 double precision,
intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
4306 double precision,
intent(in) :: my_dt
4307 logical,
intent(in) :: fix_conserve_at_step
4309 double precision,
dimension(ixI^S,1:3) :: tmp,ff
4310 double precision :: fluxall(ixi^s,1:nflux,1:
ndim)
4311 double precision :: fe(ixi^s,
sdim:3)
4312 double precision :: btot(ixi^s,1:3),tmp2(ixi^s)
4313 integer :: i, ixa^
l, ie_
4319 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,tmp)
4329 btot(ixa^s,1:3)=0.d0
4335 btot(ixa^s,1:
ndir) = w(ixa^s,mag(1:
ndir))
4338 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4339 if(fix_conserve_at_step) fluxall(ixi^s,1,1:
ndim)=ff(ixi^s,1:
ndim)
4341 wres(ixo^s,
e_)=-tmp2(ixo^s)
4347 ff(ixa^s,1) = tmp(ixa^s,2)
4348 ff(ixa^s,2) = -tmp(ixa^s,1)
4350 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4351 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4352 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4355 call update_faces_ambipolar(ixi^
l,ixo^
l,w,x,tmp,fe,btot)
4357 ixamin^
d=ixomin^
d-1;
4358 wres(ixa^s,mag(1:
ndim))=-btot(ixa^s,1:
ndim)
4367 ff(ixa^s,2) = tmp(ixa^s,3)
4368 ff(ixa^s,3) = -tmp(ixa^s,2)
4369 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4370 if(fix_conserve_at_step) fluxall(ixi^s,2,1:
ndim)=ff(ixi^s,1:
ndim)
4372 wres(ixo^s,mag(1))=-tmp2(ixo^s)
4374 ff(ixa^s,1) = -tmp(ixa^s,3)
4376 ff(ixa^s,3) = tmp(ixa^s,1)
4377 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4378 if(fix_conserve_at_step) fluxall(ixi^s,3,1:
ndim)=ff(ixi^s,1:
ndim)
4379 wres(ixo^s,mag(2))=-tmp2(ixo^s)
4383 ff(ixa^s,1) = tmp(ixa^s,2)
4384 ff(ixa^s,2) = -tmp(ixa^s,1)
4386 call get_flux_on_cell_face(ixi^
l,ixo^
l,ff,tmp2)
4387 if(fix_conserve_at_step) fluxall(ixi^s,1+
ndir,1:
ndim)=ff(ixi^s,1:
ndim)
4388 wres(ixo^s,mag(
ndir))=-tmp2(ixo^s)
4393 if(fix_conserve_at_step)
then
4394 fluxall=my_dt*fluxall
4401 end subroutine sts_set_source_ambipolar
4404 subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4407 integer,
intent(in) :: ixi^
l, ixo^
l
4408 double precision,
intent(in) :: w(ixi^s,1:nw)
4409 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4411 double precision,
intent(in) :: ecc(ixi^s,1:3)
4412 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
4413 double precision,
intent(out) :: circ(ixi^s,1:
ndim)
4415 integer :: hxc^
l,ixc^
l,ixa^
l
4416 integer :: idim1,idim2,idir,ix^
d
4422 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
4424 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
4425 ixamin^
d=ixcmin^
d+ix^
d;
4426 ixamax^
d=ixcmax^
d+ix^
d;
4427 fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4429 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4435 ixcmin^d=ixomin^d-1;
4443 hxc^l=ixc^l-kr(idim2,^d);
4445 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4446 +lvc(idim1,idim2,idir)&
4451 circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4454 end subroutine update_faces_ambipolar
4458 subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4461 integer,
intent(in) :: ixi^
l, ixo^
l
4462 double precision,
dimension(:^D&,:),
intent(inout) :: ff
4463 double precision,
intent(out) :: src(ixi^s)
4465 double precision :: ffc(ixi^s,1:
ndim)
4466 double precision :: dxinv(
ndim)
4467 integer :: idims, ix^
d, ixa^
l, ixb^
l, ixc^
l
4473 ixcmax^
d=ixomax^
d; ixcmin^
d=ixomin^
d-1;
4475 ixbmin^
d=ixcmin^
d+ix^
d;
4476 ixbmax^
d=ixcmax^
d+ix^
d;
4479 ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4481 ff(ixi^s,1:ndim)=0.d0
4483 ixb^l=ixo^l-kr(idims,^d);
4484 ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4486 if({ ix^d==0 .and. ^d==idims | .or.})
then
4487 ixbmin^d=ixcmin^d-ix^d;
4488 ixbmax^d=ixcmax^d-ix^d;
4489 ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4492 ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4495 if(slab_uniform)
then
4497 ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4498 ixb^l=ixo^l-kr(idims,^d);
4499 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4503 ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4504 ixb^l=ixo^l-kr(idims,^d);
4505 src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4507 src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4509 end subroutine get_flux_on_cell_face
4513 function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x)
result(dtnew)
4516 integer,
intent(in) :: ixi^
l, ixo^
l
4517 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim)
4518 double precision,
intent(in) :: w(ixi^s,1:nw)
4519 double precision :: dtnew
4521 double precision :: coef
4522 double precision :: dxarr(
ndim)
4523 double precision :: tmp(ixi^s)
4528 coef = maxval(abs(tmp(ixo^s)))
4535 dtnew=minval(dxarr(1:
ndim))**2.0d0*coef
4537 dtnew=minval(
block%ds(ixo^s,1:
ndim))**2.0d0*coef
4540 end function get_ambipolar_dt
4548 integer,
intent(in) :: ixi^
l, ixo^
l
4549 double precision,
intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:
ndim)
4550 double precision,
intent(inout) :: res(ixi^s)
4551 double precision :: tmp(ixi^s)
4552 double precision :: rho(ixi^s)
4561 res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4565 subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4572 integer,
intent(in) :: ixi^
l, ixo^
l
4573 double precision,
intent(in) :: qdt,dtfactor
4574 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
4575 double precision,
intent(inout) :: w(ixi^s,1:nw)
4576 logical,
intent(in) :: qsourcesplit
4577 logical,
intent(inout) :: active
4584 if (.not. qsourcesplit)
then
4588 call add_source_internal_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4592 call add_pe0_divv(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4597 call add_hypertc_source(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4603 call add_source_b0split(qdt,dtfactor,ixi^
l,ixo^
l,wctprim,w,x)
4607 if (abs(
mhd_eta)>smalldouble)
then
4609 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
4614 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
4620 call add_source_hydrodynamic_e(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4624 call add_source_semirelativistic(qdt,ixi^
l,ixo^
l,wct,w,x,wctprim)
4631 select case (type_divb)
4636 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4639 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4642 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4643 case (divb_janhunen)
4645 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4646 case (divb_lindejanhunen)
4648 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4649 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4650 case (divb_lindepowel)
4652 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4653 call add_source_powel(qdt,ixi^
l,ixo^
l,wctprim,w,x)
4654 case (divb_lindeglm)
4656 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
4657 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
4658 case (divb_multigrid)
4663 call mpistop(
'Unknown divB fix')
4670 w,x,qsourcesplit,active,
rc_fl)
4680 w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4689 if(.not.qsourcesplit)
then
4691 call mhd_update_temperature(ixi^
l,ixo^
l,wct,w,x)
4695 end subroutine mhd_add_source
4697 subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4701 integer,
intent(in) :: ixi^
l, ixo^
l
4702 double precision,
intent(in) :: qdt,dtfactor
4703 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4704 double precision,
intent(inout) :: w(ixi^s,1:nw)
4705 double precision :: divv(ixi^s)
4721 end subroutine add_pe0_divv
4723 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4725 integer,
intent(in) :: ixi^
l,ixo^
l
4726 double precision,
intent(in) :: qdt
4727 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
4728 double precision,
dimension(ixI^S,1:nw),
intent(in) :: wct,wctprim
4729 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: w
4731 double precision,
dimension(ixI^S) :: r,te,rho_loc
4732 double precision :: sigma_t5,sigma_t7,f_sat,sigmat5_bgradt,tau,bdir(
ndim),bunitvec(
ndim)
4736 call mhd_get_rfactor(wctprim,x,ixi^
l,ixi^
l,r)
4737 te(ixi^s)=wctprim(ixi^s,
p_)/(r(ixi^s)*rho_loc(ixi^s))
4741 do ix2=ixomin2,ixomax2
4742 do ix1=ixomin1,ixomax1
4749 sigma_t7=sigma_t5*te(ix^
d)
4753 sigma_t7=sigma_t5*te(ix^
d)
4756 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))+
block%B0({ix^
d},^
d,0)\
4758 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))\
4760 if(bdir(1)/=0.d0)
then
4761 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2)
4765 if(bdir(2)/=0.d0)
then
4766 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2)
4770 sigmat5_bgradt=sigma_t5*(&
4771 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)&
4772 +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))
4774 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)
4776 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^
d,
q_))/tau
4778 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(sigmat5_bgradt+wct(ix^
d,
q_))/&
4785 do ix3=ixomin3,ixomax3
4786 do ix2=ixomin2,ixomax2
4787 do ix1=ixomin1,ixomax1
4794 sigma_t7=sigma_t5*te(ix^
d)
4798 sigma_t7=sigma_t5*te(ix^
d)
4801 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))+
block%B0({ix^
d},^
d,0)\
4803 ^
d&bdir(^
d)=wct({ix^
d},mag(^
d))\
4805 if(bdir(1)/=0.d0)
then
4806 bunitvec(1)=sign(1.d0,bdir(1))/dsqrt(1.d0+(bdir(2)/bdir(1))**2+(bdir(3)/bdir(1))**2)
4810 if(bdir(2)/=0.d0)
then
4811 bunitvec(2)=sign(1.d0,bdir(2))/dsqrt(1.d0+(bdir(1)/bdir(2))**2+(bdir(3)/bdir(2))**2)
4815 if(bdir(3)/=0.d0)
then
4816 bunitvec(3)=sign(1.d0,bdir(3))/dsqrt(1.d0+(bdir(1)/bdir(3))**2+(bdir(2)/bdir(3))**2)
4820 sigmat5_bgradt=sigma_t5*(&
4821 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)&
4822 +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)&
4823 +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))
4825 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)
4827 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(f_sat*sigmat5_bgradt+wct(ix^
d,
q_))/tau
4829 w(ix^
d,
q_)=w(ix^
d,
q_)-qdt*(sigmat5_bgradt+wct(ix^
d,
q_))/&
4836 end subroutine add_hypertc_source
4839 subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4841 integer,
intent(in) :: ixi^
l, ixo^
l
4842 double precision,
intent(in) :: w(ixi^s,1:nw)
4843 double precision,
intent(inout) :: jxb(ixi^s,3)
4844 double precision :: a(ixi^s,3),
b(ixi^s,3)
4846 double precision :: current(ixi^s,7-2*
ndir:3)
4847 integer :: idir, idirmin
4852 b(ixo^s, idir) = w(ixo^s,mag(idir))+
block%B0(ixo^s,idir,0)
4856 b(ixo^s, idir) = w(ixo^s,mag(idir))
4865 a(ixo^s,idir)=current(ixo^s,idir)
4869 end subroutine get_lorentz_force
4873 subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4875 integer,
intent(in) :: ixi^
l, ixo^
l
4876 double precision,
intent(in) :: w(ixi^s, nw)
4877 double precision,
intent(out) :: gamma_a2(ixo^s)
4878 double precision :: rho(ixi^s)
4884 rho(ixo^s) = w(ixo^s,
rho_)
4887 gamma_a2(ixo^s) = 1.0d0/(1.0d0+
mhd_mag_en_all(w, ixi^
l, ixo^
l)/rho(ixo^s)*inv_squared_c)
4888 end subroutine mhd_gamma2_alfven
4892 function mhd_gamma_alfven(w, ixI^L, ixO^L)
result(gamma_A)
4894 integer,
intent(in) :: ixi^
l, ixo^
l
4895 double precision,
intent(in) :: w(ixi^s, nw)
4896 double precision :: gamma_a(ixo^s)
4898 call mhd_gamma2_alfven(ixi^
l, ixo^
l, w, gamma_a)
4899 gamma_a = sqrt(gamma_a)
4900 end function mhd_gamma_alfven
4904 integer,
intent(in) :: ixi^
l, ixo^
l
4905 double precision,
intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:
ndim)
4906 double precision,
intent(out) :: rho(ixi^s)
4911 rho(ixo^s) = w(ixo^s,
rho_)
4917 subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4920 integer,
intent(in) :: ixi^
l,ixo^
l, ie
4921 double precision,
intent(inout) :: w(ixi^s,1:nw)
4922 double precision,
intent(in) :: x(ixi^s,1:
ndim)
4923 character(len=*),
intent(in) :: subname
4925 double precision :: rho(ixi^s)
4927 logical :: flag(ixi^s,1:nw)
4931 where(w(ixo^s,ie)+
block%equi_vars(ixo^s,
equi_pe0_,0)*inv_gamma_1<small_e)&
4932 flag(ixo^s,ie)=.true.
4934 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4936 if(any(flag(ixo^s,ie)))
then
4940 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4943 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4949 w(ixo^s,
e_)=w(ixo^s,
e_)*gamma_1
4952 w(ixo^s,
mom(idir)) = w(ixo^s,
mom(idir))/rho(ixo^s)
4958 end subroutine mhd_handle_small_ei
4960 subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
4964 integer,
intent(in) :: ixi^
l, ixo^
l
4965 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4966 double precision,
intent(inout) :: w(ixi^s,1:nw)
4968 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
4977 end subroutine mhd_update_temperature
4980 subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4983 integer,
intent(in) :: ixi^
l, ixo^
l
4984 double precision,
intent(in) :: qdt, dtfactor,wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
4985 double precision,
intent(inout) :: w(ixi^s,1:nw)
4987 double precision :: a(ixi^s,3),
b(ixi^s,3), axb(ixi^s,3)
4999 a(ixo^s,idir)=
block%J0(ixo^s,idir)
5004 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5007 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5013 if(total_energy)
then
5016 b(ixo^s,:)=wct(ixo^s,mag(:))
5025 axb(ixo^s,idir)=axb(ixo^s,idir)*
block%dt(ixo^s)*dtfactor
5028 axb(ixo^s,:)=axb(ixo^s,:)*qdt
5032 w(ixo^s,
e_)=w(ixo^s,
e_)-axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5036 call mhd_get_jxbxb(wct,x,ixi^
l,ixo^
l,axb)
5041 w(ixo^s,
e_)=w(ixo^s,
e_)+axb(ixo^s,idir)*
block%J0(ixo^s,idir)
5046 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_B0')
5048 end subroutine add_source_b0split
5051 subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5055 integer,
intent(in) :: ixi^
l, ixo^
l
5056 double precision,
intent(in) :: qdt, 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),
optional :: wctprim(ixi^s,1:nw)
5060 double precision :: e(ixi^s,1:3),curle(ixi^s,1:3),dive(ixi^s)
5061 integer :: idir, idirmin, ix^
d
5065 {
do ix^db=iximin^db,iximax^db\}
5067 e(ix^
d,1)=w(ix^
d,b2_)*wctprim(ix^
d,m3_)-w(ix^
d,b3_)*wctprim(ix^
d,m2_)
5068 e(ix^
d,2)=w(ix^
d,b3_)*wctprim(ix^
d,m1_)-w(ix^
d,b1_)*wctprim(ix^
d,m3_)
5069 e(ix^
d,3)=w(ix^
d,b1_)*wctprim(ix^
d,m2_)-w(ix^
d,b2_)*wctprim(ix^
d,m1_)
5071 call divvector(e,ixi^l,ixo^l,dive)
5073 call curlvector(e,ixi^l,ixo^l,curle,idirmin,1,3)
5076 {
do ix^db=ixomin^db,ixomax^db\}
5077 w(ix^d,m1_)=w(ix^d,m1_)+qdt*(inv_squared_c0-inv_squared_c)*&
5078 (e(ix^d,1)*dive(ix^d)-e(ix^d,2)*curle(ix^d,3)+e(ix^d,3)*curle(ix^d,2))
5079 w(ix^d,m2_)=w(ix^d,m2_)+qdt*(inv_squared_c0-inv_squared_c)*&
5080 (e(ix^d,2)*dive(ix^d)-e(ix^d,3)*curle(ix^d,1)+e(ix^d,1)*curle(ix^d,3))
5081 w(ix^d,m3_)=w(ix^d,m3_)+qdt*(inv_squared_c0-inv_squared_c)*&
5082 (e(ix^d,3)*dive(ix^d)-e(ix^d,1)*curle(ix^d,2)+e(ix^d,2)*curle(ix^d,1) )
5086 end subroutine add_source_semirelativistic
5089 subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5093 integer,
intent(in) :: ixi^
l, ixo^
l
5094 double precision,
intent(in) :: qdt
5095 double precision,
intent(in) :: 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) :: wctprim(ixi^s,1:nw)
5099 double precision :: divv(ixi^s), tmp
5111 {
do ix^db=ixomin^db,ixomax^db\}
5113 w(ix^
d,
e_)=w(ix^
d,
e_)-qdt*wctprim(ix^
d,
p_)*divv(ix^
d)
5114 if(w(ix^
d,
e_)<small_e)
then
5119 call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,
e_)
5122 if(fix_small_values)
then
5123 call mhd_handle_small_ei(w,x,ixi^l,ixo^l,
e_,
'add_source_internal_e')
5125 end subroutine add_source_internal_e
5128 subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
5133 integer,
intent(in) :: ixi^
l, ixo^
l
5134 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5135 double precision,
intent(inout) :: w(ixi^s,1:nw)
5136 double precision,
intent(in),
optional :: wctprim(ixi^s,1:nw)
5138 double precision ::
b(ixi^s,3), j(ixi^s,3), jxb(ixi^s,3)
5139 double precision :: current(ixi^s,7-2*
ndir:3)
5140 double precision :: bu(ixo^s,1:
ndir), tmp(ixo^s), b2(ixo^s)
5141 double precision :: gravity_field(ixi^s,1:
ndir), vaoc
5142 integer :: idir, idirmin, idims, ix^
d
5147 b(ixo^s, idir) = wct(ixo^s,mag(idir))
5155 j(ixo^s,idir)=current(ixo^s,idir)
5173 call gradient(wctprim(ixi^s,
mom(idir)),ixi^
l,ixo^
l,idims,j(ixi^s,idims))
5179 call gradient(wctprim(ixi^s,
p_),ixi^
l,ixo^
l,idir,j(ixi^s,idir))
5186 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)
5190 b(ixo^s,idir)=wct(ixo^s,
rho_)*
b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
5194 b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=
ndim+1)
5195 tmp(ixo^s)=sqrt(b2(ixo^s))
5196 where(tmp(ixo^s)>smalldouble)
5197 tmp(ixo^s)=1.d0/tmp(ixo^s)
5203 bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
5208 {
do ix^db=ixomin^db,ixomax^db\}
5210 vaoc=b2(ix^
d)/w(ix^
d,
rho_)*inv_squared_c
5212 b2(ix^
d)=vaoc/(1.d0+vaoc)
5215 tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*
b(ixo^s,1:ndir),dim=ndim+1)
5218 j(ixo^s,idir)=b2(ixo^s)*(
b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
5222 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+qdt*j(ixo^s,idir)
5225 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*&
5226 (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
5229 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*sum(wctprim(ixo^s,
mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
5232 end subroutine add_source_hydrodynamic_e
5238 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
5243 integer,
intent(in) :: ixi^
l, ixo^
l
5244 double precision,
intent(in) :: qdt
5245 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5246 double precision,
intent(inout) :: w(ixi^s,1:nw)
5247 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
5248 integer :: lxo^
l, kxo^
l
5250 double precision :: tmp(ixi^s),tmp2(ixi^s)
5253 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5254 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
5263 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5264 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
5271 gradeta(ixo^s,1:
ndim)=zero
5277 gradeta(ixo^s,idim)=tmp(ixo^s)
5284 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
5291 tmp2(ixi^s)=bf(ixi^s,idir)
5293 lxo^
l=ixo^
l+2*
kr(idim,^
d);
5294 jxo^
l=ixo^
l+
kr(idim,^
d);
5295 hxo^
l=ixo^
l-
kr(idim,^
d);
5296 kxo^
l=ixo^
l-2*
kr(idim,^
d);
5297 tmp(ixo^s)=tmp(ixo^s)+&
5298 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
5303 tmp2(ixi^s)=bf(ixi^s,idir)
5305 jxo^
l=ixo^
l+
kr(idim,^
d);
5306 hxo^
l=ixo^
l-
kr(idim,^
d);
5307 tmp(ixo^s)=tmp(ixo^s)+&
5308 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
5313 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
5317 do jdir=1,
ndim;
do kdir=idirmin,3
5318 if (
lvc(idir,jdir,kdir)/=0)
then
5319 if (
lvc(idir,jdir,kdir)==1)
then
5320 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5322 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
5329 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
5330 if(total_energy)
then
5331 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
5337 w(ixo^s,
e_)=w(ixo^s,
e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5340 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_res1')
5342 end subroutine add_source_res1
5346 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
5351 integer,
intent(in) :: ixi^
l, ixo^
l
5352 double precision,
intent(in) :: qdt
5353 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5354 double precision,
intent(inout) :: w(ixi^s,1:nw)
5357 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
5358 double precision :: tmpvec(ixi^s,1:3),tmp(ixo^s)
5359 integer :: ixa^
l,idir,idirmin,idirmin1
5363 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5364 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
5374 tmpvec(ixa^s,idir)=current(ixa^s,idir)*
mhd_eta
5379 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
5388 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
5391 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
5396 tmp(ixo^s)=qdt*
mhd_eta*sum(current(ixo^s,:)**2,dim=
ndim+1)
5398 tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=
ndim+1)
5400 if(total_energy)
then
5403 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)-&
5404 qdt*sum(wct(ixo^s,mag(1:
ndir))*curlj(ixo^s,1:
ndir),dim=
ndim+1)
5407 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)
5411 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_res2')
5412 end subroutine add_source_res2
5416 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5420 integer,
intent(in) :: ixi^
l, ixo^
l
5421 double precision,
intent(in) :: qdt
5422 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5423 double precision,
intent(inout) :: w(ixi^s,1:nw)
5425 double precision :: current(ixi^s,7-2*
ndir:3)
5426 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
5427 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
5430 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
5431 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
5434 tmpvec(ixa^s,1:
ndir)=zero
5436 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5440 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5443 tmpvec(ixa^s,1:
ndir)=zero
5444 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
5448 tmpvec2(ixa^s,1:
ndir)=zero
5449 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
5452 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5455 if(total_energy)
then
5458 tmpvec2(ixa^s,1:
ndir)=zero
5459 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
5460 tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5461 +
lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5462 end do;
end do;
end do
5464 call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5465 w(ixo^s,
e_)=w(ixo^s,
e_)+tmp(ixo^s)*qdt
5468 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_hyperres')
5470 end subroutine add_source_hyperres
5472 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5479 integer,
intent(in) :: ixi^
l, ixo^
l
5480 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5481 double precision,
intent(inout) :: w(ixi^s,1:nw)
5483 double precision:: divb(ixi^s), gradpsi(ixi^s), ba(ixo^s,1:
ndir)
5504 ba(ixo^s,1:
ndir)=wct(ixo^s,mag(1:
ndir))
5507 if(total_energy)
then
5516 w(ixo^s,
e_) = w(ixo^s,
e_)-qdt*ba(ixo^s,idir)*gradpsi(ixo^s)
5525 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-qdt*ba(ixo^s,idir)*divb(ixo^s)
5529 if (
fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^
l,ixo^
l,
'add_source_glm')
5531 end subroutine add_source_glm
5534 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5537 integer,
intent(in) :: ixi^
l, ixo^
l
5538 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5539 double precision,
intent(inout) :: w(ixi^s,1:nw)
5541 double precision :: divb(ixi^s), ba(1:
ndir)
5542 integer :: idir, ix^
d
5548 {
do ix^db=ixomin^db,ixomax^db\}
5553 if (total_energy)
then
5559 {
do ix^db=ixomin^db,ixomax^db\}
5561 ^
c&w(ix^d,
b^
c_)=w(ix^d,
b^
c_)-qdt*wct(ix^d,
m^
c_)*divb(ix^d)\
5563 ^
c&w(ix^d,
m^
c_)=w(ix^d,
m^
c_)-qdt*wct(ix^d,
b^
c_)*divb(ix^d)\
5564 if (total_energy)
then
5566 w(ix^d,
e_)=w(ix^d,
e_)-qdt*(^
c&wct(ix^d,
m^
c_)*wct(ix^d,
b^
c_)+)*divb(ix^d)
5571 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_powel')
5573 end subroutine add_source_powel
5575 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5580 integer,
intent(in) :: ixi^
l, ixo^
l
5581 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5582 double precision,
intent(inout) :: w(ixi^s,1:nw)
5584 double precision :: divb(ixi^s)
5585 integer :: idir, ix^
d
5590 {
do ix^db=ixomin^db,ixomax^db\}
5595 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_janhunen')
5597 end subroutine add_source_janhunen
5599 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5604 integer,
intent(in) :: ixi^
l, ixo^
l
5605 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
5606 double precision,
intent(inout) :: w(ixi^s,1:nw)
5608 double precision :: divb(ixi^s),graddivb(ixi^s)
5609 integer :: idim, idir, ixp^
l, i^
d, iside
5610 logical,
dimension(-1:1^D&) :: leveljump
5618 if(i^
d==0|.and.) cycle
5619 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
5620 leveljump(i^
d)=.true.
5622 leveljump(i^
d)=.false.
5631 i^dd=kr(^dd,^d)*(2*iside-3);
5632 if (leveljump(i^dd))
then
5634 ixpmin^d=ixomin^d-i^d
5636 ixpmax^d=ixomax^d-i^d
5647 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5649 {
do i^db=ixpmin^db,ixpmax^db\}
5651 graddivb(i^d)=graddivb(i^d)*divbdiff/(^d&1.0d0/block%ds({i^d},^d)**2+)
5653 w(i^d,mag(idim))=w(i^d,mag(idim))+graddivb(i^d)
5655 if (typedivbdiff==
'all' .and. total_energy)
then
5657 w(i^d,
e_)=w(i^d,
e_)+wct(i^d,mag(idim))*graddivb(i^d)
5662 if (fix_small_values)
call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,
'add_source_linde')
5664 end subroutine add_source_linde
5671 integer,
intent(in) :: ixi^
l, ixo^
l
5672 double precision,
intent(in) :: w(ixi^s,1:nw)
5673 double precision :: divb(ixi^s), dsurface(ixi^s)
5675 double precision :: invb(ixo^s)
5676 integer :: ixa^
l,idims
5678 call get_divb(w,ixi^
l,ixo^
l,divb)
5680 where(invb(ixo^s)/=0.d0)
5681 invb(ixo^s)=1.d0/invb(ixo^s)
5684 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
5686 ixamin^
d=ixomin^
d-1;
5687 ixamax^
d=ixomax^
d-1;
5688 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
5690 ixa^
l=ixo^
l-
kr(idims,^
d);
5691 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
5693 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5694 block%dvolume(ixo^s)/dsurface(ixo^s)
5705 integer,
intent(in) :: ixo^
l, ixi^
l
5706 double precision,
intent(in) :: w(ixi^s,1:nw)
5707 integer,
intent(out) :: idirmin
5710 double precision :: current(ixi^s,7-2*
ndir:3)
5711 integer :: idir, idirmin0
5717 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5718 block%J0(ixo^s,idirmin0:3)
5722 subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5730 integer,
intent(in) :: ixi^
l, ixo^
l
5731 double precision,
intent(inout) :: dtnew
5732 double precision,
intent(in) ::
dx^
d
5733 double precision,
intent(in) :: w(ixi^s,1:nw)
5734 double precision,
intent(in) :: x(ixi^s,1:
ndim)
5736 double precision :: dxarr(
ndim)
5737 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
5738 integer :: idirmin,idim
5752 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5755 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
5781 dtnew=min(
dtdiffpar*get_ambipolar_dt(w,ixi^
l,ixo^
l,
dx^
d,x),dtnew)
5788 end subroutine mhd_get_dt
5791 subroutine mhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5796 integer,
intent(in) :: ixi^
l, ixo^
l
5797 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5798 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5800 double precision :: tmp,tmp1,invr,cot
5802 integer :: mr_,mphi_
5803 integer :: br_,bphi_
5806 br_=mag(1); bphi_=mag(1)-1+
phi_
5811 {
do ix^db=ixomin^db,ixomax^db\}
5814 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5819 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
5824 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
5825 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
5826 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5827 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5828 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
5830 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5831 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5832 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5835 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
5840 {
do ix^db=ixomin^db,ixomax^db\}
5842 if(local_timestep)
then
5843 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
5848 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
5854 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
5857 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
5858 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
5862 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
5868 cot=1.d0/tan(x(ix^d,2))
5872 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5873 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
5875 if(.not.stagger_grid)
then
5876 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5878 tmp=tmp+wprim(ix^d,
psi_)*cot
5880 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5885 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
5886 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
5887 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
5889 if(.not.stagger_grid)
then
5890 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
5892 tmp=tmp+wprim(ix^d,
psi_)*cot
5894 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
5897 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
5898 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
5899 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
5900 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
5901 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
5903 if(.not.stagger_grid)
then
5904 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
5905 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
5906 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
5907 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
5908 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
5915 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
5918 end subroutine mhd_add_source_geom
5921 subroutine mhd_add_source_geom_semirelati(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
5926 integer,
intent(in) :: ixi^
l, ixo^
l
5927 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
5928 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
5930 double precision :: tmp,tmp1,tmp2,invr,cot,e(ixo^s,1:
ndir)
5932 integer :: mr_,mphi_
5933 integer :: br_,bphi_
5936 br_=mag(1); bphi_=mag(1)-1+
phi_
5941 {
do ix^db=ixomin^db,ixomax^db\}
5944 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
5955 e(ix^
d,1)=wprim(ix^
d,b2_)*wprim(ix^
d,m3_)-wprim(ix^
d,b3_)*wprim(ix^
d,m2_)
5956 e(ix^
d,2)=wprim(ix^
d,b3_)*wprim(ix^
d,m1_)-wprim(ix^
d,b1_)*wprim(ix^
d,m3_)
5957 e(ix^
d,3)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5962 e(ix^
d,2)=wprim(ix^
d,b1_)*wprim(ix^
d,m2_)-wprim(ix^
d,b2_)*wprim(ix^
d,m1_)
5968 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+&
5969 half*((^
c&wprim(ix^
d,
b^
c_)**2+)+(^
c&e(ix^
d,^
c)**2+)*inv_squared_c) -&
5970 wprim(ix^
d,bphi_)**2+wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)**2)
5971 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
5972 -wprim(ix^
d,
rho_)*wprim(ix^
d,mphi_)*wprim(ix^
d,mr_) &
5973 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_)+e(ix^
d,
phi_)*e(ix^
d,1)*inv_squared_c)
5975 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
5976 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
5977 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
5980 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp+half*((^
c&wprim(ix^
d,
b^
c_)**2+)+&
5981 (^
c&e(ix^
d,^
c)**2+)*inv_squared_c))
5986 {
do ix^db=ixomin^db,ixomax^db\}
5988 if(local_timestep)
then
5989 invr=block%dt(ix^d)*dtfactor/x(ix^d,1)
5995 e(ix^d,1)=wprim(ix^d,b2_)*wprim(ix^d,m3_)-wprim(ix^d,b3_)*wprim(ix^d,m2_)
5996 e(ix^d,2)=wprim(ix^d,b3_)*wprim(ix^d,m1_)-wprim(ix^d,b1_)*wprim(ix^d,m3_)
5997 e(ix^d,3)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6001 e(ix^d,1)=wprim(ix^d,b1_)*wprim(ix^d,m2_)-wprim(ix^d,b2_)*wprim(ix^d,m1_)
6008 tmp1=wprim(ix^d,
p_)+half*((^
c&wprim(ix^d,
b^
c_)**2+)+(^
c&e(ix^d,^
c)**2+)*inv_squared_c)
6014 w(ix^d,m1_)=w(ix^d,m1_)+two*tmp1*invr
6017 w(ix^d,m1_)=w(ix^d,m1_)+invr*&
6018 (two*tmp1+(^ce&wprim(ix^d,
rho_)*wprim(ix^d,
m^ce_)**2-&
6019 wprim(ix^d,
b^ce_)**2-e(ix^d,^ce)**2*inv_squared_c+))
6023 w(ix^d,b1_)=w(ix^d,b1_)+invr*2.0d0*wprim(ix^d,
psi_)
6029 cot=1.d0/tan(x(ix^d,2))
6033 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_)&
6034 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c)
6036 if(.not.stagger_grid)
then
6037 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6039 tmp=tmp+wprim(ix^d,
psi_)*cot
6041 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6047 w(ix^d,m2_)=w(ix^d,m2_)+invr*(tmp1*cot-wprim(ix^d,
rho_)*wprim(ix^d,m1_)*wprim(ix^d,m2_) &
6048 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+e(ix^d,1)*e(ix^d,2)*inv_squared_c&
6049 +(wprim(ix^d,
rho_)*wprim(ix^d,m3_)**2&
6050 -wprim(ix^d,b3_)**2-e(ix^d,3)**2*inv_squared_c)*cot)
6052 if(.not.stagger_grid)
then
6053 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6055 tmp=tmp+wprim(ix^d,
psi_)*cot
6057 w(ix^d,b2_)=w(ix^d,b2_)+tmp*invr
6060 w(ix^d,m3_)=w(ix^d,m3_)+invr*&
6061 (-wprim(ix^d,m3_)*wprim(ix^d,m1_)*wprim(ix^d,
rho_) &
6062 +wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6063 +e(ix^d,3)*e(ix^d,1)*inv_squared_c&
6064 +(-wprim(ix^d,m2_)*wprim(ix^d,m3_)*wprim(ix^d,
rho_) &
6065 +wprim(ix^d,b2_)*wprim(ix^d,b3_)&
6066 +e(ix^d,2)*e(ix^d,3)*inv_squared_c)*cot)
6068 if(.not.stagger_grid)
then
6069 w(ix^d,b3_)=w(ix^d,b3_)+invr*&
6070 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6071 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6072 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6073 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6080 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6083 end subroutine mhd_add_source_geom_semirelati
6086 subroutine mhd_add_source_geom_split(qdt,dtfactor,ixI^L,ixO^L,wCT,wprim,w,x)
6091 integer,
intent(in) :: ixi^
l, ixo^
l
6092 double precision,
intent(in) :: qdt, dtfactor,x(ixi^s,1:
ndim)
6093 double precision,
intent(inout) :: wct(ixi^s,1:nw),wprim(ixi^s,1:nw),w(ixi^s,1:nw)
6095 double precision :: tmp,tmp1,tmp2,invr,cot
6097 integer :: mr_,mphi_
6098 integer :: br_,bphi_
6101 br_=mag(1); bphi_=mag(1)-1+
phi_
6106 {
do ix^db=ixomin^db,ixomax^db\}
6109 invr=
block%dt(ix^
d) * dtfactor/x(ix^
d,1)
6114 tmp=wprim(ix^
d,
p_)+half*(^
c&wprim(ix^
d,
b^
c_)**2+)
6119 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*(tmp-&
6120 wprim(ix^
d,bphi_)**2+wprim(ix^
d,mphi_)*wct(ix^
d,mphi_))
6121 w(ix^
d,mphi_)=w(ix^
d,mphi_)+invr*(&
6122 -wct(ix^
d,mphi_)*wprim(ix^
d,mr_) &
6123 +wprim(ix^
d,bphi_)*wprim(ix^
d,br_))
6125 w(ix^
d,bphi_)=w(ix^
d,bphi_)+invr*&
6126 (wprim(ix^
d,bphi_)*wprim(ix^
d,mr_) &
6127 -wprim(ix^
d,br_)*wprim(ix^
d,mphi_))
6130 w(ix^
d,mr_)=w(ix^
d,mr_)+invr*tmp
6135 {
do ix^db=ixomin^db,ixomax^db\}
6137 if(local_timestep)
then
6138 invr=block%dt(ix^d) * dtfactor/x(ix^d,1)
6142 tmp1=wprim(ix^d,
p_)+half*(^
c&wprim(ix^d,
b^
c_)**2+)
6143 if(b0field) tmp2=(^
c&block%B0(ix^d,^
c,0)*wprim(ix^d,
b^
c_)+)
6146 w(ix^d,
mom(1))=w(ix^d,
mom(1))+two*tmp1*invr
6150 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6151 (two*(tmp1+tmp2)+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+)- &
6152 (^ce&two*block%B0(ix^d,^ce,0)*wprim(ix^d,
b^ce_)+))
6154 w(ix^d,
mom(1))=w(ix^d,
mom(1))+invr*&
6155 (two*tmp1+(^ce&wprim(ix^d,
m^ce_)*wct(ix^d,
m^ce_)-wprim(ix^d,
b^ce_)**2+))
6160 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wprim(ix^d,
psi_)
6166 cot=1.d0/tan(x(ix^d,2))
6171 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6172 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6173 +wprim(ix^d,b1_)*block%B0(ix^d,2,0))
6175 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6176 +wprim(ix^d,b1_)*wprim(ix^d,b2_))
6179 if(.not.stagger_grid)
then
6181 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6182 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6184 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6187 tmp=tmp+wprim(ix^d,
psi_)*cot
6189 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6195 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*((tmp1+tmp2)*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6196 +wprim(ix^d,b1_)*wprim(ix^d,b2_)+block%B0(ix^d,1,0)*wprim(ix^d,b2_)&
6197 +wprim(ix^d,b1_)*block%B0(ix^d,2,0)&
6198 +(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)
6200 w(ix^d,
mom(2))=w(ix^d,
mom(2))+invr*(tmp1*cot-wprim(ix^d,m1_)*wct(ix^d,m2_)&
6201 +wprim(ix^d,b1_)*wprim(ix^d,b2_)&
6202 +(wprim(ix^d,m3_)*wct(ix^d,m3_)-wprim(ix^d,b3_)**2)*cot)
6205 if(.not.stagger_grid)
then
6207 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)&
6208 +wprim(ix^d,m1_)*block%B0(ix^d,2,0)-wprim(ix^d,m2_)*block%B0(ix^d,1,0)
6210 tmp=wprim(ix^d,m1_)*wprim(ix^d,b2_)-wprim(ix^d,m2_)*wprim(ix^d,b1_)
6213 tmp=tmp+wprim(ix^d,
psi_)*cot
6215 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
6219 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6220 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6221 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6222 +block%B0(ix^d,1,0)*wprim(ix^d,b3_) &
6223 +wprim(ix^d,b1_)*block%B0(ix^d,3,0) &
6224 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6225 -wprim(ix^d,b2_)*wprim(ix^d,b3_) &
6226 +block%B0(ix^d,2,0)*wprim(ix^d,b3_) &
6227 +wprim(ix^d,b2_)*block%B0(ix^d,3,0))*cot)
6229 w(ix^d,
mom(3))=w(ix^d,
mom(3))-invr*&
6230 (wprim(ix^d,m3_)*wct(ix^d,m1_) &
6231 -wprim(ix^d,b3_)*wprim(ix^d,b1_) &
6232 +(wprim(ix^d,m2_)*wct(ix^d,m3_) &
6233 -wprim(ix^d,b2_)*wprim(ix^d,b3_))*cot)
6236 if(.not.stagger_grid)
then
6238 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6239 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6240 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6241 +wprim(ix^d,m1_)*block%B0(ix^d,3,0) &
6242 -wprim(ix^d,m3_)*block%B0(ix^d,1,0) &
6243 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6244 -wprim(ix^d,m2_)*wprim(ix^d,b3_) &
6245 +wprim(ix^d,m3_)*block%B0(ix^d,2,0) &
6246 -wprim(ix^d,m2_)*block%B0(ix^d,3,0))*cot)
6248 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
6249 (wprim(ix^d,m1_)*wprim(ix^d,b3_) &
6250 -wprim(ix^d,m3_)*wprim(ix^d,b1_) &
6251 -(wprim(ix^d,m3_)*wprim(ix^d,b2_) &
6252 -wprim(ix^d,m2_)*wprim(ix^d,b3_))*cot)
6260 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
6263 end subroutine mhd_add_source_geom_split
6268 integer,
intent(in) :: ixi^
l, ixo^
l
6269 double precision,
intent(in) :: w(ixi^s, nw)
6270 double precision :: mge(ixo^s)
6273 mge = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=
ndim+1)
6275 mge = sum(w(ixo^s, mag(:))**2, dim=
ndim+1)
6279 subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
6282 integer,
intent(in) :: ixi^
l, ixo^
l
6283 double precision,
intent(in) :: w(ixi^s,nw)
6284 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6285 double precision,
intent(inout) :: vhall(ixi^s,1:
ndir)
6287 double precision :: current(ixi^s,7-2*
ndir:3)
6288 double precision :: rho(ixi^s)
6289 integer :: idir, idirmin, ix^
d
6294 do idir = idirmin,
ndir
6295 {
do ix^db=ixomin^db,ixomax^db\}
6296 vhall(ix^
d,idir)=-
mhd_etah*current(ix^
d,idir)/rho(ix^
d)
6300 end subroutine mhd_getv_hall
6302 subroutine mhd_get_jambi(w,x,ixI^L,ixO^L,res)
6305 integer,
intent(in) :: ixi^
l, ixo^
l
6306 double precision,
intent(in) :: w(ixi^s,nw)
6307 double precision,
intent(in) :: x(ixi^s,1:
ndim)
6308 double precision,
allocatable,
intent(inout) :: res(:^
d&,:)
6311 double precision :: current(ixi^s,7-2*
ndir:3)
6312 integer :: idir, idirmin
6319 res(ixo^s,idirmin:3)=-current(ixo^s,idirmin:3)
6320 do idir = idirmin, 3
6324 end subroutine mhd_get_jambi
6326 subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
6329 integer,
intent(in) :: ixi^
l, ixo^
l, idir
6330 double precision,
intent(in) :: qt
6331 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
6332 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
6335 double precision :: db(ixo^s), dpsi(ixo^s)
6339 {
do ix^db=ixomin^db,ixomax^db\}
6340 wlc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6341 wrc(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6342 wlp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6343 wrp(ix^
d,mag(idir))=s%ws(ix^
d,idir)
6352 {
do ix^db=ixomin^db,ixomax^db\}
6353 db(ix^d)=wrp(ix^d,mag(idir))-wlp(ix^d,mag(idir))
6354 dpsi(ix^d)=wrp(ix^d,
psi_)-wlp(ix^d,
psi_)
6355 wlp(ix^d,mag(idir))=half*(wrp(ix^d,mag(idir))+wlp(ix^d,mag(idir))-dpsi(ix^d)/cmax_global)
6356 wlp(ix^d,
psi_)=half*(wrp(ix^d,
psi_)+wlp(ix^d,
psi_)-db(ix^d)*cmax_global)
6357 wrp(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6359 if(total_energy)
then
6360 wrc(ix^d,
e_)=wrc(ix^d,
e_)-half*wrc(ix^d,mag(idir))**2
6361 wlc(ix^d,
e_)=wlc(ix^d,
e_)-half*wlc(ix^d,mag(idir))**2
6363 wrc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6365 wlc(ix^d,mag(idir))=wlp(ix^d,mag(idir))
6368 if(total_energy)
then
6369 wrc(ix^d,
e_)=wrc(ix^d,
e_)+half*wrc(ix^d,mag(idir))**2
6370 wlc(ix^d,
e_)=wlc(ix^d,
e_)+half*wlc(ix^d,mag(idir))**2
6375 if(
associated(usr_set_wlr))
call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
6377 end subroutine mhd_modify_wlr
6379 subroutine mhd_boundary_adjust(igrid,psb)
6381 integer,
intent(in) :: igrid
6384 integer :: ib, idims, iside, ixo^
l, i^
d
6393 i^
d=
kr(^
d,idims)*(2*iside-3);
6394 if (neighbor_type(i^
d,igrid)/=1) cycle
6395 ib=(idims-1)*2+iside
6413 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
6418 end subroutine mhd_boundary_adjust
6420 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
6423 integer,
intent(in) :: ixg^
l,ixo^
l,ib
6424 double precision,
intent(inout) :: w(ixg^s,1:nw)
6425 double precision,
intent(in) :: x(ixg^s,1:
ndim)
6427 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
6428 integer :: ix^
d,ixf^
l
6441 do ix1=ixfmax1,ixfmin1,-1
6442 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
6443 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6444 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6447 do ix1=ixfmax1,ixfmin1,-1
6448 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
6449 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
6450 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6451 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6452 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6453 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6454 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6468 do ix1=ixfmax1,ixfmin1,-1
6469 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6470 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6471 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6472 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6473 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6474 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6477 do ix1=ixfmax1,ixfmin1,-1
6478 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6479 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6480 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6481 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6482 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6483 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6484 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6485 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6486 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6487 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6488 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6489 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6490 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6491 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6492 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6493 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6494 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6495 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6509 do ix1=ixfmin1,ixfmax1
6510 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
6511 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
6512 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
6515 do ix1=ixfmin1,ixfmax1
6516 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
6517 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
6518 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
6519 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
6520 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
6521 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
6522 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
6536 do ix1=ixfmin1,ixfmax1
6537 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6538 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
6539 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
6540 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
6541 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
6542 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
6545 do ix1=ixfmin1,ixfmax1
6546 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
6547 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
6548 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
6549 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
6550 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
6551 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
6552 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
6553 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
6554 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
6555 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
6556 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
6557 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
6558 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
6559 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
6560 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6561 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
6562 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
6563 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
6577 do ix2=ixfmax2,ixfmin2,-1
6578 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
6579 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6580 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6583 do ix2=ixfmax2,ixfmin2,-1
6584 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
6585 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
6586 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6587 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6588 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6589 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6590 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6604 do ix2=ixfmax2,ixfmin2,-1
6605 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6606 ix2+1,ixfmin3:ixfmax3,mag(2)) &
6607 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6608 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6609 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6610 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6613 do ix2=ixfmax2,ixfmin2,-1
6614 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
6615 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
6616 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6617 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
6618 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6619 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6620 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6621 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6622 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6623 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6624 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6625 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6626 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6627 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6628 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6629 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6630 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
6631 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6645 do ix2=ixfmin2,ixfmax2
6646 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
6647 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
6648 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
6651 do ix2=ixfmin2,ixfmax2
6652 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
6653 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
6654 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
6655 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
6656 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
6657 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
6658 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
6672 do ix2=ixfmin2,ixfmax2
6673 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
6674 ix2-1,ixfmin3:ixfmax3,mag(2)) &
6675 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
6676 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
6677 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
6678 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
6681 do ix2=ixfmin2,ixfmax2
6682 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
6683 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
6684 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
6685 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
6686 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
6687 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6688 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
6689 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
6690 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
6691 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
6692 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
6693 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
6694 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
6695 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
6696 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
6697 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
6698 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
6699 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
6716 do ix3=ixfmax3,ixfmin3,-1
6717 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
6718 ixfmin2:ixfmax2,ix3+1,mag(3)) &
6719 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6720 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6721 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6722 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6725 do ix3=ixfmax3,ixfmin3,-1
6726 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
6727 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
6728 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6729 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
6730 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6731 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6732 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6733 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6734 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6735 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6736 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6737 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6738 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6739 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6740 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6741 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6742 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
6743 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6758 do ix3=ixfmin3,ixfmax3
6759 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
6760 ixfmin2:ixfmax2,ix3-1,mag(3)) &
6761 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
6762 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
6763 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
6764 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
6767 do ix3=ixfmin3,ixfmax3
6768 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
6769 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
6770 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
6771 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
6772 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
6773 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6774 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
6775 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
6776 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
6777 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
6778 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
6779 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
6780 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
6781 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
6782 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
6783 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
6784 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
6785 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
6791 call mpistop(
"Special boundary is not defined for this region")
6794 end subroutine fixdivb_boundary
6803 double precision,
intent(in) :: qdt
6804 double precision,
intent(in) :: qt
6805 logical,
intent(inout) :: active
6808 integer,
parameter :: max_its = 50
6809 double precision :: residual_it(max_its), max_divb
6810 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
6811 double precision :: res
6812 double precision,
parameter :: max_residual = 1
d-3
6813 double precision,
parameter :: residual_reduction = 1
d-10
6814 integer :: iigrid, igrid
6815 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
6818 mg%operator_type = mg_laplacian
6826 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6827 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6830 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
6831 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6833 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6834 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6837 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6838 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6842 write(*,*)
"mhd_clean_divb_multigrid warning: unknown boundary type"
6843 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
6844 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
6852 do iigrid = 1, igridstail
6853 igrid = igrids(iigrid);
6856 lvl =
mg%boxes(id)%lvl
6857 nc =
mg%box_size_lvl(lvl)
6863 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
6865 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
6866 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
6871 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
6874 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
6875 if (
mype == 0) print *,
"iteration vs residual"
6878 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
6879 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
6880 if (residual_it(n) < residual_reduction * max_divb)
exit
6882 if (
mype == 0 .and. n > max_its)
then
6883 print *,
"divb_multigrid warning: not fully converged"
6884 print *,
"current amplitude of divb: ", residual_it(max_its)
6885 print *,
"multigrid smallest grid: ", &
6886 mg%domain_size_lvl(:,
mg%lowest_lvl)
6887 print *,
"note: smallest grid ideally has <= 8 cells"
6888 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
6889 print *,
"note: dx/dy/dz should be similar"
6893 call mg_fas_vcycle(
mg, max_res=res)
6894 if (res < max_residual)
exit
6896 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
6901 do iigrid = 1, igridstail
6902 igrid = igrids(iigrid);
6911 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
6915 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
6917 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
6919 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
6932 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
6933 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
6936 if(total_energy)
then
6938 tmp(
ixm^t) = 0.5_dp * (sum(ps(igrid)%w(
ixm^t, &
6941 ps(igrid)%w(
ixm^t,
e_) = ps(igrid)%w(
ixm^t,
e_) + tmp(
ixm^t)
6951 subroutine mhd_update_faces_average(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
6955 integer,
intent(in) :: ixi^
l, ixo^
l
6956 double precision,
intent(in) :: qt,qdt
6958 double precision,
intent(in) :: wp(ixi^s,1:nw)
6959 type(state) :: sct, s
6960 type(ct_velocity) :: vcts
6961 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
6962 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
6964 double precision :: circ(ixi^s,1:
ndim)
6966 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
6967 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
6968 integer :: idim1,idim2,idir,iwdim1,iwdim2
6970 associate(bfaces=>s%ws,x=>s%x)
6977 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
6984 i1kr^
d=
kr(idim1,^
d);
6987 i2kr^
d=
kr(idim2,^
d);
6990 if (
lvc(idim1,idim2,idir)==1)
then
6992 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
6994 {
do ix^db=ixcmin^db,ixcmax^db\}
6995 fe(ix^
d,idir)=quarter*&
6996 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
6997 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
6999 if(
mhd_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
7004 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
7012 if(
associated(usr_set_electric_field)) &
7013 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7015 circ(ixi^s,1:ndim)=zero
7020 ixcmin^d=ixomin^d-kr(idim1,^d);
7022 ixa^l=ixc^l-kr(idim2,^d);
7025 if(lvc(idim1,idim2,idir)==1)
then
7027 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7030 else if(lvc(idim1,idim2,idir)==-1)
then
7032 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7038 {
do ix^db=ixcmin^db,ixcmax^db\}
7040 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7042 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7049 end subroutine mhd_update_faces_average
7052 subroutine mhd_update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7057 integer,
intent(in) :: ixi^
l, ixo^
l
7058 double precision,
intent(in) :: qt, qdt
7060 double precision,
intent(in) :: wp(ixi^s,1:nw)
7061 type(state) :: sct, s
7062 type(ct_velocity) :: vcts
7063 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7064 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7066 double precision :: circ(ixi^s,1:
ndim)
7068 double precision :: ecc(ixi^s,
sdim:3)
7069 double precision :: ein(ixi^s,
sdim:3)
7071 double precision :: el(ixi^s),er(ixi^s)
7073 double precision :: elc,erc
7075 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7077 double precision :: jce(ixi^s,
sdim:3)
7079 double precision :: xs(ixgs^t,1:
ndim)
7080 double precision :: gradi(ixgs^t)
7081 integer :: ixc^
l,ixa^
l
7082 integer :: idim1,idim2,idir,iwdim1,iwdim2,ix^
d,i1kr^
d,i2kr^
d
7084 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm,wcts=>sct%ws)
7087 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7093 {
do ix^db=iximin^db,iximax^db\}
7096 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_)
7097 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_)
7098 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_)
7101 ecc(ix^
d,3)=wp(ix^
d,b1_)*wp(ix^
d,m2_)-wp(ix^
d,b2_)*wp(ix^
d,m1_)
7108 {
do ix^db=iximin^db,iximax^db\}
7111 ecc(ix^d,1)=wp(ix^d,b2_)*wp(ix^d,m3_)-wp(ix^d,b3_)*wp(ix^d,m2_)
7112 ecc(ix^d,2)=wp(ix^d,b3_)*wp(ix^d,m1_)-wp(ix^d,b1_)*wp(ix^d,m3_)
7113 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7116 ecc(ix^d,3)=wp(ix^d,b1_)*wp(ix^d,m2_)-wp(ix^d,b2_)*wp(ix^d,m1_)
7130 i1kr^d=kr(idim1,^d);
7133 i2kr^d=kr(idim2,^d);
7136 if (lvc(idim1,idim2,idir)==1)
then
7138 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7141 {
do ix^db=ixcmin^db,ixcmax^db\}
7142 fe(ix^d,idir)=quarter*&
7143 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
7144 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
7149 ixamax^d=ixcmax^d+i1kr^d;
7150 {
do ix^db=ixamin^db,ixamax^db\}
7151 el(ix^d)=fc(ix^d,iwdim1,idim2)-ecc(ix^d,idir)
7152 er(ix^d)=fc(ix^d,iwdim1,idim2)-ecc({ix^d+i2kr^d},idir)
7155 do ix^db=ixcmin^db,ixcmax^db\}
7156 if(vnorm(ix^d,idim1)>0.d0)
then
7158 else if(vnorm(ix^d,idim1)<0.d0)
then
7159 elc=el({ix^d+i1kr^d})
7161 elc=0.5d0*(el(ix^d)+el({ix^d+i1kr^d}))
7163 if(vnorm({ix^d+i2kr^d},idim1)>0.d0)
then
7165 else if(vnorm({ix^d+i2kr^d},idim1)<0.d0)
then
7166 erc=er({ix^d+i1kr^d})
7168 erc=0.5d0*(er(ix^d)+er({ix^d+i1kr^d}))
7170 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7175 ixamax^d=ixcmax^d+i2kr^d;
7176 {
do ix^db=ixamin^db,ixamax^db\}
7177 el(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc(ix^d,idir)
7178 er(ix^d)=-fc(ix^d,iwdim2,idim1)-ecc({ix^d+i1kr^d},idir)
7181 do ix^db=ixcmin^db,ixcmax^db\}
7182 if(vnorm(ix^d,idim2)>0.d0)
then
7184 else if(vnorm(ix^d,idim2)<0.d0)
then
7185 elc=el({ix^d+i2kr^d})
7187 elc=0.5d0*(el(ix^d)+el({ix^d+i2kr^d}))
7189 if(vnorm({ix^d+i1kr^d},idim2)>0.d0)
then
7191 else if(vnorm({ix^d+i1kr^d},idim2)<0.d0)
then
7192 erc=er({ix^d+i2kr^d})
7194 erc=0.5d0*(er(ix^d)+er({ix^d+i2kr^d}))
7196 fe(ix^d,idir)=fe(ix^d,idir)+0.25d0*(elc+erc)
7200 if(
mhd_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
7205 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
7219 if (lvc(idim1,idim2,idir)==0) cycle
7221 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7222 ixamax^d=ixcmax^d-kr(idir,^d)+1;
7225 xs(ixa^s,:)=x(ixa^s,:)
7226 xs(ixa^s,idim2)=x(ixa^s,idim2)+half*s%dx(ixa^s,idim2)
7227 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
7228 if (lvc(idim1,idim2,idir)==1)
then
7229 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7231 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7238 ixcmin^d=ixomin^d+kr(idir,^d)-1;
7240 ein(ixc^s,idir)=ein(ixc^s,idir)*jce(ixc^s,idir)
7244 {
do ix^db=ixomin^db,ixomax^db\}
7245 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1,ix2-1,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7246 +ein(ix1,ix2-1,ix3-1,idir))
7247 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7248 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7250 else if(idir==2)
then
7251 {
do ix^db=ixomin^db,ixomax^db\}
7252 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2,ix3-1,idir)&
7253 +ein(ix1-1,ix2,ix3-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)
7258 {
do ix^db=ixomin^db,ixomax^db\}
7259 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,ix3,idir)+ein(ix1,ix2-1,ix3,idir)&
7260 +ein(ix1-1,ix2-1,ix3,idir))
7261 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7262 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7268 {
do ix^db=ixomin^db,ixomax^db\}
7269 jce(ix^d,idir)=0.25d0*(ein(ix^d,idir)+ein(ix1-1,ix2,idir)+ein(ix1,ix2-1,idir)&
7270 +ein(ix1-1,ix2-1,idir))
7271 if(jce(ix^d,idir)<0.d0) jce(ix^d,idir)=0.d0
7272 w(ix^d,
e_)=w(ix^d,
e_)+qdt*jce(ix^d,idir)
7277 block%w(ixo^s,nw)=block%w(ixo^s,nw)+jce(ixo^s,idir)
7283 if(
associated(usr_set_electric_field)) &
7284 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
7286 circ(ixi^s,1:ndim)=zero
7291 ixcmin^d=ixomin^d-kr(idim1,^d);
7293 ixa^l=ixc^l-kr(idim2,^d);
7296 if(lvc(idim1,idim2,idir)==1)
then
7298 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7301 else if(lvc(idim1,idim2,idir)==-1)
then
7303 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7309 {
do ix^db=ixcmin^db,ixcmax^db\}
7311 if(s%surfaceC(ix^d,idim1) > smalldouble)
then
7313 bfaces(ix^d,idim1)=bfaces(ix^d,idim1)-circ(ix^d,idim1)/s%surfaceC(ix^d,idim1)
7320 end subroutine mhd_update_faces_contact
7323 subroutine mhd_update_faces_hll(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
7328 integer,
intent(in) :: ixi^
l, ixo^
l
7329 double precision,
intent(in) :: qt, qdt
7331 double precision,
intent(in) :: wp(ixi^s,1:nw)
7332 type(state) :: sct, s
7333 type(ct_velocity) :: vcts
7334 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
7335 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
7337 double precision :: vtill(ixi^s,2)
7338 double precision :: vtilr(ixi^s,2)
7339 double precision :: bfacetot(ixi^s,
ndim)
7340 double precision :: btill(ixi^s,
ndim)
7341 double precision :: btilr(ixi^s,
ndim)
7342 double precision :: cp(ixi^s,2)
7343 double precision :: cm(ixi^s,2)
7344 double precision :: circ(ixi^s,1:
ndim)
7346 double precision,
dimension(ixI^S,sdim:3) :: e_resi, e_ambi
7347 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
7348 integer :: idim1,idim2,idir,ix^
d
7350 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
7351 cbarmax=>vcts%cbarmax)
7364 if(
mhd_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,wp,sct,s,e_resi)
7380 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
7384 idim2=mod(idir+1,3)+1
7386 jxc^
l=ixc^
l+
kr(idim1,^
d);
7387 ixcp^
l=ixc^
l+
kr(idim2,^
d);
7391 vtill(ixi^s,2),vtilr(ixi^s,2))
7394 vtill(ixi^s,1),vtilr(ixi^s,1))
7400 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+
block%B0(ixi^s,idim1,idim1)
7401 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+
block%B0(ixi^s,idim2,idim2)
7403 bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
7404 bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
7407 btill(ixi^s,idim1),btilr(ixi^s,idim1))
7410 btill(ixi^s,idim2),btilr(ixi^s,idim2))
7414 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
7415 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
7417 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
7418 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
7422 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
7423 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
7424 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
7425 /(cp(ixc^s,1)+cm(ixc^s,1)) &
7426 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
7427 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
7428 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
7429 /(cp(ixc^s,2)+cm(ixc^s,2))
7432 if(
mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+e_resi(ixc^s,idir)
7436 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
7450 circ(ixi^s,1:
ndim)=zero
7455 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
7459 if(
lvc(idim1,idim2,idir)/=0)
then
7460 hxc^
l=ixc^
l-
kr(idim2,^
d);
7462 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
7463 +
lvc(idim1,idim2,idir)&
7469 {
do ix^db=ixcmin^db,ixcmax^db\}
7471 if(s%surfaceC(ix^
d,idim1) > smalldouble)
then
7473 bfaces(ix^
d,idim1)=bfaces(ix^
d,idim1)-circ(ix^
d,idim1)/s%surfaceC(ix^
d,idim1)
7479 end subroutine mhd_update_faces_hll
7482 subroutine get_resistive_electric_field(ixI^L,ixO^L,wp,sCT,s,jce)
7487 integer,
intent(in) :: ixi^
l, ixo^
l
7489 double precision,
intent(in) :: wp(ixi^s,1:nw)
7490 type(state),
intent(in) :: sct, s
7492 double precision :: jce(ixi^s,
sdim:3)
7495 double precision :: jcc(ixi^s,7-2*
ndir:3)
7497 double precision :: xs(ixgs^t,1:
ndim)
7499 double precision :: eta(ixi^s)
7500 double precision :: gradi(ixgs^t)
7501 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
7503 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
7509 if (
lvc(idim1,idim2,idir)==0) cycle
7511 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7512 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
7515 xs(ixb^s,:)=x(ixb^s,:)
7516 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
7517 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi,2)
7518 if (
lvc(idim1,idim2,idir)==1)
then
7519 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
7521 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
7528 jce(ixi^s,:)=jce(ixi^s,:)*
mhd_eta
7536 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7537 jcc(ixc^s,idir)=0.d0
7539 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7540 ixamin^
d=ixcmin^
d+ix^
d;
7541 ixamax^
d=ixcmax^
d+ix^
d;
7542 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
7544 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
7545 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
7550 end subroutine get_resistive_electric_field
7553 subroutine get_ambipolar_electric_field(ixI^L,ixO^L,w,x,fE)
7556 integer,
intent(in) :: ixi^
l, ixo^
l
7557 double precision,
intent(in) :: w(ixi^s,1:nw)
7558 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7559 double precision,
intent(out) :: fe(ixi^s,
sdim:3)
7561 double precision :: jxbxb(ixi^s,1:3)
7562 integer :: idir,ixa^
l,ixc^
l,ix^
d
7565 call mhd_get_jxbxb(w,x,ixi^
l,ixa^
l,jxbxb)
7572 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
7575 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
7576 ixamin^
d=ixcmin^
d+ix^
d;
7577 ixamax^
d=ixcmax^
d+ix^
d;
7578 fe(ixc^s,idir)=fe(ixc^s,idir)+jxbxb(ixa^s,idir)
7580 fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0
7583 end subroutine get_ambipolar_electric_field
7589 integer,
intent(in) :: ixo^
l
7599 do ix^db=ixomin^db,ixomax^db\}
7601 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7602 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
7603 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7604 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
7605 s%w(ix^
d,b3_)=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
7606 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
7609 s%w(ix^
d,b1_)=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
7610 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
7611 s%w(ix^
d,b2_)=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
7612 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
7655 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
7656 double precision,
intent(inout) :: ws(ixis^s,1:nws)
7657 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7659 double precision :: adummy(ixis^s,1:3)
7665 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
7668 integer,
intent(in) :: ixi^
l, ixo^
l
7669 double precision,
intent(in) :: w(ixi^s,1:nw)
7670 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7671 double precision,
intent(out):: rfactor(ixi^s)
7673 double precision :: iz_h(ixo^s),iz_he(ixo^s)
7677 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)
7679 end subroutine rfactor_from_temperature_ionization
7681 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
7683 integer,
intent(in) :: ixi^
l, ixo^
l
7684 double precision,
intent(in) :: w(ixi^s,1:nw)
7685 double precision,
intent(in) :: x(ixi^s,1:
ndim)
7686 double precision,
intent(out):: rfactor(ixi^s)
7690 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.